TeX 62.3%
Python 37.6%
1# Author: Simon-Pierre Boucher — simon-pierre.boucher@uqo.ca2"""Forecast evaluation: point and probabilistic scores, Diebold–Mariano (HAC), Giacomini–White, Model Confidence Set, block bootstrap."""3from __future__ import annotations45import numpy as np6from scipy import stats789# ---- point scores -----------------------------------------------------------------------------10def point_scores(y: np.ndarray, f: np.ndarray, y_naive: np.ndarray | None = None) -> dict:11 e = f - y; ok = np.isfinite(e); e = e[ok]12 out = {"n": int(ok.sum()), "mae": float(np.mean(np.abs(e))), "rmse": float(np.sqrt(np.mean(e ** 2))), "medae": float(np.median(np.abs(e))), "bias": float(np.mean(e))}13 if y_naive is not None:14 en = (y_naive - y)[ok]; out["mase"] = float(np.mean(np.abs(e)) / max(np.mean(np.abs(en)), 1e-12))15 return out161718# ---- probabilistic scores ---------------------------------------------------------------------19def pinball(y: np.ndarray, q: np.ndarray, tau: float) -> np.ndarray:20 d = y - q21 return np.maximum(tau * d, (tau - 1) * d)222324def crps_from_quantiles(y: np.ndarray, Q: np.ndarray, taus: np.ndarray) -> np.ndarray:25 """CRPS approximated by the average pinball loss over a quantile grid (×2), Laio & Tamea (2007) relation."""26 return 2 * np.mean(np.stack([pinball(y, Q[:, j], t) for j, t in enumerate(taus)], axis=1), axis=1)272829def interval_scores(y: np.ndarray, lo: np.ndarray, hi: np.ndarray, alpha: float = 0.1) -> dict:30 cov = np.mean((y >= lo) & (y <= hi)); width = np.mean(hi - lo)31 wis = np.mean((hi - lo) + 2 / alpha * (lo - y) * (y < lo) + 2 / alpha * (y - hi) * (y > hi))32 return {"coverage": float(cov), "width": float(width), "interval_score": float(wis)}333435def crps_gaussian(y: np.ndarray, mu: np.ndarray, sigma: np.ndarray) -> np.ndarray:36 z = (y - mu) / sigma37 return sigma * (z * (2 * stats.norm.cdf(z) - 1) + 2 * stats.norm.pdf(z) - 1 / np.sqrt(np.pi))383940# ---- tests ------------------------------------------------------------------------------------41def hac_var(d: np.ndarray, lags: int) -> float:42 d = d - d.mean(); n = len(d); v = np.sum(d * d) / n43 for L in range(1, lags + 1):44 w = 1 - L / (lags + 1); v += 2 * w * np.sum(d[L:] * d[:-L]) / n45 return max(v, 1e-18)464748def diebold_mariano(loss_a: np.ndarray, loss_b: np.ndarray, h: int = 1, small_sample: bool = True) -> dict:49 """DM test of H0: E[loss_a − loss_b] = 0 with Newey–West variance (h−1 lags) and the Harvey–Leybourne–Newbold correction."""50 d = loss_a - loss_b; d = d[np.isfinite(d)]; n = len(d)51 if n < 10:52 return {"dm": np.nan, "p": np.nan, "n": n, "mean_diff": float(np.mean(d)) if n else np.nan}53 v = hac_var(d, max(h - 1, 0)); dm = d.mean() / np.sqrt(v / n)54 if small_sample:55 dm *= np.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n)56 p = 2 * (1 - stats.t.cdf(abs(dm), df=n - 1))57 else:58 p = 2 * (1 - stats.norm.cdf(abs(dm)))59 return {"dm": float(dm), "p": float(p), "n": n, "mean_diff": float(d.mean())}606162def giacomini_white(loss_a: np.ndarray, loss_b: np.ndarray, Z: np.ndarray | None = None, h: int = 1) -> dict:63 """Conditional predictive ability test: regress d_t on instruments Z_{t-1} (constant + lagged d by default); Wald test with HAC."""64 d = loss_a - loss_b; ok = np.isfinite(d); d = d[ok]65 if Z is None:66 Z = np.column_stack([np.ones(len(d) - 1), d[:-1]]); dd = d[1:]67 else:68 Z = Z[ok][1:]; dd = d[1:]69 n, q = Z.shape70 if n < 10 + q:71 return {"stat": np.nan, "p": np.nan, "n": n}72 Zd = Z * dd[:, None]; m = Zd.mean(axis=0)73 S = np.zeros((q, q)); lags = max(h - 1, 0)74 Zc = Zd - m75 S = Zc.T @ Zc / n76 for L in range(1, lags + 1):77 w = 1 - L / (lags + 1); G = Zc[L:].T @ Zc[:-L] / n; S += w * (G + G.T)78 try:79 stat = float(n * m @ np.linalg.solve(S, m))80 except np.linalg.LinAlgError:81 return {"stat": np.nan, "p": np.nan, "n": n}82 return {"stat": stat, "p": float(1 - stats.chi2.cdf(stat, q)), "n": n}838485def stationary_bootstrap_indices(n: int, mean_block: float, rng: np.random.Generator) -> np.ndarray:86 p = 1 / mean_block; idx = np.empty(n, int); i = rng.integers(n)87 for t in range(n):88 if t > 0 and rng.random() > p:89 i = (i + 1) % n90 else:91 i = rng.integers(n)92 idx[t] = i93 return idx949596def model_confidence_set(losses: np.ndarray, names: list[str], alpha: float = 0.10, B: int = 1000, mean_block: float = 3.0, seed: int = 0) -> dict:97 """Hansen–Lunde–Nason MCS with the T_max statistic and the stationary bootstrap. losses: (T, M)."""98 rng = np.random.default_rng(seed)99 L = losses[np.all(np.isfinite(losses), axis=1)]; T, M = L.shape100 alive = list(range(M)); pvals = {}; order = []101 boots = [stationary_bootstrap_indices(T, mean_block, rng) for _ in range(B)]102 while len(alive) > 1:103 Lk = L[:, alive]; dbar = Lk.mean(axis=0) - Lk.mean()104 # bootstrap distribution of the deviations105 bd = np.array([(Lk[b].mean(axis=0) - Lk[b].mean()) - dbar for b in boots])106 var = np.maximum(bd.var(axis=0), 1e-18)107 t_stat = dbar / np.sqrt(var); tmax = t_stat.max()108 tmax_b = (bd / np.sqrt(var)).max(axis=1)109 p = float(np.mean(tmax_b >= tmax))110 worst = int(np.argmax(t_stat)); name = names[alive[worst]]111 pvals[name] = max(p, max(pvals.values()) if pvals else 0.0); order.append(name)112 if p >= alpha:113 break114 alive.pop(worst)115 for i in alive:116 pvals.setdefault(names[i], 1.0)117 return {"mcs_pvalues": pvals, "in_mcs": [n for n, p in pvals.items() if p >= alpha], "eliminated_in_order": order}118119120def block_bootstrap_ci(values: np.ndarray, stat=np.mean, B: int = 1000, mean_block: float = 3.0, seed: int = 0, level: float = 0.95) -> tuple[float, float]:121 rng = np.random.default_rng(seed); v = values[np.isfinite(values)]; n = len(v)122 s = np.array([stat(v[stationary_bootstrap_indices(n, mean_block, rng)]) for _ in range(B)])123 return float(np.quantile(s, (1 - level) / 2)), float(np.quantile(s, 1 - (1 - level) / 2))124