diff --git a/events/20260917/.gitignore b/events/20260917/.gitignore new file mode 100644 index 0000000000000000000000000000000000000000..5acb80fc7ba6a56306f4347091cd85a9720cc256 --- /dev/null +++ b/events/20260917/.gitignore @@ -0,0 +1,6 @@ +data/raw/* +data/processed/* +!data/**/.gitkeep +output/* +!output/.gitkeep +__pycache__/ diff --git a/events/20260917/README.md b/events/20260917/README.md new file mode 100644 index 0000000000000000000000000000000000000000..bb060abccbfeeccad4b5b47f7e8c6240d45e0de9 --- /dev/null +++ b/events/20260917/README.md @@ -0,0 +1,360 @@ +--- +title: "TimesFM 3.0 on German agricultural and forestry price indices" +date: 2026-09-17 +issue: https://aidaho-edu.uni-hohenheim.de/gitlab/modellierung/modellierungs-fitnessclub/-/work_items/2 +--- + +# 1. Objective + +Test whether the zero-shot multivariate foundation model **TimesFM 3.0** +(`google/timesfm-3.0-pytorch`, released August 2026) forecasts monthly German +producer price indices for agricultural and forestry products better than +classical econometric models that are re-estimated on the same information set. + +"Better" is decided on an out-of-sample rolling-origin backtest, not on a single +train/test cut, because a single 12-month test window on four series gives an +effective sample size close to one. + +Three questions from the work item, made operational: + +| Question | Operationalisation | +|---|---| +| Which models to compare? | Random walk, seasonal naive, drift, damped-trend ETS, SARIMA, VAR (multivariate) vs. TimesFM 3.0 univariate and TimesFM 3.0 multivariate | +| Which metrics? | MASE, MAE, RMSE, sMAPE for the point forecast; pinball loss and 80 % interval coverage for the predictive distribution; Diebold-Mariano test against the random walk | +| Non-stationarity, seasonality, regime shifts? | Section 5: log transform + differencing for the econometric models, raw levels for TimesFM (it normalises internally); seasonality is explicit in SARIMA/ETS, implicit in TimesFM; regime shifts get their own named-origin stress set | + +A secondary objective is methodological: build a backtest harness where any new +model is one function, so the club can re-use it on other datasets in later +sessions. + +# 2. Input data + +## 2.1 Primary source + +Destatis GENESIS-Online, price indices in agriculture and forestry: + + +| Table | Content | Frequency | Base | +|---|---|---|---| +| `61231` | Producer price index for forestry products — oak, beech, spruce, pine stemwood | monthly | 2020 = 100 | +| `61211` | Producer price index for agricultural products — wheat, barley, sugar beet, potatoes | monthly | 2020 = 100 | + +Both are index series, so they are unit-free and directly comparable across +species and crops — convenient for a multivariate model that shares one +normalisation across channels. + +## 2.2 Access + +GENESIS closed its anonymous `GAST` account. Verified on 2026-09-17: the REST +API at `https://www-genesis.destatis.de/genesisWS/rest/2020/` answers +`helloworld/whoami` without credentials, but `catalogue/tables` and +`data/tablefile` return `401 Code 15`. Registration is free. + +`src/get_data.py` therefore supports three input paths: + +1. `--source genesis` — REST API, credentials from `GENESIS_USERNAME` / + `GENESIS_PASSWORD` (or `GENESIS_TOKEN`). +2. `--source csv --csv-file ` — a flat CSV (`ffcsv`) downloaded by hand + from the GENESIS web interface. **No account needed.** This is the fallback + for the session if nobody has registered. +3. `--source synthetic` — a generated stand-in panel with the same shape, + seasonality and two deliberate regime shifts. Lets the whole pipeline run + offline, and is the right thing to develop the harness against. + +## 2.3 Structure after ingestion + +`data/processed/prices_monthly.csv`, tidy long format: + +| column | type | meaning | +|---|---|---| +| `date` | ISO date, first of month | observation month | +| `series_id` | string | e.g. `spruce_stemwood`, `wheat` | +| `group` | string | `forestry` or `agriculture` | +| `value` | float | price index, base 2020 = 100 | + +Expected coverage: roughly 2005-01 to the most recent published month +(2026-07 or 2026-08 given the one-to-two month publication lag), i.e. about +250 monthly observations per series. The loader prints the coverage it actually +found; every date in the design below is derived from the last observation `T`, +not hard-coded. + +# 3. Output data + +Everything lands in `output/`. + +| File | Schema | +|---|---| +| `forecasts.csv` | `origin, series_id, model, horizon, target_date, y_true, y_pred, q10 … q90` — one row per (origin, series, model, horizon) | +| `metrics_by_horizon.csv` | `model, series_id, horizon, n, mae, rmse, smape, mase, pinball, coverage_80` | +| `metrics_overall.csv` | `model, series_id, split` aggregated over horizons, plus `dm_stat`, `dm_pvalue` against the random walk | +| `fanchart_.png` | last origin, median + 80 % band, TimesFM vs SARIMA vs truth | +| `mase_by_horizon.png` | MASE against horizon, one line per model | +| `summary.md` | the table the group reads out loud in the wrap-up phase | + +`forecasts.csv` is the atomic artefact: metrics and plots are pure functions of +it, so a metric can be added after the fact without re-running the model. + +## 3.1 Visualising the forecasts against the test set + +`src/plots.py` reads `forecasts.csv` and replots without touching a model. Note +that with rolling origins there is **no single out-of-sample path to draw**: +each series carries 36 overlapping 12-month trajectories. That is why there are +several figures rather than one, each answering a different question. + +| File | Question it answers | +|---|---| +| `fixed_horizon_h.png` | The literal one. Fix h, and each point comes from a different origin, so it collapses to a genuine out-of-sample series: what each model said h months ahead against what happened. One panel per series. | +| `trajectories_.png` | All 36 forecast paths laid over the realised series, one panel per model. Shows *how* a model fails - lagging a turning point looks like a comb of paths peeling off in the wrong direction, which no summary metric conveys. | +| `scatter_pred_actual.png` | Predicted against realised, coloured by horizon. Systematic bias is a cloud off the diagonal; a model that only fails far ahead shows as colour separation. | +| `calibration.png` | Reliability over the 9 quantiles. A curve sitting on one side of the diagonal is bias; a curve *crossing* it is the wrong interval width. | +| `coverage_by_horizon.png` | Where the 80 % interval breaks down as the horizon grows. | +| `error_distribution.png` | Boxplot of absolute scaled error. The mean MASE in the summary hides its own tail, and for a price forecast the worst month is often what matters. | + +Read `fixed_horizon` and `trajectories` together: the first says how big the +errors are, the second says what kind of error it is. + +# 4. Training set and test set + +TimesFM is **zero-shot** — it is never fitted on these data. So "training set" +means two different things and the split must be fair to both: + +- for TimesFM: the **context window** handed to the model at an origin; +- for the econometric models: the **estimation sample**, re-fitted from scratch + at every origin. + +Both see exactly the same information: all observations up to and including the +origin month, nothing after. This is enforced in one place, `iter_origins()` in +`src/backtest.py`. + +## 4.1 The split, anchored on the last observation `T` + +With `T = 2026-07` (adjust automatically if the loader finds a later month): + +| Split | Origins | Count | Purpose | +|---|---|---|---| +| **Context / estimation pool** | 2005-01 … origin | expanding | never evaluated, only consumed | +| **Development** | `T-95` … `T-48` → 2018-08 … 2022-07 | 48 | choose context length, log vs. level, covariate set, SARIMA order | +| **Test (hold-out)** | `T-47` … `T-12` → 2022-08 … 2025-07 | 36 | reported once, after development is frozen | + +Horizon `h = 1 … 12` at every origin. The test window stops at `T-12` so every +origin has a full 12 months of realised values — no ragged edge, no horizon +where the model count differs. + +36 origins × 12 horizons × 8 series ≈ 3 456 forecast points per model. Enough +for a Diebold-Mariano test at short horizons; at `h = 12` the overlapping +windows make the errors strongly autocorrelated, so the DM test uses a +Newey-West correction with lag `h-1`. + +Rolling origins are **expanding**, not sliding: a 2005 observation still helps a +2025 SARIMA. For TimesFM the context is truncated to the last `--context-length` +months (default 512, well inside the model's 15 360-step limit). + +## 4.2 Two extra evaluation sets + +**Regime-shift stress set.** Six named origins, each reported on its own rather +than averaged away: + +| Origin | Event | +|---|---| +| 2018-06 | onset of the spruce bark-beetle calamity — supply shock, price collapse | +| 2020-03 | COVID-19 | +| 2021-09 | construction-timber price surge | +| 2022-02 | Russian invasion of Ukraine — grain and energy | +| 2022-08 | peak of the energy-price spike | +| 2023-06 | disinflation / mean reversion | + +A foundation model that has seen millions of series may extrapolate a shock +better — or may smooth it into a trend. Either result is interesting, and both +are invisible in an average over 36 origins. + +**Leakage-aware honest set.** TimesFM 3.0 was released in August 2026 and its +pretraining corpus is not documented at series level. Destatis indices are +public and widely mirrored, so origins before ~2026 may be **contaminated**: the +model may have memorised the realised path. Origins from `T-11` onward +(2026-01 …) with `h = 1 … 6` are almost certainly after the pretraining cutoff. +Few origins, weak statistics, but it is the only honest zero-shot comparison and +it must be reported next to the main result, not instead of it. + +This is the single biggest threat to the validity of the whole exercise and +should be stated in the wrap-up. + +## 4.3 What counts as a win + +TimesFM beats a baseline if its MASE is lower at the same horizon **and** the DM +test rejects equal predictive accuracy at 5 %. On monthly price indices the +random walk is famously hard to beat, so "TimesFM ties the random walk at +`h = 1` and beats seasonal naive at `h = 12`" is a plausible, publishable and +entirely acceptable outcome. + +# 5. Non-stationarity, seasonality, regime shifts + +| Issue | Econometric models | TimesFM 3.0 | +|---|---|---| +| Non-stationarity | log transform, then `d=1`; unit-root behaviour is assumed rather than tested per origin (an ADF test at every origin would be a pre-test bias) | raw levels; the model applies reversible instance normalisation and optional linear detrending internally (`use_linear_detrending=True`) | +| Seasonality | explicit: `D=1, s=12` in SARIMA, additive seasonal in ETS, and the seasonal naive baseline | implicit; no frequency indicator since 2.5, the model infers periodicity from the context | +| Regime shifts | none of the baselines model breaks; that is the point — they are the honest status quo | unknown behaviour, which is what the stress set measures | + +Forecasts from log-space models are converted back with `exp()` on the median, +which is the median of the level — appropriate here, since MASE and pinball loss +are both median-oriented. No smearing correction is applied; if the group wants +mean forecasts this has to change and should be flagged. + +# 6. Models + +| Key | Type | Specification | +|---|---|---| +| `rw` | baseline | last observed value, carried forward | +| `snaive` | baseline | value 12 months earlier | +| `drift` | baseline | random walk with drift estimated over the whole estimation sample | +| `ets` | econometric | `statsmodels` Holt-Winters, damped additive trend, additive seasonal, on logs | +| `sarima` | econometric | `statsmodels` SARIMAX(1,1,1)(0,1,1)₁₂ on logs | +| `var` | econometric, multivariate | VAR on log first differences, lag order by AIC (max 12) — the fair multivariate comparator | +| `timesfm_uni` | foundation | one call per series, levels, 9 quantiles | +| `timesfm_mv` | foundation | one call per group, all four channels jointly, levels, 9 quantiles | + +`timesfm_uni` vs `timesfm_mv` isolates what the multivariate machinery of 3.0 +actually buys on these data — the headline feature of the release. + +Prediction intervals: TimesFM returns quantiles natively. SARIMA and ETS get +theirs from `get_forecast().conf_int()`. The naive baselines get an empirical +interval from the in-sample residual quantiles at each horizon, so that pinball +loss and coverage are defined for every model. + +# 7. Environment + +`environment.yml` builds a conda environment `timesfm-fitness` on Python 3.11 +(PyTorch wheels are most reliable there). TimesFM itself is pip-only: + +```bash +conda env create -f environment.yml +conda activate timesfm-fitness +python -c "import timesfm3, torch; print(torch.__version__)" +``` + +CPU is enough. The 3.0 checkpoint is ~330 M parameters; a 512-step context with +horizon 12 takes well under a second per series on a laptop, and the whole +backtest is about 300 model calls. + +**Licence warning.** The TimesFM *source code* is Apache-2.0, but the **3.0 +pretrained weights are released under `timesfm-non-commercial-license-v1.0`** +and are restricted to non-commercial, non-production use. A university workshop +is fine. Anything that ends up in a contract-funded deliverable is not. Weights +up to 2.5 remain Apache-2.0, which is the fallback if this becomes a problem. + +# 8. How to run + +```bash +conda activate timesfm-fitness +cd events/20260917 + +# 1. Data. Pick one: +python src/get_data.py --source synthetic # works offline, right now +python src/get_data.py --source csv --csv-file data/raw/61231.csv +GENESIS_USERNAME=... GENESIS_PASSWORD=... python src/get_data.py --source genesis + +# 2. Backtest. Start without TimesFM to check the harness (seconds): +python src/backtest.py --models rw,snaive,drift,ets,sarima,var --split test + +# 3. Add the foundation model (downloads ~1.3 GB on first run): +python src/backtest.py --models all --split test +python src/backtest.py --models all --split stress +python src/backtest.py --models all --split honest --horizon 6 +``` + +```bash +# 4. Replot from the saved forecasts - no model is re-run, takes seconds: +python src/plots.py --split test +python src/plots.py --split test --models timesfm_mv,sarima,rw --horizons 1,3,12 +python src/plots.py --split stress --series spruce_stemwood,wheat +``` + +`--split dev` during the specification phase, `--split test` once and once only. + +## 8.1 Measured runtime + +The whole pipeline was run end to end on the synthetic panel (259 months, +8 series, 36 origins, h = 1…12) on CPU, Python 3.11, single-threaded BLAS: + +| Step | Time | +|---|---| +| `rw`, `snaive`, `drift` | 0.6 s each, both groups | +| `ets` | ~20 s per group | +| `sarima` | ~90 s per group (warm-started from the previous origin) | +| `var` | ~1.5 s per group | +| TimesFM checkpoint load | ~30 s, once | +| `timesfm_uni` + `timesfm_mv` | ~2 min per group | +| **full `--models all --split test`** | **~8 min** | +| `--models baselines --split test` | ~4 min | + +So a full re-run fits inside the implementation hour, and the baseline-only run +is fast enough to iterate on. SARIMA dominates the cost: it is 288 maximum- +likelihood fits per group. It is warm-started from the previous origin's +estimates, which roughly triples its speed, and the warm start is discarded +whenever a fit fails to converge so a bad optimum cannot propagate down the +backtest. + +Two behaviours already visible on the synthetic panel, both worth watching for +on the real data: + +- **ETS is over-confident in the wrong direction**: 80 % coverage came out at + 0.92 against SARIMA's 0.88. That is the `sqrt(h)` interval approximation in + `forecast_ets` being too generous, exactly as flagged in the code comment. +- **Every model's intervals collapse on the stress set**: coverage fell to + 0.37-0.63 across the six regime-shift origins, against a 0.80 target. If that + reproduces on the real series it is the most decision-relevant result of the + session — a forecast interval that fails precisely when prices move is worse + than no interval at all. + +# 9. Open questions for the specification phase + +1. **Data access** — does anyone have a GENESIS account, or do we download the + two CSVs by hand at the start of the session? +2. **Scope** — forestry (61231) only, agriculture (61211) only, or both? Both is + eight series and a richer multivariate test; one table is faster. +3. **Covariates** — TimesFM 3.0 supports past-only and past-and-future + covariates. Obvious candidates: energy price index, diesel, fertiliser, + HICP, harvest volumes. Worth the extra ingestion effort today, or a follow-up + session? +4. **Horizon** — is 12 months the decision-relevant horizon for the people in + the room, or is 1-3 months (contract negotiation) what actually matters? +5. **Leakage** — do we accept the contaminated main result with a caveat, or + make the post-cutoff honest set the headline? +6. **Higher frequency** — the work item welcomes higher-frequency data. Weekly + timber auction or grain spot prices would make the foundation model much more + interesting, but needs a source. + +# 10. Status + +The pipeline runs end to end on the synthetic panel: all three splits +(`test`, `stress`, `honest`) produce forecasts, metrics, plots and a summary. +The test split yields exactly 8 models x 8 series x 36 origins x 12 horizons = +27 648 forecast rows. + +The TimesFM 3.0 calls are verified by execution, not only by reading the +package: `TimesFM3Forecaster.from_pretrained("google/timesfm-3.0-pytorch", +device="cpu")` loads in ~30 s, univariate `predict_batch` returns +`forecast (12,)` and `quantiles (12, 9)`, multivariate returns +`forecast (k, 12)` and `quantiles (k, 12, 9)` — the shapes the harness assumes. + +**Not yet verified:** `parse_genesis_csv` has never seen a real GENESIS export. +The column sniffing is deliberately tolerant and prints what it found, but +budget five minutes at the start of the session to adjust it. + +**Meaningless so far:** every number in `output/` comes from the synthetic +panel. It demonstrates that the harness works, and nothing about German price +formation. + +# 11. Files + +``` +events/20260917/ +├── README.md this document +├── environment.yml conda environment +├── src/ +│ ├── get_data.py GENESIS API / CSV / synthetic → tidy panel +│ ├── metrics.py MASE, pinball, coverage, Diebold-Mariano +│ ├── backtest.py rolling-origin harness, all models +│ └── plots.py replots forecasts.csv, no model re-run +├── data/ gitignored, except this structure +└── output/ gitignored +``` diff --git a/events/20260917/data/processed/.gitkeep b/events/20260917/data/processed/.gitkeep new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/events/20260917/data/raw/.gitkeep b/events/20260917/data/raw/.gitkeep new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/events/20260917/environment.yml b/events/20260917/environment.yml new file mode 100644 index 0000000000000000000000000000000000000000..7f333183c7df136ff43b54575f9500267d0c5b61 --- /dev/null +++ b/events/20260917/environment.yml @@ -0,0 +1,53 @@ +# Modellierungs-Fitnessclub, session 2026-09-17 +# TimesFM 3.0 on German agricultural and forestry producer price indices +# +# conda env create -f environment.yml +# conda activate timesfm-fitness +# +# Python 3.11: the version with the least friction for PyTorch wheels. +name: timesfm-fitness + +channels: + - conda-forge + +dependencies: + - python=3.11 + + # data handling + - numpy>=1.26.4 + - pandas>=2.2 + - pyarrow # parquet, and a faster csv reader + - requests # GENESIS REST API + + # econometric baselines + - statsmodels>=0.14 # SARIMAX, ExponentialSmoothing, VAR + - scipy>=1.11 + - scikit-learn>=1.4 + + # plots + - matplotlib>=3.8 + + # notebooks, optional but handy in a workshop + - jupyterlab + - ipykernel + + - pip>=24 + - pip: + # TimesFM is not on conda-forge. The [torch] extra pulls a CPU-capable + # torch wheel; on a machine with CUDA, install torch from the PyTorch + # index first and pip will keep it. + - timesfm[torch]==3.0.2 + + # checkpoint download, pinned by timesfm but listed for clarity + - huggingface_hub>=0.28.0 + - safetensors>=0.5.3 + +# Checkpoint (~1.3 GB, cached in ~/.cache/huggingface after the first run): +# google/timesfm-3.0-pytorch +# +# Licence: the TimesFM source is Apache-2.0, but the 3.0 *weights* are under +# timesfm-non-commercial-license-v1.0 — non-commercial, non-production use only. +# Weights up to 2.5 are Apache-2.0 if that ever becomes a constraint. +# +# Smoke test after creating the environment: +# python -c "import torch, statsmodels, timesfm3; print(torch.__version__)" diff --git a/events/20260917/output/.gitkeep b/events/20260917/output/.gitkeep new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/events/20260917/src/backtest.py b/events/20260917/src/backtest.py new file mode 100644 index 0000000000000000000000000000000000000000..922e6ae5a4d21271e0f4510e0ba69afe05360338 --- /dev/null +++ b/events/20260917/src/backtest.py @@ -0,0 +1,655 @@ +"""Rolling-origin backtest: TimesFM 3.0 against econometric baselines. + +Every model sees exactly the same information set at every origin: all +observations up to and including the origin month, nothing after. That rule is +enforced in one place, `iter_origins` + the `history` slice in `run_split`, so +there is a single line to audit for look-ahead bias. + + python src/backtest.py --models rw,snaive,drift,ets,sarima,var --split test + python src/backtest.py --models all --split test + python src/backtest.py --models all --split stress + python src/backtest.py --models all --split honest + +Adding a model is one function with the signature + + fit_predict(history, horizon, warm) -> (point[H], quantiles[H, 9], warm_out) + +registered in UNIVARIATE_MODELS. `warm` carries the previous origin's fitted +parameters so that iterative models can warm-start; return None if unused. +""" + +from __future__ import annotations + +import argparse +import contextlib +import os +import sys +import time +import warnings +from pathlib import Path + +# Hundreds of small SARIMAX fits oversubscribe the BLAS thread pool and end up +# slower than single-threaded. Must be set before numpy is imported. +for _v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS", + "NUMEXPR_NUM_THREADS", "VECLIB_MAXIMUM_THREADS"): + os.environ.setdefault(_v, "1") + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +from scipy import stats + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +from metrics import ( # noqa: E402 + QUANTILE_COLS, QUANTILE_LEVELS, metrics_by_horizon, metrics_overall, + seasonal_mae_scale, +) + +warnings.filterwarnings("ignore") # statsmodels convergence noise, by the hundred + + +@contextlib.contextmanager +def quiet(): + """Silence warnings inside a fit. + + A module-level filterwarnings is not enough: statsmodels re-arms the warning + filter inside its own fit routines, so hundreds of ConvergenceWarnings end up + on stdout and bury the progress output. Non-convergence is handled explicitly + below by refusing to warm-start from a failed fit, so suppressing the message + does not suppress the problem. + """ + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + yield + +HERE = Path(__file__).resolve().parent.parent +PANEL = HERE / "data" / "processed" / "prices_monthly.csv" +OUTPUT = HERE / "output" + +SEASON = 12 +Z = np.array([stats.norm.ppf(q) for q in QUANTILE_LEVELS]) + +# Named origins for the regime-shift stress set. Each is reported on its own: +# averaging a bark-beetle collapse together with 35 quiet months hides exactly +# the behaviour we want to see. +STRESS_ORIGINS = { + "2018-06-01": "bark-beetle calamity onset", + "2020-03-01": "COVID-19", + "2021-09-01": "construction timber surge", + "2022-02-01": "Russian invasion of Ukraine", + "2022-08-01": "energy price peak", + "2023-06-01": "disinflation", +} + + +# --------------------------------------------------------------------------- # +# Splits +# --------------------------------------------------------------------------- # +def iter_origins(dates: pd.DatetimeIndex, split: str, horizon: int) -> list: + """Forecast origins for a split. + + Every date is derived from the last observation T, so the design survives a + data update without anyone editing a constant. + + dev T-95 .. T-48 48 origins, for tuning + test T-47 .. T-12 36 origins, reported once + stress six named origins, reported individually + honest T-11 .. T-h after the TimesFM 3.0 pretraining cutoff + """ + dates = pd.DatetimeIndex(sorted(dates)) + n = len(dates) + if split == "dev": + lo, hi = n - 96, n - 48 + elif split == "test": + lo, hi = n - 48, n - horizon + elif split == "honest": + lo, hi = n - 12, n - horizon + elif split == "stress": + wanted = pd.DatetimeIndex([pd.Timestamp(d) for d in STRESS_ORIGINS]) + return [d for d in wanted if d in set(dates) and d <= dates[-1] - pd.DateOffset(months=horizon)] + else: + raise ValueError(f"unknown split {split!r}") + + lo = max(lo, SEASON * 4) # need enough history to fit a seasonal model + if hi <= lo: + hint = (" The honest split only has the last 12 months of origins, so it " + "needs a short horizon: --split honest --horizon 6." + if split == "honest" else "") + raise SystemExit( + f"Split {split!r} is empty: {n} observations is too short for " + f"horizon {horizon}.{hint}" + ) + return list(dates[lo:hi]) + + +# --------------------------------------------------------------------------- # +# Baselines: point forecast plus empirical quantiles from in-sample errors +# --------------------------------------------------------------------------- # +def _empirical_quantiles(point: np.ndarray, errors_by_h: list) -> np.ndarray: + """Quantiles from the in-sample distribution of h-step forecast errors. + + Gives the naive baselines an honest predictive distribution, so pinball loss + and interval coverage are defined for them too — without that, TimesFM would + win the distributional comparison by walkover. + """ + out = np.full((len(point), len(QUANTILE_LEVELS)), np.nan) + for i, errs in enumerate(errors_by_h): + errs = np.asarray(errs, float) + errs = errs[np.isfinite(errs)] + if len(errs) >= 10: + out[i] = point[i] + np.quantile(errs, QUANTILE_LEVELS) + elif len(errs) >= 2: # too few for empirical quantiles, fall back to normal + out[i] = point[i] + Z * errs.std(ddof=1) + return out + + +def forecast_rw(history: pd.Series, horizon: int, warm=None): + y = history.to_numpy(float) + point = np.repeat(y[-1], horizon) + errors = [y[h:] - y[:-h] if len(y) > h else np.array([]) for h in range(1, horizon + 1)] + return point, _empirical_quantiles(point, errors), None + + +def forecast_snaive(history: pd.Series, horizon: int, warm=None): + y = history.to_numpy(float) + if len(y) < SEASON: + return forecast_rw(history, horizon) + # Reference value: the same month of the year, stepping back in whole years + # until we land inside the history. Index -1 is the origin month itself. + point = np.array( + [y[-1 + h - SEASON * int(np.ceil(h / SEASON))] for h in range(1, horizon + 1)], + float, + ) + errors = [] + for h in range(1, horizon + 1): + lag = SEASON * int(np.ceil(h / SEASON)) + errors.append(y[lag:] - y[:-lag] if len(y) > lag else np.array([])) + return point, _empirical_quantiles(point, errors), None + + +def forecast_drift(history: pd.Series, horizon: int, warm=None): + y = history.to_numpy(float) + slope = (y[-1] - y[0]) / (len(y) - 1) if len(y) > 1 else 0.0 + point = y[-1] + slope * np.arange(1, horizon + 1) + errors = [] + for h in range(1, horizon + 1): + errors.append(y[h:] - (y[:-h] + slope * h) if len(y) > h else np.array([])) + return point, _empirical_quantiles(point, errors), None + + +# --------------------------------------------------------------------------- # +# Econometric models. Fitted on logs: prices are positive and multiplicative, +# and a monotone transform maps quantiles to quantiles, so exp() of a log-space +# quantile is the level quantile. No smearing correction — these are median +# forecasts, which is what MASE and pinball loss ask for. +# --------------------------------------------------------------------------- # +def forecast_sarima(history: pd.Series, horizon: int, warm=None): + from statsmodels.tsa.statespace.sarimax import SARIMAX + + ylog = np.log(history.to_numpy(float)) + model = SARIMAX( + ylog, order=(1, 1, 1), seasonal_order=(0, 1, 1, SEASON), + trend="n", enforce_stationarity=False, enforce_invertibility=False, + ) + # Consecutive origins differ by one observation, so the previous origin's + # estimates are an excellent starting point. Same specification, same + # likelihood, just fewer optimiser iterations — roughly a 3x speedup over a + # 36-origin backtest, with identical optima up to optimiser tolerance. + with quiet(): + try: + res = (model.fit(disp=False, start_params=warm) if warm is not None + else model.fit(disp=False)) + except Exception: + res = model.fit(disp=False) + + fc = res.get_forecast(steps=horizon) + mean, se = np.asarray(fc.predicted_mean, float), np.asarray(fc.se_mean, float) + quant = np.exp(mean[:, None] + Z[None, :] * se[:, None]) + + # Never warm-start the next origin from a fit that did not converge: that + # would carry a bad optimum forward through the whole backtest instead of + # letting the next origin retry from the default starting values. + converged = bool(getattr(res, "mle_retvals", {}).get("converged", False)) + return np.exp(mean), quant, (np.asarray(res.params, float) if converged else None) + + +def forecast_ets(history: pd.Series, horizon: int, warm=None): + from statsmodels.tsa.holtwinters import ExponentialSmoothing + + ylog = np.log(history.to_numpy(float)) + seasonal = "add" if len(ylog) >= 2 * SEASON + 4 else None + with quiet(): + res = ExponentialSmoothing( + ylog, trend="add", damped_trend=True, seasonal=seasonal, + seasonal_periods=SEASON if seasonal else None, + initialization_method="estimated", + ).fit(optimized=True) + mean = np.asarray(res.forecast(horizon), float) + # Holt-Winters has no closed-form predictive variance here. Approximate it by + # the residual sd growing with sqrt(h) — a random-walk-ish assumption that is + # generous at long horizons. Flagged rather than hidden. + sigma = float(np.std(res.resid[np.isfinite(res.resid)], ddof=1)) + se = sigma * np.sqrt(np.arange(1, horizon + 1)) + quant = np.exp(mean[:, None] + Z[None, :] * se[:, None]) + return np.exp(mean), quant, None + + +UNIVARIATE_MODELS = { + "rw": forecast_rw, + "snaive": forecast_snaive, + "drift": forecast_drift, + "ets": forecast_ets, + "sarima": forecast_sarima, +} + + +# --------------------------------------------------------------------------- # +# VAR: the econometric multivariate comparator for timesfm_mv +# --------------------------------------------------------------------------- # +def forecast_var(history: pd.DataFrame, horizon: int) -> dict: + """VAR on log first differences, lag order by AIC. Returns {series_id: (point, quant)}.""" + from statsmodels.tsa.api import VAR + + logs = np.log(history.to_numpy(float)) + dlog = np.diff(logs, axis=0) + k = dlog.shape[1] + maxlags = min(12, max(1, (len(dlog) - k - 2) // (k + 1))) + with quiet(): + res = VAR(dlog).fit(maxlags=maxlags, ic="aic") + + steps_dlog = res.forecast(dlog[-res.k_ar:], horizon) # (H, k) + log_point = logs[-1] + np.cumsum(steps_dlog, axis=0) + + # Forecast error variance of the *cumulated* differences. res.mse(h) is the + # covariance of the h-step error on the differences; summing the diagonals + # ignores the cross-horizon correlation, so this understates the true level + # uncertainty. Adequate for ranking coverage, not for a risk statement. + mse = res.mse(horizon) # (H, k, k) + var_cum = np.cumsum(np.array([np.diag(m) for m in mse]), axis=0) + se = np.sqrt(var_cum) # (H, k) + + out = {} + for j, sid in enumerate(history.columns): + mean = log_point[:, j] + quant = np.exp(mean[:, None] + Z[None, :] * se[:, j][:, None]) + out[sid] = (np.exp(mean), quant) + return out + + +# --------------------------------------------------------------------------- # +# TimesFM 3.0 +# --------------------------------------------------------------------------- # +_FORECASTER = None + + +def get_forecaster(device: str, checkpoint: str, batch_size: int): + """Load the checkpoint once and keep it. ~1.3 GB on the first call.""" + global _FORECASTER + if _FORECASTER is None: + try: + from timesfm3 import TimesFM3Forecaster + except ImportError: + raise SystemExit( + "timesfm is not installed. `pip install 'timesfm[torch]==3.0.2'`, " + "or run with --models rw,snaive,drift,ets,sarima,var to check the " + "harness without it." + ) + print(f"Loading {checkpoint} on {device} ...") + _FORECASTER = TimesFM3Forecaster.from_pretrained( + checkpoint, device=device, per_core_batch_size=batch_size + ) + return _FORECASTER + + +def run_timesfm_univariate(wide: pd.DataFrame, origins, horizon: int, + context_length: int, forecaster) -> list: + """One context per (origin, series). All of them go through in one batch.""" + contexts, keys = [], [] + for origin in origins: + hist = wide.loc[:origin] + for sid in wide.columns: + y = hist[sid].dropna().to_numpy(np.float32) + if len(y) < SEASON * 2: + continue + contexts.append(y[-context_length:]) + keys.append((origin, sid)) + + print(f" timesfm_uni: {len(contexts)} series-origins in one batch") + outputs = list(forecaster.predict_batch( + contexts=contexts, horizon=horizon, return_quantiles=True, + )) + return [(k, np.asarray(o.forecast, float), np.asarray(o.quantiles, float)) + for k, o in zip(keys, outputs)] + + +def run_timesfm_multivariate(wide: pd.DataFrame, origins, horizon: int, + context_length: int, forecaster) -> list: + """All channels of a group forecast jointly — the headline feature of 3.0. + + Context shape is (num_variates, context_length); the output is (k, H) and + (k, H, 9). Comparing this against timesfm_uni isolates what the multivariate + attention actually buys on these series. + """ + contexts, keys = [], [] + for origin in origins: + hist = wide.loc[:origin].dropna() + if len(hist) < SEASON * 2: + continue + block = hist.to_numpy(np.float32).T[:, -context_length:] # (k, L) + contexts.append(block) + keys.append(origin) + + print(f" timesfm_mv: {len(contexts)} origins x {wide.shape[1]} variates") + outputs = list(forecaster.predict_batch( + contexts=contexts, horizon=horizon, return_quantiles=True, + )) + + results = [] + for origin, out in zip(keys, outputs): + point = np.asarray(out.forecast, float) # (k, H) + quant = np.asarray(out.quantiles, float) # (k, H, 9) + for j, sid in enumerate(wide.columns): + results.append(((origin, sid), point[j], quant[j])) + return results + + +# --------------------------------------------------------------------------- # +# The harness +# --------------------------------------------------------------------------- # +def collect_rows(key, model, point, quant, wide, horizon) -> list: + """Attach realised values to a forecast and flatten to tidy rows.""" + origin, sid = key + rows = [] + for h in range(1, horizon + 1): + target = origin + pd.DateOffset(months=h) + if target not in wide.index: + continue + row = { + "origin": origin, "series_id": sid, "model": model, "horizon": h, + "target_date": target, "y_true": float(wide.at[target, sid]), + "y_pred": float(point[h - 1]), + } + if quant is not None and np.ndim(quant) == 2 and quant.shape[0] >= h: + row.update(dict(zip(QUANTILE_COLS, quant[h - 1].astype(float)))) + rows.append(row) + return rows + + +def run_split(panel: pd.DataFrame, models: list, split: str, horizon: int, + context_length: int, device: str, checkpoint: str, + batch_size: int): + all_rows, scales = [], {} + + for group, gdf in panel.groupby("group"): + wide = gdf.pivot(index="date", columns="series_id", values="value").sort_index() + wide.index = pd.DatetimeIndex(wide.index) + origins = iter_origins(wide.index, split, horizon) + if not origins: + print(f"[{group}] no usable origins for split {split!r}, skipping") + continue + print(f"\n[{group}] {len(wide)} months, {wide.shape[1]} series, " + f"{len(origins)} origins: {origins[0]:%Y-%m} .. {origins[-1]:%Y-%m}") + + # MASE scale, per origin, from the estimation sample only. + for origin in origins: + for sid in wide.columns: + hist = wide.loc[:origin, sid].dropna().to_numpy(float) + scales[(origin, sid)] = seasonal_mae_scale(hist, SEASON) + + # ---- univariate baselines and econometric models -------------------- # + for name in [m for m in models if m in UNIVARIATE_MODELS]: + t0 = time.time() + print(f" {name:<8}", end="", flush=True) + warm_cache = {} # (name, sid) -> previous origin's fitted parameters + for i, origin in enumerate(origins): + for sid in wide.columns: + hist = wide.loc[:origin, sid].dropna() + if len(hist) < SEASON * 3: + continue + try: + point, quant, warm = UNIVARIATE_MODELS[name]( + hist, horizon, warm_cache.get(sid) + ) + except Exception as exc: # a single bad fit must not kill the run + print(f"\n {name} failed at {origin:%Y-%m} {sid}: {exc}") + continue + if warm is not None: + warm_cache[sid] = warm + all_rows += collect_rows((origin, sid), name, point, quant, + wide, horizon) + if (i + 1) % 6 == 0: + print(".", end="", flush=True) + print(f" {time.time() - t0:5.1f}s") + + # ---- VAR ------------------------------------------------------------ # + if "var" in models and wide.shape[1] > 1: + t0 = time.time() + print(f" {'var':<8}", end="", flush=True) + for origin in origins: + hist = wide.loc[:origin].dropna() + if len(hist) < SEASON * 4: + continue + try: + per_series = forecast_var(hist, horizon) + except Exception as exc: + print(f"\n var failed at {origin:%Y-%m}: {exc}") + continue + for sid, (point, quant) in per_series.items(): + all_rows += collect_rows((origin, sid), "var", point, quant, + wide, horizon) + print(f" {time.time() - t0:5.1f}s") + + # ---- TimesFM -------------------------------------------------------- # + tfm_wanted = [m for m in models if m.startswith("timesfm")] + if tfm_wanted: + forecaster = get_forecaster(device, checkpoint, batch_size) + if "timesfm_uni" in tfm_wanted: + for key, point, quant in run_timesfm_univariate( + wide, origins, horizon, context_length, forecaster + ): + all_rows += collect_rows(key, "timesfm_uni", point, quant, + wide, horizon) + if "timesfm_mv" in tfm_wanted and wide.shape[1] > 1: + for key, point, quant in run_timesfm_multivariate( + wide, origins, horizon, context_length, forecaster + ): + all_rows += collect_rows(key, "timesfm_mv", point, quant, + wide, horizon) + + return pd.DataFrame(all_rows), scales + + +# --------------------------------------------------------------------------- # +# Reporting +# --------------------------------------------------------------------------- # +def plot_mase_by_horizon(by_h: pd.DataFrame, path: Path) -> None: + agg = by_h.groupby(["model", "horizon"])["mase"].mean().reset_index() + fig, ax = plt.subplots(figsize=(8, 5)) + for model, g in agg.groupby("model"): + style = "-o" if model.startswith("timesfm") else "--" + ax.plot(g.horizon, g.mase, style, label=model, linewidth=2 if style == "-o" else 1.4) + ax.axhline(1.0, color="grey", linewidth=0.8) + ax.annotate("seasonal naive, in-sample", (0.02, 1.02), xycoords=("axes fraction", "data"), + fontsize=8, color="grey") + ax.set_xlabel("forecast horizon (months)") + ax.set_ylabel("MASE, averaged over series and origins") + ax.set_title("Forecast accuracy by horizon") + ax.legend(fontsize=8) + fig.tight_layout() + fig.savefig(path, dpi=140) + plt.close(fig) + + +def plot_fanchart(forecasts: pd.DataFrame, panel: pd.DataFrame, sid: str, + path: Path, models=("timesfm_uni", "sarima")) -> None: + g = forecasts[forecasts.series_id == sid] + if g.empty: + return + origin = g.origin.max() + hist = (panel[panel.series_id == sid].set_index("date")["value"] + .sort_index().loc[:origin].tail(60)) + + fig, ax = plt.subplots(figsize=(10, 5)) + ax.plot(hist.index, hist.values, color="black", linewidth=1.2, label="observed") + + truth = g[(g.origin == origin) & (g.model == g.model.iloc[0])] + ax.plot(truth.target_date, truth.y_true, color="black", linewidth=1.2, linestyle=":", + label="realised") + + colours = {"timesfm_uni": "tab:red", "timesfm_mv": "tab:orange", "sarima": "tab:blue", + "ets": "tab:green", "rw": "tab:grey", "var": "tab:purple"} + for model in models: + m = g[(g.origin == origin) & (g.model == model)] + if m.empty: + continue + colour = colours.get(model, None) + ax.plot(m.target_date, m.y_pred, color=colour, linewidth=2, label=model) + if {"q10", "q90"} <= set(m.columns) and m.q10.notna().any(): + ax.fill_between(m.target_date, m.q10, m.q90, color=colour, alpha=0.18, + label=f"{model} 80 %") + + ax.axvline(origin, color="grey", linestyle="--", linewidth=0.8) + ax.set_title(f"{sid} — forecast from origin {origin:%Y-%m}") + ax.set_ylabel("price index (2020 = 100)") + ax.legend(fontsize=8) + fig.tight_layout() + fig.savefig(path, dpi=140) + plt.close(fig) + + +def write_summary(by_h, overall, split, horizon, path: Path) -> None: + ranked = (overall.groupby("model")[["mase", "rmse", "pinball", "coverage_80"]] + .mean().sort_values("mase")) + lines = [ + f"# Backtest summary — split `{split}`, horizon 1..{horizon}", + "", + "## Ranking, averaged over series, origins and horizons", + "", + "| model | MASE | RMSE | pinball | 80 % coverage |", + "|---|---:|---:|---:|---:|", + ] + for model, r in ranked.iterrows(): + cov = "" if not np.isfinite(r.coverage_80) else f"{r.coverage_80:.2f}" + pin = "" if not np.isfinite(r.pinball) else f"{r.pinball:.3f}" + lines.append(f"| {model} | {r.mase:.3f} | {r.rmse:.2f} | {pin} | {cov} |") + + lines += ["", "MASE below 1 beats the in-sample seasonal naive. Coverage should", + "sit near 0.80 — a model that is accurate but overconfident is not", + "usable for decisions, and this column is where that shows up.", ""] + + h_max = int(by_h.horizon.max()) + lines += [f"## MASE at h=1 and h={h_max}", "", + f"| model | h=1 | h={h_max} |", "|---|---:|---:|"] + piv = by_h.groupby(["model", "horizon"])["mase"].mean().unstack() + for model in ranked.index: + a = piv.loc[model, 1] if 1 in piv.columns else np.nan + b = piv.loc[model, h_max] if h_max in piv.columns else np.nan + lines.append(f"| {model} | {a:.3f} | {b:.3f} |") + + dm_cols = [c for c in overall.columns if c.startswith("dm_pvalue")] + if dm_cols: + lines += ["", "## Diebold-Mariano against the random walk", "", + "Negative statistic: the model beats the random walk. " + "Origins overlap, so p-values are optimistic — read them as a " + "rough screen, not as proof.", "", + "| model | series | " + " | ".join( + c.replace("dm_pvalue_", "p ") for c in dm_cols) + " |", + "|---|---|" + "---:|" * len(dm_cols)] + for _, r in overall[overall.model != "rw"].iterrows(): + ps = " | ".join( + "" if not np.isfinite(r[c]) else f"{r[c]:.3f}" for c in dm_cols) + lines.append(f"| {r.model} | {r.series_id} | {ps} |") + + if split == "stress": + lines += ["", "## Regime-shift origins", "", + "| origin | event |", "|---|---|"] + for d, label in STRESS_ORIGINS.items(): + lines.append(f"| {pd.Timestamp(d):%Y-%m} | {label} |") + + if split != "honest": + lines += ["", "## Caveat", "", + "TimesFM 3.0 was released in August 2026 and its pretraining " + "corpus is undocumented at series level. These Destatis indices " + "are public, so origins before 2026 may be contaminated. Run " + "`--split honest` for the post-cutoff comparison and report both."] + + path.write_text("\n".join(lines) + "\n", encoding="utf-8") + + +# --------------------------------------------------------------------------- # +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--panel", type=Path, default=PANEL) + ap.add_argument("--models", default="rw,snaive,drift,ets,sarima,var", + help="comma separated, or 'all', or 'baselines'") + ap.add_argument("--split", choices=["dev", "test", "stress", "honest"], default="test") + ap.add_argument("--horizon", type=int, default=12) + ap.add_argument("--context-length", type=int, default=512, + help="months of context handed to TimesFM (max 15360)") + ap.add_argument("--device", default="cpu", help="cpu or cuda") + ap.add_argument("--checkpoint", default="google/timesfm-3.0-pytorch") + ap.add_argument("--batch-size", type=int, default=32) + ap.add_argument("--group", choices=["forestry", "agriculture"], default=None) + ap.add_argument("--out", type=Path, default=OUTPUT) + args = ap.parse_args() + + if not args.panel.exists(): + sys.exit(f"No panel at {args.panel}. Run: python src/get_data.py --source synthetic") + + all_models = list(UNIVARIATE_MODELS) + ["var", "timesfm_uni", "timesfm_mv"] + if args.models == "all": + models = all_models + elif args.models == "baselines": + models = list(UNIVARIATE_MODELS) + ["var"] + else: + models = [m.strip() for m in args.models.split(",") if m.strip()] + unknown = [m for m in models if m not in all_models] + if unknown: + sys.exit(f"Unknown model(s) {unknown}. Known: {all_models}") + + panel = pd.read_csv(args.panel, parse_dates=["date"]) + if args.group: + panel = panel[panel.group == args.group] + print(f"Panel: {panel.date.min():%Y-%m} .. {panel.date.max():%Y-%m}, " + f"{panel.series_id.nunique()} series") + print(f"Split: {args.split}, horizon 1..{args.horizon}, models: {models}") + + forecasts, scales = run_split( + panel, models, args.split, args.horizon, args.context_length, + args.device, args.checkpoint, args.batch_size, + ) + if forecasts.empty: + sys.exit("No forecasts were produced — check the panel length and the split.") + + args.out.mkdir(parents=True, exist_ok=True) + suffix = f"_{args.split}" + forecasts.to_csv(args.out / f"forecasts{suffix}.csv", index=False) + + by_h = metrics_by_horizon(forecasts, scales) + overall = metrics_overall(forecasts, scales, reference_model="rw") + by_h.to_csv(args.out / f"metrics_by_horizon{suffix}.csv", index=False) + overall.to_csv(args.out / f"metrics_overall{suffix}.csv", index=False) + + plot_mase_by_horizon(by_h, args.out / f"mase_by_horizon{suffix}.png") + plotted = [m for m in ("timesfm_uni", "timesfm_mv", "sarima", "rw") if m in models] + for sid in forecasts.series_id.unique(): + plot_fanchart(forecasts, panel, sid, + args.out / f"fanchart_{sid}{suffix}.png", models=plotted) + + write_summary(by_h, overall, args.split, args.horizon, + args.out / f"summary{suffix}.md") + + print("\n" + "=" * 66) + print(f"Ranking on split {args.split!r} (lower MASE is better)") + print("=" * 66) + ranked = (by_h.groupby("model")[["mase", "rmse", "pinball", "coverage_80"]] + .mean().sort_values("mase")) + print(ranked.to_string(float_format=lambda v: f"{v:8.3f}")) + print(f"\nWrote forecasts, metrics, plots and summary{suffix}.md to {args.out}") + + +if __name__ == "__main__": + main() diff --git a/events/20260917/src/get_data.py b/events/20260917/src/get_data.py new file mode 100644 index 0000000000000000000000000000000000000000..e5aa60f3832166ef02b1d93f3e39ff5cb1983629 --- /dev/null +++ b/events/20260917/src/get_data.py @@ -0,0 +1,321 @@ +"""Build the monthly price-index panel for the TimesFM 3.0 experiment. + +Three input paths, because GENESIS closed its anonymous account: + + --source synthetic generated stand-in panel, no network, no account + --source csv a flat CSV downloaded by hand from GENESIS + --source genesis REST API, needs GENESIS_USERNAME / GENESIS_PASSWORD + +All three write the same tidy long table to data/processed/prices_monthly.csv: + + date,series_id,group,value + +Usage: + python src/get_data.py --source synthetic + python src/get_data.py --source csv --csv-file data/raw/61231.csv + GENESIS_USERNAME=... GENESIS_PASSWORD=... python src/get_data.py --source genesis +""" + +from __future__ import annotations + +import argparse +import io +import os +import sys +from pathlib import Path + +import numpy as np +import pandas as pd + +HERE = Path(__file__).resolve().parent.parent +RAW = HERE / "data" / "raw" +PROCESSED = HERE / "data" / "processed" + +GENESIS_BASE = "https://www-genesis.destatis.de/genesisWS/rest/2020" + +# Destatis GENESIS tables named in the work item. +TABLES = { + "61231-0002": "forestry", # producer price index, forestry products + "61211-0002": "agriculture", # producer price index, agricultural products +} + +# The series we want out of those tables, and the label fragments that identify +# them. GENESIS labels are verbose, so we match on a fragment, case-insensitively. +SERIES = { + "forestry": { + "oak_stemwood": ["oak", "eiche"], + "beech_stemwood": ["beech", "buche"], + "spruce_stemwood": ["spruce", "fichte"], + "pine_stemwood": ["pine", "kiefer"], + }, + "agriculture": { + "wheat": ["wheat", "weizen"], + "barley": ["barley", "gerste"], + "sugar_beet": ["sugar beet", "zuckerrüben", "zuckerrueben"], + "potatoes": ["potato", "kartoffel"], + }, +} + + +# --------------------------------------------------------------------------- # +# 1. Synthetic panel +# --------------------------------------------------------------------------- # +def make_synthetic(start="2005-01-01", end="2026-07-01", seed=20260917) -> pd.DataFrame: + """A stand-in panel with the shape, seasonality and shocks of the real data. + + Not a simulation of German price formation. Its job is to let the backtest + harness be written, debugged and demonstrated before anyone has downloaded + anything, and to give the regime-shift stress set something to bite on. + """ + rng = np.random.default_rng(seed) + dates = pd.date_range(start, end, freq="MS") + n = len(dates) + t = np.arange(n) + year = dates.year + (dates.month - 1) / 12.0 + + def window(lo, hi): + """Smooth 0->1->0 bump over the calendar window [lo, hi].""" + mid, half = (lo + hi) / 2, (hi - lo) / 2 + return np.exp(-(((year - mid) / (half / 2.0)) ** 2)) + + # Shocks shared across the panel, with different loadings per series. + energy_spike = window(2021.5, 2023.0) * 45.0 # gas, diesel, fertiliser + covid = window(2020.1, 2020.8) * -6.0 + beetle = window(2018.4, 2021.0) * -28.0 # bark-beetle calamity + common = np.cumsum(rng.normal(0.0, 0.45, n)) # slow common drift + + spec = { + # series_id: (group, level_2020, trend/yr, seas_amp, seas_phase, + # energy_load, beetle_load, sigma, ar1) + "oak_stemwood": ("forestry", 100, 1.6, 2.0, 0.1, 0.35, 0.20, 1.3, 0.55), + "beech_stemwood": ("forestry", 100, 1.2, 2.5, 0.2, 0.40, 0.45, 1.5, 0.55), + "spruce_stemwood": ("forestry", 100, 0.6, 3.0, 0.0, 0.90, 1.00, 2.6, 0.60), + "pine_stemwood": ("forestry", 100, 0.8, 2.2, 0.0, 0.70, 0.75, 2.0, 0.58), + "wheat": ("agriculture", 100, 1.1, 4.0, 0.55, 0.80, 0.0, 2.4, 0.50), + "barley": ("agriculture", 100, 1.0, 4.5, 0.55, 0.75, 0.0, 2.6, 0.50), + "sugar_beet": ("agriculture", 100, 1.3, 1.5, 0.75, 0.45, 0.0, 1.8, 0.45), + "potatoes": ("agriculture", 100, 1.4, 12.0, 0.30, 0.35, 0.0, 5.5, 0.40), + } + + frames = [] + for sid, (grp, lvl, slope, amp, phase, e_load, b_load, sigma, ar1) in spec.items(): + seasonal = amp * np.sin(2 * np.pi * ((dates.month - 1) / 12.0 + phase)) + # Second harmonic: real price seasonality is not a clean sine. + seasonal += 0.4 * amp * np.sin(4 * np.pi * ((dates.month - 1) / 12.0 + phase)) + + # AR(1) idiosyncratic noise, so first differences are not white. + eps = rng.normal(0.0, sigma, n) + noise = np.zeros(n) + for i in range(1, n): + noise[i] = ar1 * noise[i - 1] + eps[i] + + value = ( + lvl + + slope * (year - 2020.0) + + seasonal + + e_load * energy_spike + + b_load * beetle + + covid + + 0.6 * common + + noise + ) + frames.append( + pd.DataFrame({"date": dates, "series_id": sid, "group": grp, "value": value}) + ) + + out = pd.concat(frames, ignore_index=True) + out["value"] = out["value"].clip(lower=5.0).round(1) + return out + + +# --------------------------------------------------------------------------- # +# 2. GENESIS REST API +# --------------------------------------------------------------------------- # +def fetch_genesis(table: str, startyear: int = 2005) -> str: + """Download one GENESIS table as a flat CSV (ffcsv) and return the text. + + Verified 2026-09-17: the anonymous GAST account no longer works, the API + answers 401 Code 15 without credentials. Registration at + https://www-genesis.destatis.de is free. + """ + import requests + + user = os.environ.get("GENESIS_USERNAME", "") + password = os.environ.get("GENESIS_PASSWORD", "") + token = os.environ.get("GENESIS_TOKEN", "") + if not token and not (user and password): + sys.exit( + "No GENESIS credentials. Set GENESIS_USERNAME and GENESIS_PASSWORD " + "(or GENESIS_TOKEN), or use --source csv / --source synthetic." + ) + + headers = {"Authorization": f"Bearer {token}"} if token else { + "username": user, "password": password + } + params = { + "name": table, + "area": "all", + "compress": "false", + "transpose": "false", + "startyear": str(startyear), + "format": "ffcsv", + "language": "en", + } + resp = requests.post( + f"{GENESIS_BASE}/data/tablefile", headers=headers, data=params, timeout=180 + ) + resp.raise_for_status() + if resp.text.lstrip().startswith("{"): # GENESIS reports errors as JSON 200s + sys.exit(f"GENESIS error for {table}: {resp.text[:300]}") + return resp.text + + +# --------------------------------------------------------------------------- # +# 3. Parse a GENESIS flat CSV +# --------------------------------------------------------------------------- # +def parse_genesis_csv(text: str, group: str) -> pd.DataFrame: + """Parse a GENESIS ffcsv export into the tidy panel. + + GENESIS exports are semicolon separated with German decimals, and the exact + column names depend on the table and the export options. This parser sniffs + the time and value columns rather than hard-coding them, and prints what it + found so a mismatch is visible immediately instead of silently producing an + empty panel. + """ + df = pd.read_csv(io.StringIO(text), sep=";", decimal=",", dtype=str) + df.columns = [c.strip() for c in df.columns] + print(f" columns: {list(df.columns)}") + + time_col = next( + (c for c in df.columns if c.lower() in ("time", "zeit", "time_label")), None + ) + value_col = next( + (c for c in df.columns if c.lower() in ("value", "wert", "value_variable_code")), + None, + ) + # The month sits in a "1_Auspraegung_Label"-style column, or in the time column. + month_col = next( + (c for c in df.columns if "auspraegung_label" in c.lower() or "month" in c.lower()), + None, + ) + label_cols = [c for c in df.columns if "label" in c.lower() and c != time_col] + + if time_col is None or value_col is None: + sys.exit( + "Could not identify the time and value columns in this export.\n" + f"Columns seen: {list(df.columns)}\n" + "Adjust parse_genesis_csv() — this is a five-minute fix once we see " + "the real file." + ) + + month_map = {m: i + 1 for i, m in enumerate( + ["january", "february", "march", "april", "may", "june", "july", + "august", "september", "october", "november", "december"] + )} + month_map.update({m: i + 1 for i, m in enumerate( + ["januar", "februar", "märz", "april", "mai", "juni", "juli", + "august", "september", "oktober", "november", "dezember"] + )}) + + def to_date(row): + year = str(row[time_col]).strip()[:4] + raw = str(row[month_col]).strip().lower() if month_col else "" + month = month_map.get(raw) + if month is None: # "MONAT01", "M01", "01" + digits = "".join(ch for ch in raw if ch.isdigit()) + month = int(digits) if digits else None + if month is None or not year.isdigit(): + return pd.NaT + return pd.Timestamp(int(year), int(month), 1) + + df["date"] = df.apply(to_date, axis=1) + df["value"] = pd.to_numeric( + df[value_col].astype(str).str.replace(",", ".", regex=False), errors="coerce" + ) + + # Map the verbose GENESIS labels onto our short series ids. + def to_series_id(row): + haystack = " ".join(str(row[c]).lower() for c in label_cols) + for sid, fragments in SERIES[group].items(): + if any(f in haystack for f in fragments): + return sid + return None + + df["series_id"] = df.apply(to_series_id, axis=1) + df["group"] = group + + out = df.dropna(subset=["date", "value", "series_id"]) + out = out[["date", "series_id", "group", "value"]].drop_duplicates( + subset=["date", "series_id"] + ) + print(f" kept {len(out)} rows, series: {sorted(out.series_id.unique())}") + return out.reset_index(drop=True) + + +# --------------------------------------------------------------------------- # +def report(panel: pd.DataFrame) -> None: + print("\nPanel summary") + print("-" * 60) + for (grp, sid), g in panel.groupby(["group", "series_id"]): + print( + f" {grp:<12} {sid:<18} {g.date.min():%Y-%m} .. {g.date.max():%Y-%m} " + f"n={len(g):>4} mean={g.value.mean():7.1f}" + ) + print("-" * 60) + print(f" last observation T = {panel.date.max():%Y-%m}") + print(" the backtest derives every split date from T, nothing is hard-coded") + + +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--source", choices=["synthetic", "csv", "genesis"], + default="synthetic") + ap.add_argument("--csv-file", type=Path, nargs="*", default=[], + help="one or more GENESIS csv exports (--source csv)") + ap.add_argument("--group", choices=["forestry", "agriculture"], default=None, + help="group of the --csv-file(s); inferred from the name if absent") + ap.add_argument("--startyear", type=int, default=2005) + ap.add_argument("--out", type=Path, + default=PROCESSED / "prices_monthly.csv") + args = ap.parse_args() + + RAW.mkdir(parents=True, exist_ok=True) + PROCESSED.mkdir(parents=True, exist_ok=True) + + if args.source == "synthetic": + print("Generating the synthetic stand-in panel (no network, no account).") + panel = make_synthetic() + + elif args.source == "genesis": + frames = [] + for table, group in TABLES.items(): + print(f"GENESIS {table} ({group}) ...") + text = fetch_genesis(table, args.startyear) + (RAW / f"{table}.csv").write_text(text, encoding="utf-8") + frames.append(parse_genesis_csv(text, group)) + panel = pd.concat(frames, ignore_index=True) + + else: # csv + if not args.csv_file: + sys.exit("--source csv needs at least one --csv-file") + frames = [] + for path in args.csv_file: + group = args.group + if group is None: + group = "agriculture" if "61211" in path.name else "forestry" + print(f"Parsing {path} as {group} ...") + frames.append( + parse_genesis_csv(path.read_text(encoding="utf-8", errors="replace"), group) + ) + panel = pd.concat(frames, ignore_index=True) + + panel = panel.sort_values(["group", "series_id", "date"]).reset_index(drop=True) + args.out.parent.mkdir(parents=True, exist_ok=True) + panel.to_csv(args.out, index=False, date_format="%Y-%m-%d") + report(panel) + print(f"\nWrote {len(panel)} rows to {args.out}") + + +if __name__ == "__main__": + main() diff --git a/events/20260917/src/metrics.py b/events/20260917/src/metrics.py new file mode 100644 index 0000000000000000000000000000000000000000..afafdedad8de8a501acf19c5c72b3c28a3b76ce2 --- /dev/null +++ b/events/20260917/src/metrics.py @@ -0,0 +1,224 @@ +"""Forecast accuracy metrics for the rolling-origin backtest. + +Everything here is a pure function of the forecast table, so a metric can be +added after the fact without re-running a single model. + + forecasts.csv -> metrics_by_horizon.csv, metrics_overall.csv +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from scipy import stats + +QUANTILE_LEVELS = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9] +QUANTILE_COLS = [f"q{int(round(q * 100))}" for q in QUANTILE_LEVELS] + + +# --------------------------------------------------------------------------- # +# Point-forecast metrics +# --------------------------------------------------------------------------- # +def mae(y, yhat): + return float(np.mean(np.abs(np.asarray(y) - np.asarray(yhat)))) + + +def rmse(y, yhat): + return float(np.sqrt(np.mean((np.asarray(y) - np.asarray(yhat)) ** 2))) + + +def smape(y, yhat): + """Symmetric MAPE in percent, on the 0-200 convention.""" + y, yhat = np.asarray(y, float), np.asarray(yhat, float) + denom = (np.abs(y) + np.abs(yhat)) / 2.0 + with np.errstate(divide="ignore", invalid="ignore"): + ratio = np.where(denom > 0, np.abs(y - yhat) / denom, 0.0) + return float(100.0 * np.mean(ratio)) + + +def seasonal_mae_scale(history: np.ndarray, season_length: int = 12) -> float: + """In-sample seasonal-naive MAE — the denominator of MASE. + + Computed on the estimation sample available at the origin, never on the + test values, so the scale is itself an out-of-sample-safe quantity. + """ + history = np.asarray(history, float) + if len(history) <= season_length: + return float("nan") + diffs = np.abs(history[season_length:] - history[:-season_length]) + scale = float(np.mean(diffs)) + return scale if scale > 0 else float("nan") + + +def mase(y, yhat, scale: float) -> float: + if not np.isfinite(scale) or scale <= 0: + return float("nan") + return mae(y, yhat) / scale + + +# --------------------------------------------------------------------------- # +# Distributional metrics +# --------------------------------------------------------------------------- # +def pinball_loss(y, quantile_matrix, levels=QUANTILE_LEVELS) -> float: + """Mean pinball (quantile) loss averaged over levels and observations. + + quantile_matrix: (n_obs, n_levels) + """ + y = np.asarray(y, float).reshape(-1, 1) + q = np.asarray(quantile_matrix, float) + if q.ndim != 2 or q.shape[1] != len(levels): + return float("nan") + levels = np.asarray(levels, float).reshape(1, -1) + diff = y - q + loss = np.maximum(levels * diff, (levels - 1.0) * diff) + return float(np.nanmean(loss)) + + +def interval_coverage(y, lower, upper) -> float: + """Empirical coverage of a prediction interval, as a fraction.""" + y, lower, upper = (np.asarray(a, float) for a in (y, lower, upper)) + ok = np.isfinite(lower) & np.isfinite(upper) + if not ok.any(): + return float("nan") + return float(np.mean((y[ok] >= lower[ok]) & (y[ok] <= upper[ok]))) + + +# --------------------------------------------------------------------------- # +# Diebold-Mariano +# --------------------------------------------------------------------------- # +def diebold_mariano(errors_a, errors_b, horizon: int = 1, power: int = 1): + """Test equal predictive accuracy of model A and model B. + + Returns (statistic, two-sided p-value). Negative statistic: A is better. + + Rolling origins with horizon h produce overlapping forecast windows, so the + loss differential is autocorrelated up to lag h-1. The long-run variance is + therefore estimated with a Newey-West kernel truncated at h-1, and the + Harvey-Leybourne-Newbold small-sample correction is applied before the + statistic is compared to a t distribution. + + Caveat worth stating out loud in the wrap-up: the standard DM test assumes + the forecasts are not nested and the models are not estimated on the same + data in a way that makes the differential degenerate. Comparing a zero-shot + foundation model against a re-estimated SARIMA is exactly the setting DM was + built for, but reading a p-value of 0.04 on 36 overlapping origins as strong + evidence would be overclaiming. + """ + ea = np.asarray(errors_a, float) + eb = np.asarray(errors_b, float) + ok = np.isfinite(ea) & np.isfinite(eb) + ea, eb = ea[ok], eb[ok] + n = len(ea) + if n < 8: + return float("nan"), float("nan") + + d = np.abs(ea) ** power - np.abs(eb) ** power + d_bar = d.mean() + d_dev = d - d_bar + + max_lag = max(0, horizon - 1) + gamma0 = float(np.mean(d_dev ** 2)) + long_run = gamma0 + for lag in range(1, max_lag + 1): + if lag >= n: + break + gamma = float(np.mean(d_dev[lag:] * d_dev[:-lag])) + long_run += 2.0 * (1.0 - lag / (max_lag + 1.0)) * gamma # Bartlett weight + if long_run <= 0: + return float("nan"), float("nan") + + stat = d_bar / np.sqrt(long_run / n) + + # Harvey, Leybourne & Newbold (1997) small-sample correction. + h = horizon + correction = (n + 1 - 2 * h + h * (h - 1) / n) / n + if correction <= 0: + return float("nan"), float("nan") + stat *= np.sqrt(correction) + + p = 2.0 * (1.0 - stats.t.cdf(abs(stat), df=n - 1)) + return float(stat), float(p) + + +# --------------------------------------------------------------------------- # +# Aggregation over the forecast table +# --------------------------------------------------------------------------- # +def metrics_by_horizon(forecasts: pd.DataFrame, scales: dict) -> pd.DataFrame: + """One row per (model, series_id, horizon). + + `scales` maps (origin, series_id) -> in-sample seasonal-naive MAE. MASE is + computed on origin-specific scales and then averaged, which keeps every + origin on its own footing instead of letting the 2022 price level dominate. + """ + rows = [] + for (model, sid, h), g in forecasts.groupby(["model", "series_id", "horizon"]): + g = g.dropna(subset=["y_true", "y_pred"]) + if g.empty: + continue + key = list(zip(g["origin"], g["series_id"])) + sc = np.array([scales.get(k, np.nan) for k in key], float) + with np.errstate(invalid="ignore", divide="ignore"): + per_origin_mase = np.abs(g["y_true"].to_numpy() - g["y_pred"].to_numpy()) / sc + + qm = g[QUANTILE_COLS].to_numpy(float) if set(QUANTILE_COLS) <= set(g.columns) else None + rows.append({ + "model": model, + "series_id": sid, + "horizon": int(h), + "n": len(g), + "mae": mae(g["y_true"], g["y_pred"]), + "rmse": rmse(g["y_true"], g["y_pred"]), + "smape": smape(g["y_true"], g["y_pred"]), + "mase": float(np.nanmean(per_origin_mase)), + "pinball": pinball_loss(g["y_true"], qm) if qm is not None else np.nan, + "coverage_80": ( + interval_coverage(g["y_true"], g["q10"], g["q90"]) + if {"q10", "q90"} <= set(g.columns) else np.nan + ), + }) + return pd.DataFrame(rows).sort_values(["series_id", "model", "horizon"]) + + +def metrics_overall(forecasts: pd.DataFrame, scales: dict, + reference_model: str = "rw") -> pd.DataFrame: + """One row per (model, series_id), aggregated over all horizons. + + The Diebold-Mariano column tests each model against `reference_model` at + horizon 1 and at the longest horizon, which is where the two stories usually + differ: a foundation model that ties the random walk one month ahead can + still win decisively twelve months ahead. + """ + by_h = metrics_by_horizon(forecasts, scales) + agg = ( + by_h.groupby(["model", "series_id"]) + .agg(n=("n", "sum"), mae=("mae", "mean"), rmse=("rmse", "mean"), + smape=("smape", "mean"), mase=("mase", "mean"), + pinball=("pinball", "mean"), coverage_80=("coverage_80", "mean")) + .reset_index() + ) + + h_max = int(forecasts["horizon"].max()) + for h_test, suffix in ((1, "h1"), (h_max, f"h{h_max}")): + stats_rows = {} + ref = forecasts[ + (forecasts.model == reference_model) & (forecasts.horizon == h_test) + ].set_index(["origin", "series_id"]) + for (model, sid), g in forecasts[forecasts.horizon == h_test].groupby( + ["model", "series_id"] + ): + g = g.set_index(["origin", "series_id"]) + r = ref.reindex(g.index) + e_model = (g["y_true"] - g["y_pred"]).to_numpy(float) + e_ref = (r["y_true"] - r["y_pred"]).to_numpy(float) + stat, p = diebold_mariano(e_model, e_ref, horizon=h_test) + stats_rows[(model, sid)] = (stat, p) + agg[f"dm_stat_{suffix}"] = [ + stats_rows.get((m, s), (np.nan, np.nan))[0] + for m, s in zip(agg.model, agg.series_id) + ] + agg[f"dm_pvalue_{suffix}"] = [ + stats_rows.get((m, s), (np.nan, np.nan))[1] + for m, s in zip(agg.model, agg.series_id) + ] + + return agg.sort_values(["series_id", "mase"]) diff --git a/events/20260917/src/plots.py b/events/20260917/src/plots.py new file mode 100644 index 0000000000000000000000000000000000000000..7f09bfd90a8979549c9c92e0e7cf78e61d9a9281 --- /dev/null +++ b/events/20260917/src/plots.py @@ -0,0 +1,359 @@ +"""Visualise the out-of-sample forecasts against the test set. + +Reads output/forecasts_.csv and replots without re-running a single +model, which is the point of keeping the forecast table as the atomic artefact. + + python src/plots.py --split test + python src/plots.py --split test --models timesfm_mv,sarima,rw --horizons 1,6,12 + +A caveat about the phrase "out-of-sample forecast versus test set": with a +rolling-origin design there is no single out-of-sample path to draw. Each series +has 36 overlapping 12-month trajectories, one per origin. That is why there are +five different figures here rather than one — each answers a different question: + + 1. trajectories every forecast path laid over the realised series + 2. fixed_horizon the literal "prediction vs actual" line chart, at fixed h + 3. scatter predicted against actual, coloured by horizon + 4. calibration are the 9 quantiles honest? + 5. coverage where does the 80 % interval break down? +""" + +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + +sys.path.insert(0, str(Path(__file__).resolve().parent)) +from metrics import QUANTILE_COLS, QUANTILE_LEVELS, seasonal_mae_scale # noqa: E402 + +HERE = Path(__file__).resolve().parent.parent +OUTPUT = HERE / "output" +PANEL = HERE / "data" / "processed" / "prices_monthly.csv" + +# Foundation models get a solid line and a marker, everything else a dashed one, +# so the comparison of interest is readable in greyscale and on a projector. +COLOURS = { + "timesfm_uni": "tab:red", "timesfm_mv": "tab:brown", "sarima": "tab:blue", + "ets": "tab:orange", "var": "tab:purple", "rw": "tab:green", + "drift": "tab:cyan", "snaive": "tab:pink", +} + + +def style(model: str) -> dict: + if model.startswith("timesfm"): + return {"color": COLOURS.get(model), "linewidth": 2.0, "linestyle": "-"} + return {"color": COLOURS.get(model), "linewidth": 1.3, "linestyle": "--"} + + +def grid(n: int, ncols: int = 2, size=(7.0, 2.6)): + nrows = int(np.ceil(n / ncols)) + fig, axes = plt.subplots(nrows, ncols, figsize=(size[0] * ncols, size[1] * nrows), + squeeze=False) + return fig, axes.ravel() + + +# --------------------------------------------------------------------------- # +# 1. Every forecast trajectory, laid over the realised series +# --------------------------------------------------------------------------- # +def plot_trajectories(fc: pd.DataFrame, panel: pd.DataFrame, models: list, + out: Path, split: str) -> None: + """One panel per model: all origins' 12-month paths over the truth. + + This is the honest picture of a rolling-origin backtest. A model that tracks + the level but lags turning points shows up as a comb of paths each starting + on the series and peeling away in the wrong direction — something no summary + metric communicates. + """ + for sid, g in fc.groupby("series_id"): + obs = (panel[panel.series_id == sid].set_index("date")["value"].sort_index()) + lo = g.origin.min() - pd.DateOffset(months=6) + obs = obs.loc[lo:g.target_date.max()] + + fig, axes = grid(len(models), ncols=2) + for ax, model in zip(axes, models): + m = g[g.model == model] + ax.plot(obs.index, obs.values, color="black", linewidth=1.4, zorder=3) + for _, path in m.groupby("origin"): + path = path.sort_values("horizon") + ax.plot(path.target_date, path.y_pred, + color=COLOURS.get(model, "grey"), linewidth=0.8, alpha=0.55) + ax.set_title(f"{model} ({m.origin.nunique()} origins)", fontsize=10) + ax.tick_params(labelsize=8) + for ax in axes[len(models):]: + ax.axis("off") + fig.suptitle(f"{sid} — every forecast trajectory vs the realised series " + f"(split: {split})", fontsize=12) + fig.tight_layout() + fig.savefig(out / f"trajectories_{sid}_{split}.png", dpi=140) + plt.close(fig) + + +# --------------------------------------------------------------------------- # +# 2. The literal answer: prediction against actual, at a fixed horizon +# --------------------------------------------------------------------------- # +def plot_fixed_horizon(fc: pd.DataFrame, models: list, horizons: list, + out: Path, split: str) -> None: + """Fix the horizon, and the backtest collapses to one comparable line per model. + + At a fixed h every point comes from a different origin, so this is a genuine + out-of-sample series: what each model said h months ahead, against what + happened. The 80 % band is drawn for the first model only, otherwise the + shading is unreadable. + """ + series = sorted(fc.series_id.unique()) + for h in horizons: + sub = fc[fc.horizon == h] + if sub.empty: + continue + fig, axes = grid(len(series), ncols=2, size=(7.0, 2.4)) + for ax, sid in zip(axes, series): + s = sub[sub.series_id == sid] + truth = s[s.model == models[0]].sort_values("target_date") + ax.plot(truth.target_date, truth.y_true, color="black", linewidth=1.8, + label="realised", zorder=3) + for j, model in enumerate(models): + m = s[s.model == model].sort_values("target_date") + if m.empty: + continue + ax.plot(m.target_date, m.y_pred, label=model, **style(model)) + if j == 0 and {"q10", "q90"} <= set(m.columns) and m.q10.notna().any(): + ax.fill_between(m.target_date, m.q10, m.q90, alpha=0.15, + color=COLOURS.get(model, "grey"), + label=f"{model} 80 %") + ax.set_title(sid, fontsize=10) + ax.tick_params(labelsize=8) + for ax in axes[len(series):]: + ax.axis("off") + axes[0].legend(fontsize=7, ncol=2) + fig.suptitle(f"Forecast {h} month(s) ahead vs the realised test set " + f"(split: {split})", fontsize=12) + fig.tight_layout() + fig.savefig(out / f"fixed_horizon_h{h:02d}_{split}.png", dpi=140) + plt.close(fig) + + +# --------------------------------------------------------------------------- # +# 3. Predicted against actual +# --------------------------------------------------------------------------- # +def plot_scatter(fc: pd.DataFrame, models: list, out: Path, split: str) -> None: + """Points off the 45-degree line are errors; colour shows which horizon. + + Systematic bias shows up as a cloud sitting above or below the diagonal, and + a model that only fails at long horizons shows up as colour separation. + """ + fig, axes = grid(len(models), ncols=3, size=(3.6, 3.7)) + lo = float(min(fc.y_true.min(), fc.y_pred.min())) + hi = float(max(fc.y_true.max(), fc.y_pred.max())) + for ax, model in zip(axes, models): + m = fc[fc.model == model] + sc = ax.scatter(m.y_true, m.y_pred, c=m.horizon, cmap="viridis", + s=5, alpha=0.5) + ax.plot([lo, hi], [lo, hi], color="black", linewidth=1.0) + ax.set_title(model, fontsize=10) + ax.set_xlabel("realised", fontsize=8) + ax.set_ylabel("predicted", fontsize=8) + ax.tick_params(labelsize=7) + for ax in axes[len(models):]: + ax.axis("off") + # tight_layout does not apply once the colorbar steals space from the axes, + # so the row spacing is set by hand: without it the top row's x-label lands + # on the bottom row's title. + fig.subplots_adjust(hspace=0.42, top=0.88) + fig.colorbar(sc, ax=axes.tolist(), label="horizon (months)", shrink=0.7) + fig.suptitle(f"Predicted vs realised, all origins and horizons (split: {split})", + fontsize=12) + fig.savefig(out / f"scatter_pred_actual_{split}.png", dpi=140, + bbox_inches="tight") + plt.close(fig) + + +# --------------------------------------------------------------------------- # +# 4. Are the quantiles honest? +# --------------------------------------------------------------------------- # +def plot_calibration(fc: pd.DataFrame, models: list, out: Path, split: str) -> None: + """Reliability diagram over the 9 predicted quantiles. + + For each nominal level q, the fraction of realised values falling below the + predicted q-quantile. On the diagonal the distribution is calibrated. Two + different faults show up here and they read differently: + + - the whole curve sits on one side of the diagonal: the predictive + distribution is *shifted*. Above the diagonal means too many realised + values landed below the forecasts, i.e. the model forecast too high. + - the curve *crosses* the diagonal: the distribution is the wrong *width*. + Above at low levels and below at high levels means the quantiles are + squeezed towards the median, so the intervals are too narrow. The other + way round means they are too wide. + + This is where a model that wins on MASE can still be unusable for a decision. + """ + fig, ax = plt.subplots(figsize=(6.4, 5.6)) + for model in models: + m = fc[fc.model == model] + if not set(QUANTILE_COLS) <= set(m.columns) or m[QUANTILE_COLS].isna().all().all(): + continue + empirical = [float((m.y_true <= m[c]).mean()) for c in QUANTILE_COLS] + ax.plot(QUANTILE_LEVELS, empirical, marker="o", markersize=4, **style(model), + label=model) + ax.plot([0, 1], [0, 1], color="black", linewidth=1.0, zorder=0) + ax.set_xlabel("nominal quantile level") + ax.set_ylabel("empirical fraction of realised values below") + ax.set_title( + f"Quantile calibration (split: {split})\n" + "curve on one side = biased · crossing down = intervals too narrow, " + "crossing up = too wide", fontsize=10) + ax.legend(fontsize=8) + fig.tight_layout() + fig.savefig(out / f"calibration_{split}.png", dpi=140) + plt.close(fig) + + +# --------------------------------------------------------------------------- # +# 5. Where does the interval break down? +# --------------------------------------------------------------------------- # +def plot_coverage(fc: pd.DataFrame, models: list, out: Path, split: str) -> None: + """Empirical 80 % coverage against horizon. The target is the grey line.""" + fig, ax = plt.subplots(figsize=(7.5, 5.0)) + lowest = 1.0 + for model in models: + m = fc[fc.model == model] + if not {"q10", "q90"} <= set(m.columns) or m.q10.isna().all(): + continue + cov = (m.assign(hit=(m.y_true >= m.q10) & (m.y_true <= m.q90)) + .groupby("horizon")["hit"].mean()) + lowest = min(lowest, float(cov.min())) + ax.plot(cov.index, cov.values, marker="o", markersize=4, **style(model), + label=model) + ax.axhline(0.8, color="grey", linewidth=1.2) + # Right-aligned so it does not sit on top of the curves at h=1. + ax.annotate("nominal 80 %", (0.995, 0.807), xycoords=("axes fraction", "data"), + fontsize=8, color="grey", ha="right") + # Zoom to the data instead of always showing 0..1: on the test split every + # model sits above 0.75 and the interesting differences are a few points wide. + ax.set_ylim(min(0.75, lowest - 0.05), 1.02) + ax.set_xlabel("forecast horizon (months)") + ax.set_ylabel("empirical coverage of the 80 % interval") + ax.set_title(f"Interval coverage by horizon (split: {split})", fontsize=11) + ax.legend(fontsize=8) + fig.tight_layout() + fig.savefig(out / f"coverage_by_horizon_{split}.png", dpi=140) + plt.close(fig) + + +# --------------------------------------------------------------------------- # +# 6. Distribution of scaled errors +# --------------------------------------------------------------------------- # +def plot_error_distribution(fc: pd.DataFrame, panel: pd.DataFrame, models: list, + out: Path, split: str) -> None: + """Boxplot of absolute scaled error per model. + + The mean MASE in the summary table hides its own tail. A model can win on + average and still own the worst single month of the test window, which for a + price forecast is often what people actually care about. + """ + scales = {} + for sid, g in panel.groupby("series_id"): + s = g.set_index("date")["value"].sort_index() + for origin in fc.origin.unique(): + hist = s.loc[:origin].to_numpy(float) + scales[(origin, sid)] = seasonal_mae_scale(hist, 12) + + fc = fc.copy() + fc["scale"] = [scales.get((o, s), np.nan) + for o, s in zip(fc.origin, fc.series_id)] + with np.errstate(invalid="ignore", divide="ignore"): + fc["ase"] = np.abs(fc.y_true - fc.y_pred) / fc["scale"] + + data, labels = [], [] + for model in models: + v = fc.loc[fc.model == model, "ase"].to_numpy(float) + v = v[np.isfinite(v)] + if len(v): + data.append(v) + labels.append(model) + + fig, ax = plt.subplots(figsize=(8.5, 5.0)) + # Not `labels=`: matplotlib renamed it to `tick_labels` and then removed it. + # Setting the ticks directly works on every version. + bp = ax.boxplot(data, showfliers=False, patch_artist=True, + medianprops={"color": "black"}) + ax.set_xticks(range(1, len(labels) + 1)) + ax.set_xticklabels(labels) + for patch, model in zip(bp["boxes"], labels): + patch.set_facecolor(COLOURS.get(model, "grey")) + patch.set_alpha(0.55) + for i, v in enumerate(data, start=1): + ax.scatter([i], [v.mean()], marker="D", color="black", s=22, zorder=4) + ax.axhline(1.0, color="grey", linewidth=1.0) + ax.set_ylabel("absolute scaled error (1.0 = in-sample seasonal naive)") + ax.set_title(f"Error distribution, not just the mean (split: {split})\n" + "box = quartiles, diamond = mean (the MASE in the summary table)", + fontsize=11) + ax.tick_params(axis="x", rotation=20, labelsize=9) + fig.tight_layout() + fig.savefig(out / f"error_distribution_{split}.png", dpi=140) + plt.close(fig) + + +# --------------------------------------------------------------------------- # +def main() -> None: + ap = argparse.ArgumentParser(description=__doc__, + formatter_class=argparse.RawDescriptionHelpFormatter) + ap.add_argument("--split", default="test") + ap.add_argument("--forecasts", type=Path, default=None) + ap.add_argument("--panel", type=Path, default=PANEL) + ap.add_argument("--out", type=Path, default=OUTPUT) + ap.add_argument("--models", default=None, + help="comma separated; default: all models in the file") + ap.add_argument("--horizons", default="1,6,12") + ap.add_argument("--series", default=None, help="comma separated subset") + args = ap.parse_args() + + path = args.forecasts or (args.out / f"forecasts_{args.split}.csv") + if not path.exists(): + sys.exit(f"No forecasts at {path}. Run src/backtest.py --split {args.split} first.") + + fc = pd.read_csv(path, parse_dates=["origin", "target_date"]) + panel = pd.read_csv(args.panel, parse_dates=["date"]) + if args.series: + keep = [s.strip() for s in args.series.split(",")] + fc = fc[fc.series_id.isin(keep)] + panel = panel[panel.series_id.isin(keep)] + + present = list(fc.model.unique()) + models = ([m.strip() for m in args.models.split(",")] if args.models + else sorted(present, key=lambda m: (not m.startswith("timesfm"), m))) + missing = [m for m in models if m not in present] + if missing: + sys.exit(f"Not in {path.name}: {missing}. Available: {sorted(present)}") + + horizons = [int(h) for h in args.horizons.split(",") + if int(h) in set(fc.horizon.unique())] + + args.out.mkdir(parents=True, exist_ok=True) + print(f"{len(fc)} forecast rows, {fc.series_id.nunique()} series, " + f"{fc.origin.nunique()} origins, models: {models}") + + plot_fixed_horizon(fc, models, horizons, args.out, args.split) + print(f" fixed_horizon_h*_{args.split}.png {horizons}") + plot_trajectories(fc, panel, models, args.out, args.split) + print(f" trajectories__{args.split}.png {fc.series_id.nunique()} files") + plot_scatter(fc, models, args.out, args.split) + print(f" scatter_pred_actual_{args.split}.png") + plot_calibration(fc, models, args.out, args.split) + print(f" calibration_{args.split}.png") + plot_coverage(fc, models, args.out, args.split) + print(f" coverage_by_horizon_{args.split}.png") + plot_error_distribution(fc, panel, models, args.out, args.split) + print(f" error_distribution_{args.split}.png") + + +if __name__ == "__main__": + main()