Model card — Forecasting and forecast evaluation¶

Family: backtest, theta_forecast, accuracy, dm_test, cw_test, gw_test, var_backtest

Producing a forecast is easy; knowing whether it is any good is the hard part. This family covers a strong, hard-to-beat benchmark (theta_forecast), an honest pseudo-out-of-sample evaluation engine (backtest), the standard accuracy measures (accuracy), the three formal tests for comparing two forecasters (dm_test, cw_test, gw_test), and the VaR backtest battery (var_backtest) for grading tail-risk forecasts. The discipline they enforce: never report an accuracy number from a single split, and never claim one model beats another without a test that accounts for dependent forecast errors.

Function	Role
`theta_forecast`	The Theta method — a benchmark that wins forecasting competitions
`backtest`	Rolling / expanding pseudo-out-of-sample evaluation
`accuracy`	ME / RMSE / MAE / MAPE / sMAPE / MASE / RMSSE
`dm_test`	Diebold-Mariano test of equal predictive accuracy (any two forecasters)
`cw_test`	Clark-West test for nested models
`gw_test`	Giacomini-White unconditional test of equal predictive ability
`var_backtest`	Kupiec + Christoffersen + Engle-Manganelli DQ backtests of a VaR path

What it estimates¶

theta_forecast(y, steps) — the Theta method (Assimakopoulos-Nikolopoulos 2000): deseasonalize, decompose into two "theta lines," extrapolate, and recombine. It is equivalent to simple exponential smoothing with drift and is a notoriously strong benchmark.
backtest(y, ...) — walks an origin forward through the series, re-fits a chosen forecaster — a built-in name or your own Python callable (see the callable contract) — on each training window, forecasts horizon steps ahead, and tabulates accuracy by horizon. This is the correct way to estimate out-of-sample error; a single train/test split is not.
accuracy(actual, forecast) — the standard scale-dependent (RMSE, MAE), percentage (MAPE, sMAPE), and scaled (MASE, RMSSE) error measures. The scaled measures divide by a naïve in-sample benchmark, so a value below 1 means you beat that benchmark.
dm_test / cw_test / gw_test — turn "model A looks better" into a hypothesis test on the loss differential, with variance estimates that account for the autocorrelation multi-step forecast errors always carry.
var_backtest(returns, var_forecasts, alpha=0.05, dq_lags=4) — grades a Value-at-Risk forecast path on the realized hit sequence hit_t = 1{return_t < VaR_t} along two separate axes: is the violation rate right (Kupiec 1995 LR_uc ~ χ²(1)), and are the violations unclustered (Christoffersen 1998 LR_ind ~ χ²(1), jointly LR_cc = LR_uc + LR_ind ~ χ²(2)) — plus the strictest check, the Engle-Manganelli (2004) dynamic quantile regression of the demeaned hits on lagged hits and the VaR forecast itself (DQ ~ χ²(k)). Also accepts a pre-computed 0/1 violation series when you only have hits.

Assumptions¶

Diebold-Mariano tests equal unconditional expected loss. Its variance is a HAC estimate to lag h−1 (multi-step errors are MA(h−1)-correlated); the reported hln_stat further applies the Harvey-Leybourne-Newbold (1997) small-sample correction and a t-reference. DM is designed for non-nested models — applied to nested ones it is undersized and loses power.
Clark-West is the nested-model fix: under the null the larger model equals the smaller, but its extra parameters add estimation noise that inflates its MSE; CW adjusts the loss differential for exactly that noise. Use cw_test when one model is a special case of the other.
Giacomini-White (unconditional here) tests equal finite-sample predictive ability of forecasting methods, and is agnostic about nesting; it operates on two loss series you supply.
MASE/RMSSE require an in-sample series to scale against; without insample they cannot be computed. The scaling benchmark is the seasonal naïve at frequency period.
var_backtest assumes the VaR path was produced one step ahead from information available at the forecast origin, and judges it by χ² asymptotics on the hit sequence. Sign convention (fixed once): returns and VaR forecasts live on the same return scale; var_forecasts[t] is the model's α-quantile of the conditional return distribution — a negative number for small α — and a violation is return < VaR. Working in positive-loss space, negate both series first. alpha is the VaR coverage level (0.05 for a 95% VaR), not a test size.
All of these assume the forecast errors were generated by a genuine out-of-sample procedure. Feeding in in-sample residuals invalidates every number here.

When to use¶

theta_forecast as the benchmark every other model must beat, and as a fast production forecaster for seasonal series.
backtest whenever you report an accuracy figure. Expanding window for a stable process, rolling window when you suspect the data-generating process drifts.
accuracy to summarize a set of forecasts; prefer MASE/RMSSE for cross-series comparison (they are unit-free and benchmark-relative) over MAPE (undefined near zero, asymmetric).
dm_test to compare two distinct models (e.g. ARIMA vs Theta).
cw_test to compare a restricted model against the model that nests it (e.g. "add this predictor or not").
gw_test to compare two forecasting procedures including estimation uncertainty, or when you simply have two loss series.
var_backtest whenever you report a VaR number — from garch_fit variance forecasts with a distributional quantile, from quantile_regression/growth_at_risk conditional quantiles, or from an EVT tail (gpd_fit's McNeil-Frey POT VaR, shipped in 0.3.0) — before anyone trades or reports against it. The three tests answer different questions: Kupiec catches the wrong level of risk, Christoffersen/DQ catch a model that is right on average but too slow to update after a breach (violations arriving in bursts).

Key arguments and defaults¶

Call	Argument	Default	Notes
`theta_forecast`	`steps`	— (required)	forecast horizon
	`period`	`1`	seasonal period; set to 12/4 for monthly/quarterly
`backtest`	`window`	`"expanding"`	or `"rolling"`
	`train`	`20`	initial (or fixed, if rolling) training length
	`horizon`	`1`	steps forecast at each origin
	`refit_every`	`1`	re-fit cadence (speed vs freshness)
	`forecaster`	`None` (→ `naive`)	`naive`, `drift`, `mean`, `seasonal_naive`, `theta` — or any Python callable `f(train, horizon)` (see the callable-contract section)
	`period`	`None` (→ 1)	seasonal period of `seasonal_naive` / `theta` only: passing it explicitly with any other forecaster (callables included) raises (0.7.0); the MASE/RMSSE scale period is the separate `insample_period`
`accuracy`	`insample`	`None`	required for MASE / RMSSE
	`period`	`1`	seasonal frequency of the scaling benchmark
`dm_test`	`h`	`1`	forecast horizon; sets the HAC lag to h−1
	`loss`	`"squared"`	or `"absolute"`
`cw_test`	`lrv_lags`	`0`	long-run-variance lags (0 = no HAC correction)
`gw_test`	`lrv_lags`	`0`	long-run-variance lags
`var_backtest`	`alpha`	`0.05`	VaR coverage level (expected violation rate)
	`dq_lags`	`4`	lagged hits in the DQ regression (the Engle-Manganelli choice)
	`input`	`"auto"`	with `var_forecasts` the first argument is returns; without, 0/1 hits; `"hits"` combines pre-computed hits with VaR forecasts

How to read the output¶

theta_forecast → a bare array of length steps.
backtest → {"origins", "n_origins", "horizon", "forecasts", "targets", "accuracy"}. forecasts and targets are horizon × n_origins; accuracy is a list of dicts, one per horizon, each with name ("h=1", …) and the full measure set (rmse, mae, mase, rmsse, …). Read bt["accuracy"][0]["rmse"] for one-step RMSE.
accuracy → {"me", "rmse", "mae", "mape", "smape", "mase", "rmsse"}. me is mean error (bias); MASE/RMSSE below 1 beat the naïve benchmark.
dm_test → {"dm_stat", "hln_stat", "p_value", "mean_loss_diff"}. Report hln_stat and p_value. The sign of mean_loss_diff (loss₁ − loss₂) says which model won: negative ⇒ the first forecaster had lower loss.
cw_test → {"cw_stat", "p_value", "mean_adj_diff"}. A large positive cw_stat favors the larger model; the test is one-sided.
gw_test → {"gw_stat", "p_value", "df"}. gw_stat is χ²(df); small p_value ⇒ the two methods differ.
var_backtest → the counts (n_violations, expected_violations, hit_rate, the Markov transition cells n00…n11 with pi01/pi11), the three statistics with p-values (lr_uc/p_uc, lr_ind/p_ind, lr_cc/p_cc, dq_stat/p_dq with dq_df), and a plain-language verdict that names which property failed and in which direction (e.g. "Reject unconditional coverage at 5% (Kupiec LR_uc = 6.071, p = 0.014): the VaR is too conservative"). dq_df is the rank of the DQ design: a constant VaR path (an unconditional VaR) is collinear with the intercept, gets dropped (dq_var_dropped=True), and the df honestly shrinks — the verdict says so.

Failure modes¶

Single-split accuracy. One train/test cut gives a noisy, often optimistic error estimate. Use backtest and read accuracy by horizon.
DM on nested models. The classic error: comparing "model" vs "model + extra regressor" with dm_test. It is undersized there — use cw_test.
MAPE on near-zero data. MAPE explodes and becomes meaningless when actuals approach zero, and it penalizes over- and under-prediction asymmetrically. Prefer MASE / RMSSE.
Forgetting insample. accuracy cannot return MASE/RMSSE without it.
Wrong h. dm_test sets its HAC lag to h−1; passing h=1 for a 12-step forecast under-corrects the variance and overstates significance.
In-sample errors. Every test here assumes genuine out-of-sample errors; in-sample residuals make everything look significant.
VaR sign slips. Passing positive-loss VaR against raw returns makes (almost) every observation a "violation"; var_backtest errors on the all-violations case and its verdict carries a warning whenever the violation rate exceeds 50% at small α. Zero violations is also an error — a teaching one, reporting the Kupiec continuity limit −2n ln(1−α): at 250 obs and α = 5% zero violations rejects (LR = 25.6), because a VaR that is never breached is mis-calibrated too.
Too few violations. At α = 1% on 250 observations only ~2.5 violations are expected and the χ² asymptotics of all three VaR tests are unreliable (Kupiec 1995); the verdict says so whenever fewer than five violations were expected. Prefer a longer window or a larger α (exact/Monte Carlo p-values per Dufour 2006 are scoped in Module 03).

Bring your own model — Python-callable forecasters¶

backtest(forecaster=...) and the split/ACI conformal base= accept any Python callable alongside the built-in names, so a real model — statsmodels, sklearn, a hand-rolled ensemble — is evaluated by the library's own leakage-safe engine instead of by a hand-written loop that has to re-derive the alignment, refit, and no-peeking rules.

The contract. The callable is invoked as forecaster(train, horizon):

train is a read-only float64 ndarray holding only the training window for the current forecast origin t — expanding: y[0..=t]; rolling: the train most recent observations ending at t. The engine never hands the callable an observation after the origin, and every forecast is scored against targets strictly after the window that produced it. The leakage discipline is the engine's, not the callable's — but everything else that could peek at the future (scaling, hyperparameter tuning, transformation choice) belongs inside the callable, where the training slice is all it can see.
horizon is the number of steps requested: return an array-like of exactly that many finite point forecasts for steps 1..=horizon counted from the end of train (a bare scalar is rejected — return a length-1 sequence for one step). With refit_every > 1 the callable runs only at refit origins and is asked for up to refit_every - 1 + horizon steps, so its one multi-step path covers the whole block (the documented refit contract).
A Python exception raised inside the callable aborts the run and is re-raised naming the failing origin and training window, with the original exception chained as __cause__ (its type, message, and traceback survive). Wrong-length, non-coercible, or non-finite returns raise teaching errors naming the callable, the step, the origin, and the window.
The split/ACI conformal base= takes the same callable (the engines drive expanding calibration windows plus one forward call on the full sample); the result's base key then reports "<callable NAME>". EnbPI trains its own AR ensemble and refuses a callable base with a teaching error.

Performance, honestly. The Rust engine drives origins sequentially (refit blocks in order, under the GIL) for string and callable forecasters alike — there is no parallel path to lose — so a callable costs one Python call, plus your model's fit, per refit origin. The built-in string forecasters are untouched and bit-identical to the pre-callable build (pinned by a float-hex snapshot test against fixtures/backtest_string_snapshot.json). For expensive models, refit_every trades refit freshness for speed without touching the alignment.

Worked example — a statsmodels model through the engine.

import numpy as np
import tsecon
from statsmodels.tsa.ar_model import AutoReg

rng = np.random.default_rng(7)
n = 160
y = np.zeros(n)
for t in range(2, n):
    y[t] = 0.6 * y[t - 1] - 0.2 * y[t - 2] + rng.standard_normal()
y += 10.0

def ar2(train, horizon):
    # Sees ONLY the training window; returns `horizon` steps ahead.
    fit = AutoReg(np.asarray(train), lags=2, trend="c").fit()
    return fit.predict(start=len(train), end=len(train) + horizon - 1)

bt = tsecon.backtest(y, window="expanding", train=80, horizon=4,
                     forecaster=ar2, refit_every=1)
nv = tsecon.backtest(y, window="expanding", train=80, horizon=4)  # naive
print("h=1 RMSE  AR(2):", round(bt["accuracy"][0]["rmse"], 3),
      " naive:", round(nv["accuracy"][0]["rmse"], 3))

# Same origins, so the h=1 error streams feed Diebold-Mariano directly.
e_ar = np.asarray(bt["targets"][0]) - np.asarray(bt["forecasts"][0])
e_nv = np.asarray(nv["targets"][0]) - np.asarray(nv["forecasts"][0])
dm = tsecon.dm_test(e_ar, e_nv, h=1)
print("DM (HLN) p:", round(dm["p_value"], 4))

# Distribution-free intervals around the same model (split conformal).
ci = tsecon.conformal_forecast(y, horizon=4, base=ar2, alpha=0.1)
print("90%% split-conformal h=1 interval: [%.2f, %.2f]  (base: %s)"
      % (ci["lower"][0], ci["upper"][0], ci["base"]))

Expected output:

h=1 RMSE  AR(2): 0.876  naive: 1.012
DM (HLN) p: 0.0089
90% split-conformal h=1 interval: [7.66, 10.56]  (base: <callable ar2>)

Validated how. In test_backtest_callable.py: a spy callable asserts every window it is handed is exactly the documented slice (both schemes, and the refit-block walk with its enlarged step requests), ending at the origin — strictly before every target it is scored against; a perturbation test pins the same claim without trusting the spy (perturbing y[k] moves no forecast whose origin precedes k, and y[-1] — never in any training window — moves none at all); a Python reimplementation of naive is bit-identical to forecaster="naive"; and statsmodels AutoReg(lags=2, trend="c") — the same OLS-AR-with-constant iterated multi-step spec as the Rust "ar" conformal base — reproduces the "ar" results through conformal_forecast/conformal_backtest at cross-implementation tolerance (1e-6 relative on the point path; the two sides solve least squares by different routes).

Validated against¶

theta_forecast — statsmodels ThetaModel(deseasonalize=True, use_test=False), matched numerically (realgdp 8-step golden at 1e-6 relative). Note the qualifier: statsmodels' default additionally runs a seasonality pre-test (use_test=True) and skips deseasonalization when it fails, so on weakly- or non-seasonal data declared with period > 1 the two defaults diverge (measured up to a few percent on iid data with period=12); pass use_test=False on the statsmodels side to compare.
dm_test — the Harvey-Leybourne-Newbold (1997) small-sample-corrected statistic, computed from the documented formula and pinned as a golden.
cw_test, gw_test — documented-formula goldens from Clark-West (2007) and Giacomini-White (2006), cross-checked against an independent NumPy reference.
accuracy — the documented error-measure definitions.
var_backtest — Kupiec/Christoffersen LR statistics against first-principles NumPy/SciPy closed forms (including a hand-derived n=250/5-violation case, LR_uc = 6.0715, worked digit-by-digit in the fixture generator); the DQ statistic against a statsmodels-OLS construction on identical hit sequences; the published J.P. Morgan 1998 example (20 breaches in 252 days at 95% VaR, Jorion Value at Risk ch. 6: LR_uc = 3.91, borderline rejection vs 3.84); and seeded Monte Carlo size/power suites (all three tests ≈ nominal size on iid Bernoulli(α) hits; LR_ind/DQ reject > 90% on Markov-clustered hits with the correct unconditional rate, where Kupiec stays far lower — the separation that is the battery's reason to exist).

Fixtures: fixtures/forecast.json, fixtures/forecast_eval2.json, and fixtures/var_backtest.json.

References¶

Assimakopoulos, V. & Nikolopoulos, K. (2000). "The theta model." Int. J. Forecasting 16.
Diebold, F. & Mariano, R. (1995). "Comparing Predictive Accuracy." JBES 13.
Harvey, D., Leybourne, S. & Newbold, P. (1997). "Testing the equality of prediction mean squared errors." Int. J. Forecasting 13.
Clark, T. & West, K. (2007). "Approximately normal tests for equal predictive accuracy in nested models." J. Econometrics 138.
Giacomini, R. & White, H. (2006). "Tests of Conditional Predictive Ability." Econometrica 74.
Hyndman, R. & Koehler, A. (2006). "Another look at measures of forecast accuracy." Int. J. Forecasting 22.
Kupiec, P. (1995). "Techniques for Verifying the Accuracy of Risk Measurement Models." J. Derivatives 3(2).
Christoffersen, P. (1998). "Evaluating Interval Forecasts." Int. Economic Review 39(4).
Engle, R. & Manganelli, S. (2004). "CAViaR: Conditional Autoregressive Value at Risk by Regression Quantiles." JBES 22(4).
Jorion, P. (2007). Value at Risk, 3rd ed., ch. 6 (the J.P. Morgan 1998 worked example pinned in the goldens).

See the guide: Forecasting: Practice and Evaluation.

Runnable example¶

import numpy as np
import tsecon

rng = np.random.default_rng(1)
# monthly series: random-walk trend + a 12-period seasonal
y = 100 + np.cumsum(rng.standard_normal(120)) + 5 * np.sin(np.arange(120) * 2 * np.pi / 12)

# 1. A point forecast from the Theta method (a hard-to-beat benchmark).
fc = tsecon.theta_forecast(y[:-12], steps=12, period=12)

# 2. Score it against the held-out tail. MASE/RMSSE are scaled against an
#    in-sample seasonal-naive benchmark: < 1 beats it, > 1 loses to it.
acc = tsecon.accuracy(y[-12:], fc, insample=y[:-12], period=12)
print("RMSE:", round(acc["rmse"], 3), " MASE:", round(acc["mase"], 3))

# 3. An honest pseudo-out-of-sample backtest instead of a single split.
bt = tsecon.backtest(y, window="expanding", train=60, horizon=6,
                     forecaster="theta", period=12)
print("origins:", bt["n_origins"], " h=1 RMSE:", round(bt["accuracy"][0]["rmse"], 3))

# 4. Diebold-Mariano: do two forecasters' errors differ significantly?
e_theta = y[-12:] - fc
e_naive = y[-12:] - y[-13]                       # last training value carried forward
dm = tsecon.dm_test(e_theta, e_naive, h=1, loss="squared")
print("DM (HLN) stat:", round(dm["hln_stat"], 3), " p:", round(dm["p_value"], 3))

# 5. Giacomini-White unconditional test on any two loss series.
gw = tsecon.gw_test(e_theta ** 2, e_naive ** 2)
print("GW stat:", round(gw["gw_stat"], 3), " p:", round(gw["p_value"], 3))

# 6. Clark-West for NESTED models. Build a genuine nested pair by expanding-window
#    pseudo-OOS: small = random walk (last value), large = drift (RW + mean step).
h = 1
origin0 = 60
tgt, ys, yl = [], [], []
for t in range(origin0, len(y) - h):
    train = y[: t + 1]
    f_small = train[-1]                          # random walk
    f_large = train[-1] + np.mean(np.diff(train))  # random walk + drift (nests RW)
    tgt.append(y[t + h]); ys.append(f_small); yl.append(f_large)
tgt, ys, yl = map(np.asarray, (tgt, ys, yl))
cw = tsecon.cw_test(tgt - ys, tgt - yl, ys, yl)
print("CW stat:", round(cw["cw_stat"], 3), " p:", round(cw["p_value"], 3))

Expected output:

RMSE: 2.131  MASE: 1.268
origins: 55  h=1 RMSE: 0.936
DM (HLN) stat: -2.932  p: 0.014
GW stat: 9.377  p: 0.002
CW stat: -0.726  p: 0.766

VaR backtesting example¶

The instructive comparison: a conditional (GARCH) VaR against an unconditional (flat) VaR on volatility-clustered returns. The flat VaR gets the violation count exactly right — Kupiec cannot tell them apart — but its violations arrive in bursts during the volatile spells, which Christoffersen and the DQ test are built to catch. (The VaR paths here are filtered in-sample for brevity; a production backtest uses genuinely one-step-ahead forecasts, e.g. from a rolling re-fit.)

import numpy as np, tsecon

rng = np.random.default_rng(1)
n = 1500
r = np.empty(n); s2 = 1.0
for t in range(n):
    r[t] = np.sqrt(s2) * rng.standard_normal()
    s2 = 0.02 + 0.13 * r[t]**2 + 0.85 * s2       # GARCH(1,1) returns

fit = tsecon.garch_fit(r, vol="garch", p=1, q=1)
sigma = np.asarray(fit["conditional_volatility"])
z05 = -1.6448536269514722                        # N(0,1) 5% quantile
var_garch = z05 * sigma                          # conditional 95% VaR (negative!)
var_flat = np.full(n, z05 * r.std())             # unconditional VaR

good = tsecon.var_backtest(r, var_garch, alpha=0.05)
flat = tsecon.var_backtest(r, var_flat, alpha=0.05)
print("GARCH VaR:", good["n_violations"], "violations —",
      "p_uc=%.3f  p_ind=%.3f  p_dq=%.3f" % (good["p_uc"], good["p_ind"], good["p_dq"]))
print("flat VaR: ", flat["n_violations"], "violations —",
      "p_uc=%.3f  p_ind=%.3f  p_dq=%.3f" % (flat["p_uc"], flat["p_ind"], flat["p_dq"]))
print(flat["verdict"].split(". ")[2] + ".")

Expected output:

GARCH VaR: 82 violations — p_uc=0.414  p_ind=0.805  p_dq=0.992
flat VaR:  75 violations — p_uc=1.000  p_ind=0.014  p_dq=0.000
Reject independence at 5% (Christoffersen LR_ind = 6.064, p = 0.014): violations cluster — P(violation | violation yesterday) = 0.120 vs 0.05 under independence — so the model is too slow to update after a breach even if the overall rate is right.

75 violations where 75 were expected: the flat VaR's Kupiec p-value is 1.000. Only the dependence tests expose it — which is exactly why the battery has three legs.

Conformal prediction intervals — `conformal_forecast` / `conformal_backtest`¶

What they do. Distribution-free prediction intervals wrapped around any of the library's point forecasters: split conformal (residual-quantile calibration on held-out forecast origins with the finite-sample ceil((m+1)(1−α)) correction, symmetric or per-tail), EnbPI (Xu & Xie, ICML 2021 / IEEE TPAMI 2023 — a bootstrap ensemble of AR learners with leave-one-out out-of-bag residuals and the paper's width-minimizing β search), and ACI — adaptive conformal inference (Gibbs & Candès, NeurIPS 2021 — the online update α_{t+1} = α_t + γ(α − err_t), the go-to under distribution shift). conformal_backtest replays any of the three over a rolling evaluation window and returns the realized per-origin errors, so coverage claims about your own series are measurable rather than assumed.

Assumptions, honestly. The split guarantee (coverage ≥ 1−α in finite samples) holds under exchangeable calibration scores; h-step forecast residuals from a time series are not exchangeable, so on real data all three are approximate — that is why every claim below is a measured number, not a theorem. Calibration respects the backtest engine's leakage discipline: expanding training windows, refit at every origin, and a regression test pins that perturbing the last observation can move only its own score.

Key arguments. method="split" (default) with base= any of "theta", "naive", "drift", "mean", "seasonal_naive", "ar", "arima" — or, for split/ACI, any Python callable with the backtest callable contract; alpha=0.1; calib=n//4 residuals per horizon. method="enbpi" (AR base only, n_boot, lags, seeded and bit-reproducible). method="aci" (gamma=0.005, the paper's step size; raise it for faster adaptation — Setting B below shows why that matters).

Measured coverage (nominal 90%, one-step-ahead). Exchangeable-iid anchor at calib=20, where the +1 correction targets 19/21 ≈ 0.905: measured within 2 MC standard errors of the guarantee across 2000 reps (asserted in CI). AR(1) with iid noise (250 reps, T=300): split 0.916, EnbPI 0.888, ACI 0.912. AR(1) with GARCH noise: split 0.936, EnbPI 0.924, ACI 0.936 — marginal coverage holds, mildly conservative, though none of the methods conditions on volatility, so conditional coverage in calm vs turbulent stretches is not claimed. Under a variance shift inside the evaluation window (the published ACI scenario, 100 reps): post-shift coverage split 0.705, EnbPI 0.510, ACI at the default γ 0.794, ACI at γ=0.05 0.892 — the ACI recursion recovers what the fixed-level methods lose, exactly as Gibbs-Candès claim, and EnbPI is the most exposed because its residual window adapts slowest (lab/experiments/results/exp06.md holds the full tables; a CI-sized subset is asserted in test_conformal.py).

Failure modes. An interval is only as good as its base's residuals: a badly misspecified base gives wide-but-valid marginal intervals, not useful ones. ACI with an aggressive γ oscillates (its alpha_trajectory is returned — look at it); when α_t collapses to 0 the interval is infinite by the paper's convention, and past 1 it is the degenerate point interval. EnbPI's guarantee arguments lean on its stationarity assumptions; under shifts it degrades fastest of the three (measured above).

Validated against. Property-MC (the tables above, seeded and asserted in CI) plus the exchangeable-case exactness anchor, and — for the split leg — a runnable third-party cross-check: on identical residuals our finite-sample-corrected quantile reproduces mapie's SplitConformalRegressor interval half-width to 1e-12 relative (test_mapie_split_quantile_cross_check, skipped automatically where mapie is absent). EnbPI and ACI are cross-checked against mapie 1.5.0's TimeSeriesRegressor (0.10.0; fixtures/conformal_mapie.json, generate_conformal_mapie_fixtures.py, pinned by test_conformal_mapie.py), graded per leg. ACI is an exact cross-check — with the same prefit linear point forecaster (handed to conformal_backtest as a Python callable), the same 75-residual sliding window, the same step size and mapie's AbsoluteConformityScore(sym=True), tsecon reproduces mapie's per-origin bounds at 1e-12 relative and its miss indicators and α_t trajectory exactly, on both γ = 0.05 (realized coverage 0.8933, final α_t 0.0750) and γ = 0.005 (0.9067, 0.1025). The sym=True is essential and is the finding worth recording: TimeSeriesRegressor defaults to sym=False, which builds the interval from a pair of signed-residual quantiles (β = α_t/2 below, 1 − α_t + β above) and is asymmetric about the point forecast — not the absolute-score construction the ACI paper specifies and this library implements. Three differences remain and none fires on these runs: mapie clips α_t to [0, 1] (the recursion here is unclipped), mapie counts a target exactly on a bound as a miss (here it is covered), and when ceil((m+1)(1−α_t)) runs past the window mapie returns an infinite bound while conformal_backtest refuses the call naming the level and the residuals it would need (on both stored runs the worst order index is 74 of 75). EnbPI is a statistical cross-check only: same algorithm, different bootstrap generators (mapie's BlockBootstrap(length=1) on a NumPy RandomState vs Philox), the +1 finite-sample correction in mapie versus the paper's empirical quantile here, and different β grids — so on the same AR(1) series and design the two online runs are compared in distribution: realized coverage within 0.06, mean width within 10%, mean absolute centre gap below 0.05. Measured (printed by the test): with the β line search off, coverage 0.9067 both sides, mean width 3.4024 here vs 3.4830 in mapie (ratio 0.977), mean absolute centre gap 0.0106; with it on, coverage 0.9067 vs 0.9333, width 3.3117 vs 3.6106 (ratio 0.917), same centre gap. No exact EnbPI pin is possible without one side adopting the other's random stream and quantile convention, and the test says so.

References. Vovk, Gammerman & Shafer (2005); Xu & Xie (2021, 2023); Gibbs & Candès (2021).

Many models at once — `spa_test` / `stepm_test` / `model_confidence_set`¶

What they do. Every test above this section compares two forecasters. These three compare many, and correct for the search that a column of pairwise p-values silently invites.

Function	Question it answers
`spa_test`	Does the best of \(m\) models beat the benchmark, once the search over all \(m\) is accounted for? (White's 2000 Reality Check; Hansen's 2005 SPA)
`stepm_test`	Which models beat the benchmark, at a controlled family-wise error rate? (Romano-Wolf 2005)
`model_confidence_set`	Which models are indistinguishable from the best? — no benchmark needed (Hansen-Lunde-Nason 2011)

All three consume the same object: a \(T \times m\) table of per-origin losses, which is what a backtest loop already produces (see the runnable example below, and guide 5 §"Many models at once").

spa_test. With \(d_{t,k} = L_t(\text{benchmark}) - L_t(\text{model } k)\) (positive favours the model), the null is \(H_0: \max_k \mathbb{E}[d_k] \le 0\) and the statistic is \(\sqrt{n}\max_k \bar d_k/\omega_k\) (studentize=True, Hansen's SPA) or \(\sqrt{n}\max_k \bar d_k\) (studentize=False, White's Reality Check). The null distribution is a block bootstrap of the whole loss-differential panel with rows resampled together, so the cross-model dependence that shapes a maximum survives. Hansen's contribution is the re-centring: p_value_upper re-centres every model (White's original — conservative, because junk models inflate the simulated maximum), p_value_consistent leaves models worse than the benchmark by more than \(\sqrt{2\omega_k^2\log\log n / n}\) un-centred (the recommended p-value, also returned as p_value), and p_value_lower re-centres none with a negative sample mean (the liberal bound). Always \(p_\text{lower} \le p_\text{consistent} \le p_\text{upper}\), and the crate asserts it.

model_confidence_set. Sequential elimination: test equal predictive ability across the models still in the set with \(T_R\) (method="R", the maximum standardized pairwise mean difference — HLN's recommended default) or \(T_{\max}\) (method="max", each model against the cross-sectional mean); if the test rejects at size, eliminate the worst and repeat. Each model's MCS p-value is the running maximum of the step p-values along the elimination path, so one call describes every size at once — the set at any \(\alpha\) is \(\{k : p_\text{MCS}(k) > \alpha\}\) and the sets are nested in \(\alpha\) (asserted).

Assumptions. Stationary, index-aligned loss series over a common evaluation sample (the same origins, the same horizon, the same scheme — backtest gives you this for free, and the example below asserts it). The block bootstrap needs the dependence to be short-memory relative to the block length. Nothing here knows where the losses came from: nested models are as invisible to these tests as they are to dm_test, so the cw_test warning still applies column by column. Non-finite losses are refused, never skipped.

When to use, and when not. Use spa_test the moment you have compared more than two or three specifications against a benchmark; use model_confidence_set when there is no natural benchmark and the honest answer is a set rather than a winner; use stepm_test when you must name the winners and defend each claim. Do not use them to rescue a single DM comparison (that is dm_test's job and its p-value is the right one), and do not use them on \(m\) models that are \(m\) tunings of one model — every column should be a forecast you would have been willing to ship.

Key arguments and defaults, and why. block_size=None takes the Politis-White (2004) / Patton-Politis-White (2009) optimal length of each column and averages, reporting it in block_size with block_size_auto True; it is a reasonable default but the block length is the single setting that moves these p-values, so report it and check a hand-set value near \(n^{1/3}\). reps=1000 is the conventional count — raise it when a p-value sits near your decision boundary, because its own Monte Carlo standard error is \(\sqrt{p(1-p)/\text{reps}}\). bootstrap="stationary" (geometric block lengths) with "circular" and "moving_block" available; studentize=True is Hansen's statistic (see the size numbers below before relying on it in a short sample); nested=False uses Hansen's eq. 9 kernel for \(\omega_k^2\) rather than a bootstrap-of-the-bootstrap; size=0.10 for the MCS is HLN's own illustration level; seed=0 — every path is seeded, and bit-identical at any thread count (asserted at 1–4 threads).

How to read the output. spa_test returns statistic (on the \(\sqrt{n}\) scale), best_model, the three p-values, crit_lower / crit_consistent / crit_upper at crit_levels \(=[0.90, 0.95, 0.99]\), mean_loss_diff, loss_diff_var, recentered (which models the consistent p-value left un-centred) and the full replicate vectors. A large gap between p_value_consistent and p_value_upper is informative, not noise: it says your search contained models far worse than the benchmark, which is exactly what the Reality Check over-penalizes. model_confidence_set returns included / excluded, mcs_p_values (per model), elimination_order with the step_p_values that removed each model, and mean_losses. Read elimination_order as a diagnostic, not a ranking — \(T_R\) eliminates the worse member of the most standardized pair, so a model with an enormous loss variance can survive longer than one with a larger mean loss. The ranking is mean_losses; the evidence is mcs_p_values.

Failure modes. A model whose losses equal the benchmark's in every period has a loss differential of exactly zero — it is the benchmark — and spa_test refuses it by name. For the MCS, a duplicated loss column is degenerate under method="R" only: the pairwise bootstrap variance of two identical columns is exactly zero, so \(T_R\) between them is 0/0 and the panel is refused. Under method="max" nothing is 0/0 — the duplicates share a statistic and leave the set together in one step — so the panel is accepted, and the refusal fires only when a remaining model's bootstrap standard deviation really does collapse. Both behaviours beat the reference's, which the fixture generator measures rather than assumes: on a panel with two identical loss columns arch 8.0.0's MCS warns about the 0/0 division and then raises IndexError under method="R", and under method="max" returns nothing within a 45-second, 2 GiB budget. The MCS set is never empty, but a set containing everything means the evaluation sample is too short to separate anything — look at n before writing it up. A p-value of exactly 0.000 means no replicate exceeded the observed statistic, i.e. \(p < 1/\text{reps}\); it is not evidence stronger than your replication count.

Validated against, and the honest grade. arch 8.0.0 (arch.bootstrap.{SPA, RealityCheck, StepM, MCS}) — an independent package — at four graded legs, plus seeded Monte Carlo for what no package can pin.

Exact, on the reference's own draws. arch draws from NumPy's generator, which the Rust core cannot reproduce, so fixtures/spa.json and fixtures/mcs.json store arch's resample index arrays and the crate goldens replay them through internal *_with_indices entry points. On those resamples tsecon reproduces arch bit for bit (every case with \(m \ge 2\)): the mean loss differentials, all 200 replicate statistics under each of the three re-centrings, the observed statistic, the critical values and the three p-values for SPA/RC; the mean losses, the pairwise variance matrix, the elimination order, the included/excluded sets and the MCS p-values for the MCS; the superior sets for StepM at sizes 0.05 and 0.10. Hansen's kernel variances are pinned at 1e-13 (one pow in the kernel weights may differ by an ulp across platforms) and the single-model case at 1e-13 (NumPy coalesces a \(T \times 1\) panel into a pairwise 1-D reduction; measured and asserted in the generator). The generator certifies a min_gap \(> 10^{-9}\) between every replicate and the observed statistic, so the p-values are exact by construction rather than by luck.
Documented-formula, for the studentized path. No package studentizes — arch 8.0's studentize flag is inert, which the generator asserts by comparing its output with the flag on and off (RealityCheck is literally class RealityCheck(SPA): pass) — so studentize=True is pinned at 1e-12 relative against a NumPy transcription of Hansen's (2005) formulas on the same resamples, with p-values exact. A second recorded finding: arch's StepM raises ValueError whenever every model is declared superior over two or more steps; tsecon stops the loop with all \(m\) models in the set, and the fixture stores which cases hit it.
Resampling conventions, measured not assumed. The generator replays arch's own raw draws through the rule tsecon-bootstrap documents and gets arch's index arrays back element for element, for all three schemes at three \((n, \text{block\_size})\) settings. One difference, documented: arch restarts a stationary block on \(u \le p\), tsecon on \(u < p\) — they disagree only on the null event \(u = p\) (probability \(2^{-53}\) per step; the exact ties are counted and were zero). Block length, consecutive layout, the modulo-\(n\) wrap for circular, its absence for moving-block, the start ranges and the truncation to \(n\) are identical.
Monte Carlo, for the seeded public path. spa_test and model_confidence_set at 4000 replications land within 0.05 of arch's 4000-replication p-values (two independent bootstraps differ by ~0.011 standard deviations at \(p = 0.5\), so 0.05 is ~4.5 sd), and reproduce arch's MCS set exactly on designs the generator certified are separated from size.

Measured size, power and coverage (seeded, asserted in CI). Design: six exchangeable squared-error loss columns whose forecast errors share an AR(1) common component — the least favourable null, every model exactly as good as the benchmark. \(n = 200\), \(m = 5\), \(B = 300\), 400 Monte Carlo replications.

Rejection of `p_value_consistent` at 5%	tsecon	arch, same design
iid losses, `block_size=2`	0.038	0.038
AR(0.5) losses, `block_size=3`	0.080	0.058
AR(0.5) losses, `block_size=8`	0.068	0.050
AR(0.5) losses, `block_size=14`	0.080	0.070

These are the un-studentized rates, because that is all arch computes; they agree with the reference at every level and the test asserts it within 4 Monte Carlo standard errors. The method is mildly over-sized under dependence at \(n = 200\) — in both libraries.

The studentized default costs more size, and we say so. With studentize=True the same design rejects at 0.123 (AR(0.5), Politis-White block) and 0.135 (block 8) against a nominal 0.05. The mechanism is Hansen's own construction, not a defect here: eqs. 5–8 divide the observed statistic and every bootstrap replicate by the same \(\omega_k\), so the bootstrap maximum carries none of \(\omega_k\)'s sampling error while the data maximum does. It is a finite-sample effect and shrinks with the sample — 0.123 → 0.093 and 0.135 → 0.088 at \(n = 800\) (asserted) — but in a short evaluation sample prefer studentize=False, or read p_value_upper alongside. No package computes the studentized statistic, so this leg is measured, not validated: the numbers above are the claim.

Power. Against a benchmark one of whose competitors has a 0.6× error scale (\(n = 200\), 200 replications): the consistent p-value rejects at 5% in 0.985 of samples, White's upper p-value in 0.985, and best_model identifies the right column in 0.995.

MCS coverage. Design: four models, two with exactly equal expected loss (so the best set \(M^*\) has two elements), one slightly worse, one clearly worst; size=0.10, \(B = 300\), 1000 replications.

	\(P(M^* \subseteq\) set\()\)	\(P\)(a best model in set)	\(P\)(worst excluded)	mean set size
\(T_R\), \(n=150\)	0.870 (arch 0.881)	0.930 (0.947)	1.000 (1.000)	2.18 (2.23)
\(T_{\max}\), \(n=150\)	0.895 (0.909)	0.944 (0.960)	1.000 (1.000)	2.21 (2.27)
\(T_R\), \(n=600\)	0.872 (0.879)	0.944 (0.944)	1.000 (1.000)	1.88 (1.89)
\(T_{\max}\), \(n=600\)	0.873 (0.879)	0.945 (0.944)	1.000 (1.000)	1.88 (1.89)

Read that table honestly. HLN's Theorem 1 is asymptotic and about the whole set \(M^*\); containing a two-element \(M^*\) runs at ~0.87 against a nominal 0.90 and does not improve from \(n = 150\) to \(n = 600\). Two finite-sample effects cause it — an early step can drop one of the two best models while the inferior ones are still in the set, and the final test between two identical models rejects at its own size — and arch shows the same frequencies to within Monte Carlo error, which is what the crate test asserts. The easier event, "a best model is in the set", does hold at the nominal level (0.93–0.945). In practice: an MCS containing two models does not promise 90% that both belong there.

References. White, H. (2000), "A Reality Check for Data Snooping," Econometrica 68(5), 1097–1126. Hansen, P. R. (2005), "A Test for Superior Predictive Ability," JBES 23(4), 365–380. Romano, J. P. and M. Wolf (2005), "Stepwise Multiple Testing as Formalized Data Snooping," Econometrica 73(4), 1237–1282. Hansen, P. R., A. Lunde and J. M. Nason (2011), "The Model Confidence Set," Econometrica 79(2), 453–497. Politis, D. N. and J. P. Romano (1994), "The Stationary Bootstrap," JASA 89, 1303–1313. Politis and White (2004) and Patton, Politis and White (2009) for the block length.

Runnable example — a `backtest` loss table into the model confidence set¶

import numpy as np
import tsecon

rng = np.random.default_rng(7)                    # quarterly: trend + season + AR noise
n = 160
t_idx = np.arange(n)
season = 4.0 * np.array([1.0, -0.4, 0.6, -1.2])[t_idx % 4]
noise = np.zeros(n)
e = rng.standard_normal(n)
for i in range(1, n):
    noise[i] = 0.6 * noise[i - 1] + 1.5 * e[i]
y = 50 + 0.3 * t_idx + season + noise

# 1. One backtest per forecaster, under ONE scheme, so the origins align.
names = ["naive", "drift", "mean", "seasonal_naive", "theta"]
losses, origins = [], None
for fc in names:
    seasonal = {"period": 4} if fc in ("seasonal_naive", "theta") else {}
    bt = tsecon.backtest(y, window="expanding", train=80, horizon=1,
                         forecaster=fc, insample_period=4, **seasonal)
    assert origins is None or bt["origins"] == origins
    origins = bt["origins"]
    err = np.array(bt["targets"][0]) - np.array(bt["forecasts"][0])
    losses.append(err ** 2)
L = np.column_stack(losses)                       # 80 origins x 5 models

# 2. Which models are indistinguishable from the best?
mcs = tsecon.model_confidence_set(L, size=0.10, reps=2000, seed=0)
print(f"{'forecaster':16s} mean loss   MCS p   in the 90% set?")
for k, nm in enumerate(names):
    print(f"{nm:16s} {mcs['mean_losses'][k]:8.2f}   {mcs['mcs_p_values'][k]:5.3f}   "
          f"{'yes' if k in mcs['included'] else 'no'}")
print("eliminated, worst first:", [names[k] for k in mcs["elimination_order"]])
print("block length used:", mcs["block_size"], "(automatic:", mcs["block_size_auto"], ")")
# forecaster       mean loss   MCS p   in the 90% set?
# naive               45.73   0.000   no
# drift               46.06   0.000   no
# mean               351.11   0.000   no
# seasonal_naive       5.93   0.024   no
# theta                4.03   1.000   yes
# eliminated, worst first: ['drift', 'naive', 'mean', 'seasonal_naive', 'theta']
# block length used: 8 (automatic: True )

# 3. Against a named benchmark: does anything beat the seasonal naive?
spa = tsecon.spa_test(L[:, 3], L[:, [0, 1, 2, 4]], reps=2000, seed=0)
print(f"SPA statistic {spa['statistic']:.3f}; p = {spa['p_value']:.3f} consistent, "
      f"{spa['p_value_upper']:.3f} upper (White); best column {spa['best_model']}")
print("StepM superior at FWER 5%:",
      tsecon.stepm_test(L[:, 3], L[:, [0, 1, 2, 4]], size=0.05, reps=2000,
                        seed=0)["superior_models"])
# SPA statistic 2.151; p = 0.019 consistent, 0.046 upper (White); best column 3
# StepM superior at FWER 5%: [3]

The 90% model confidence set is a single model — the strongest verdict an MCS can give — and seasonal_naive leaves on a p-value rather than on eyeballing a 5.93-against-4.03 gap. The gap between the consistent p-value (0.019) and White's (0.046) is Hansen's re-centring earning its keep: three of the four candidates are far worse than the benchmark, and the Reality Check pays for including them.

Model card — Forecasting and forecast evaluation¶

What it estimates¶

Assumptions¶

When to use¶

Key arguments and defaults¶

How to read the output¶

Failure modes¶

Bring your own model — Python-callable forecasters¶

Validated against¶

References¶

Runnable example¶

VaR backtesting example¶

Conformal prediction intervals — conformal_forecast / conformal_backtest¶

Many models at once — spa_test / stepm_test / model_confidence_set¶

Runnable example — a backtest loss table into the model confidence set¶

Conformal prediction intervals — `conformal_forecast` / `conformal_backtest`¶

Many models at once — `spa_test` / `stepm_test` / `model_confidence_set`¶

Runnable example — a `backtest` loss table into the model confidence set¶