Univariate Forecasting

63.1 Univariate Forecasting: Why Bother? (20 min)

A univariate model forecasts a series using only its own history. Why start there when we have a whole warehouse of other indicators?

1. It's the benchmark. Any multivariate model (Day 68) must beat the best univariate one, or the extra data isn't helping.
2. It's robust. It needs no other series to be published on time, so it can't be broken by another source's outage or revision.
3. It's often surprisingly good. Inflation is persistent, and its own past carries much of the predictable signal at short horizons. Forecasting research has repeatedly found simple univariate benchmarks hard to beat (e.g. Atkeson & Ohanian, 2001; Stock & Watson, 2007, Journal of Money, Credit and Banking).

Choosing the transformation

Our target, yoy CPI inflation, is seasonally neutral. Features based on mom changes carry seasonal patterns, which is why month_of_year is a feature (Day 61), letting trees learn month-specific behaviour.

63.2 Statistical Benchmarks: AR Models with statsmodels (25 min)

An autoregression AR(p) forecasts the series from its own last p values with a linear equation. It's the workhorse univariate benchmark.

statsmodels fits AR models and selects the lag order by an information criterion:

from statsmodels.tsa.ar_model import AutoReg, ar_select_order

y = load_cpi()["cpi_yoy"]
sel = ar_select_order(y.loc[:"2011-12"], maxlag=13, ic="aic", old_names=False)
print("AIC-selected lags:", sel.ar_lags)

res = AutoReg(y.loc[:"2011-12"], lags=sel.ar_lags, old_names=False).fit()
print(res.summary().tables[1])
print(res.forecast(steps=12).round(2))

AR models forecast recursively (each step feeds the next). To compare them fairly with our direct ML models, run them through a rolling-origin loop on the same origins and horizons. Add to forecasting/evaluate.py:

from statsmodels.tsa.ar_model import AutoReg, ar_select_order


def rolling_origin_ar(y: pd.Series, horizons, test_start: str, maxlag: int = 13,
                      model_name: str = "ar_aic", reselect_every: int = 12) -> pd.DataFrame:
    """AR(p) refitted at every origin (expanding window); p re-selected by AIC periodically."""
    rows, lags = [], None
    origins = y.index[y.index >= pd.Timestamp(test_start)]
    H = max(horizons)
    for i, t in enumerate(origins):
        hist = y.loc[:t]
        if lags is None or i % reselect_every == 0:
            lags = ar_select_order(hist, maxlag=maxlag, ic="aic", old_names=False).ar_lags or [1]
        path = AutoReg(hist, lags=lags, old_names=False).fit().forecast(steps=H)
        for h in horizons:
            t_target = t + pd.DateOffset(months=h)
            if t_target > y.index.max():
                continue
            rows.append({"model": model_name, "h": h, "origin": t, "target_date": t_target,
                         "forecast": float(path.iloc[h - 1]), "actual": float(y.loc[t_target])})
    return pd.DataFrame(rows)

Run it and add it to your Day 62 results:

ar = rolling_origin_ar(y, H, "2012-01-01")
results = pd.concat([pd.read_parquet("data/forecasts/day62_results.parquet"), ar])
print(score(results).round(3).to_string(index=False))

How does AR(AIC) compare with XGBoost at each horizon?

Other classic univariate models (ARIMA, exponential smoothing) are available in statsmodels.tsa and fit the same rolling pattern. An optional exercise: add ExponentialSmoothing or an ARIMA(1,0,1) and see whether either changes the ranking.

63.3 Direct vs Recursive Multi-Step Forecasts (23 min)

Day 61 introduced both strategies. Let's compare them with the same model class (XGBoost) and pure lag features, which recursive forecasting needs: every feature must itself be forecastable from predictions.

from forecasting.features import lag_features
from forecasting.models import Naive, xgb_factory

LAGS = 12
XL = lag_features(y, LAGS, "yoy")               # yoy_lag0 … yoy_lag11 only


def recursive_forecast(model, last_values: list[float], steps: int) -> list[float]:
    """Iterate a 1-step model: each prediction becomes yoy_lag0 for the next step."""
    window, preds = list(last_values), []          # window[0] = most recent value
    for _ in range(steps):
        x = pd.DataFrame([window[:LAGS]], columns=[f"yoy_lag{k}" for k in range(LAGS)])
        p = float(model.predict(x)[0])
        preds.append(p)
        window = [p] + window[:-1]
    return preds


rows = []
origins = y.index[y.index >= "2012-01-01"]
target1 = y.shift(-1)
for i, t in enumerate(origins):
    if i % 3:                                       # refit every 3 months for speed
        continue
    train = XL.notna().all(axis=1) & target1.notna() & (XL.index.shift(1, freq="MS") <= t)
    m = xgb_factory().fit(XL[train], target1[train])
    path = recursive_forecast(m, XL.loc[t].tolist(), max(H))
    for h in H:
        t_target = t + pd.DateOffset(months=h)
        if t_target <= y.index.max():
            rows.append({"model": "xgb_recursive", "h": h, "origin": t,
                         "target_date": t_target, "forecast": path[h - 1],
                         "actual": float(y.loc[t_target])})
recursive = pd.DataFrame(rows)
direct = rolling_origin(XL, y, xgb_factory, H, "2012-01-01", "xgb_direct", refit_every=3)
naive = rolling_origin(XL, y, Naive, H, "2012-01-01", "naive", refit_every=3)
print(score(pd.concat([recursive, direct, naive])).round(3).to_string(index=False))

Typical findings, which you should verify on your data:

• At h = 1 the two strategies are identical by construction (both use a 1-step model).
• At longer horizons, recursive XGBoost forecasts often flatten or drift: errors compound, and the tree's limited range (62.5) means iterated predictions converge towards typical training values.
• Direct models avoid compounding, but each horizon has slightly fewer usable training rows (61.4).

63.4 Is the Difference Real? The Diebold–Mariano Test (25 min)

Model A has RMSE 0.95 and model B 1.00 at h = 6. Is A really better, or is that luck over ~150 test months? The Diebold–Mariano (DM) test (Diebold & Mariano, 1995, Journal of Business & Economic Statistics) asks whether the average difference in losses between two forecasts is zero.

1. For each test period, compute the loss differential d = (error_A)² − (error_B)².
2. Test whether the mean of d is zero, using a HAC (Newey–West) standard error. Forecast errors at horizon h overlap and are serially correlated, so ordinary standard errors would be too small.

This is a regression of d on a constant with HAC standard errors, a tool you know from econometrics:

import statsmodels.api as sm


def diebold_mariano(results: pd.DataFrame, model_a: str, model_b: str, h: int) -> dict:
    """DM test of equal squared-error accuracy; negative stat favours model_a."""
    a = results[(results.model == model_a) & (results.h == h)].set_index("origin")
    b = results[(results.model == model_b) & (results.h == h)].set_index("origin")
    common = a.index.intersection(b.index)
    ea = a.loc[common, "forecast"] - a.loc[common, "actual"]
    eb = b.loc[common, "forecast"] - b.loc[common, "actual"]
    d = (ea ** 2 - eb ** 2).to_numpy()
    fit = sm.OLS(d, np.ones(len(d))).fit(cov_type="HAC", cov_kwds={"maxlags": max(h - 1, 0)})
    return {"model_a": model_a, "model_b": model_b, "h": h, "n": len(d),
            "mean_loss_diff": float(d.mean()), "dm_stat": float(fit.tvalues[0]),
            "p_value": float(fit.pvalues[0])}

Add it to forecasting/evaluate.py and test the best model at each horizon against the naive benchmark:

for h in H:
    print(diebold_mariano(results, "xgboost", "naive", h))

Caveats worth knowing:

• With few test periods, DM can over-reject. Harvey, Leybourne and Newbold (1997) proposed a small-sample correction, which you may add as a refinement.
• Running many pairwise tests invites false discoveries. Pre-specify the comparisons that matter (best model vs naive, vs AR) rather than testing everything against everything.

The humbling lesson: you'll often find that the "winning" model's advantage over the naive or AR benchmark is not statistically significant, especially at long horizons. Reporting that honestly is a mark of a good forecaster, not a failure.

63.5 Prediction Intervals: Quantile Models and Coverage (23 min)

A point forecast of "2.4%" hides the uncertainty. Policymakers want a range: "we expect inflation between 1.6% and 3.3% with 80% probability."

Method 1: empirical intervals from past errors

Use the backtest errors themselves. For each horizon, take the 10th and 90th percentiles of past forecast errors and add them to the point forecast. It's simple, model-agnostic and easy to explain.

Method 2: quantile regression with XGBoost

Instead of minimising squared error (which targets the mean), minimise the quantile (pinball) loss, which targets a chosen quantile. Recent XGBoost versions support this directly:

def xgb_quantile_factory(alpha: float):
    return lambda: XGBRegressor(objective="reg:quantileerror", quantile_alpha=alpha,
                                n_estimators=300, learning_rate=0.03, max_depth=3,
                                subsample=0.8, colsample_bytree=0.8, random_state=0)

H6 = [6]
q10 = rolling_origin(X, y, xgb_quantile_factory(0.10), H6, "2012-01-01", "q10", refit_every=3)
q90 = rolling_origin(X, y, xgb_quantile_factory(0.90), H6, "2012-01-01", "q90", refit_every=3)

(If your XGBoost version rejects reg:quantileerror, use scikit-learn's GradientBoostingRegressor(loss="quantile", alpha=...), which plugs into the same factory pattern.)

Checking calibration: coverage

An 80% interval should contain the actual outcome about 80% of the time. That's testable:

iv = (q10.set_index("origin")[["forecast", "actual"]].rename(columns={"forecast": "lo"})
        .join(q90.set_index("origin")[["forecast"]].rename(columns={"forecast": "hi"})))
inside = (iv["actual"] >= iv["lo"]) & (iv["actual"] <= iv["hi"])
print(f"80% interval coverage at h=6: {inside.mean():.0%}  (n={len(iv)})")
print(f"average width: {(iv['hi'] - iv['lo']).mean():.2f} pp")
print("coverage by period:")
print(inside.groupby(iv.index.year).mean().round(2))

Look at coverage by year. Intervals calibrated on calm years typically under-cover in 2021–22: the model's sense of uncertainty came from a period of stable inflation. That's an important, honest finding to report.

63.6 Capstone: A Univariate Leaderboard for CPI Inflation (20 min)

Produce a reproducible leaderboard script, forecasting/univariate.py, that:

1. Loads CPI from the mart (load_cpi, start 1992).
2. Runs the rolling-origin backtest from 2012-01 for horizons 1, 3, 6, 12 for: naive, historical mean, 2% target, AR(AIC), Ridge, random forest, XGBoost (direct) and XGBoost (recursive).
3. Computes MAE, RMSE and relative RMSE on common origins.
4. Runs Diebold–Mariano tests of every model against the naive benchmark at each horizon.
5. Computes 80% empirical-interval coverage for the best model at h = 6 and h = 12.
6. Writes data/forecasts/univariate_leaderboard.parquet and a short Markdown report, data/forecasts/univariate_report.md, that records the package versions, sample, test period and refit frequency.

"""Univariate CPI inflation leaderboard (Day 63)."""
import platform
from datetime import datetime, timezone
from pathlib import Path

import pandas as pd
import sklearn
import statsmodels
import xgboost
from sklearn.linear_model import RidgeCV
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from forecasting.data import load_cpi
from forecasting.evaluate import diebold_mariano, rolling_origin, rolling_origin_ar, score
from forecasting.features import univariate_cpi_features
from forecasting.models import Constant, HistoricalMean, Naive, rf_factory, xgb_factory

OUT = Path("data/forecasts")
H, TEST_START, REFIT = [1, 3, 6, 12], "2012-01-01", 3


def main():
    OUT.mkdir(parents=True, exist_ok=True)
    df = load_cpi()
    X, y = univariate_cpi_features(df), df["cpi_yoy"]

    factories = {
        "naive": Naive, "hist_mean": HistoricalMean, "target_2pct": Constant,
        "ridge": lambda: make_pipeline(StandardScaler(), RidgeCV(alphas=[0.1, 1, 10, 100])),
        "random_forest": rf_factory, "xgboost": xgb_factory,
    }
    results = pd.concat(
        [rolling_origin(X, y, f, H, TEST_START, n, refit_every=REFIT) for n, f in factories.items()]
        + [rolling_origin_ar(y, H, TEST_START)], ignore_index=True)

    board = score(results)
    dm = pd.DataFrame([diebold_mariano(results, m, "naive", h)
                       for m in results.model.unique() if m != "naive" for h in H])
    board = board.merge(dm.rename(columns={"model_a": "model"})[["model", "h", "dm_stat", "p_value"]],
                        on=["model", "h"], how="left")
    board.to_parquet(OUT / "univariate_leaderboard.parquet", index=False)
    results.to_parquet(OUT / "univariate_forecasts.parquet", index=False)

    report = [
        "# Univariate CPI inflation leaderboard",
        f"Generated {datetime.now(timezone.utc):%Y-%m-%d %H:%M UTC} · Python {platform.python_version()} · "
        f"xgboost {xgboost.__version__} · scikit-learn {sklearn.__version__} · "
        f"statsmodels {statsmodels.__version__}",
        f"Sample from {df.index.min():%Y-%m} · test origins from {TEST_START[:7]} · "
        f"refit every {REFIT} months · direct forecasts (AR: recursive)",
        "", board.round(3).to_markdown(index=False),
    ]
    (OUT / "univariate_report.md").write_text("\n".join(report))
    print("\n".join(report))


if __name__ == "__main__":
    main()

python -m forecasting.univariate

Interpreting your leaderboard

Write five sentences in notebooks/day63.ipynb answering:

1. Which model is best at each horizon, by relative RMSE?
2. Which advantages over the naive benchmark are statistically significant (DM)?
3. How did all models perform during the 2021–22 surge, and why?
4. Are the 80% intervals well calibrated? In which years do they fail?
5. What would you tell a central bank economist about univariate ML forecasts of Canadian inflation?

Recap: Days 61–63

Looking ahead: Days 64–67

• Hyperparameter tuning without leakage: time-series cross-validation inside the training window, then a single honest test.
• SHAP values: explaining why XGBoost forecasts what it does, globally and for each forecast (e.g. "which lags drove the 2022 forecast?").
• Zero-shot forecasting with Chronos: a pretrained time-series foundation model that forecasts CPI without being trained on it, evaluated in the same harness against today's leaderboard.

Day 63 checklist

• I can explain why univariate models are the essential benchmark, and choose stationary, seasonally appropriate targets
• AR(AIC) is evaluated in the same rolling-origin framework as the ML models
• I compared direct and recursive multi-step strategies and can explain the differences
• I test forecast-accuracy differences with Diebold–Mariano (HAC), and report non-significance honestly
• I produce prediction intervals and check their coverage over time
• A reproducible leaderboard script records results, tests, coverage and environment details

git add forecasting notebooks && git commit -m "Day 63: AR benchmarks, direct vs recursive, DM tests, intervals, leaderboard" && git push