Repository Public Research preview
1"""Forecast-comparison inference: Diebold–Mariano with HAC variance (Newey–West Bartlett kernel, bandwidth2h-1+floor(4(n/100)^(2/9)) and the Harvey–Leybourne–Newbold small-sample correction), the Model Confidence Set3(Hansen, Lunde & Nason 2011) and the Superior Predictive Ability test (Hansen 2005) through arch.bootstrap (stationary4bootstrap with expected block length 22 days), and a stationary block bootstrap for cross-asset pooled loss5differentials that resamples DATE blocks jointly across assets (preserves cross-sectional dependence)."""6from __future__ import annotations78import math910import numpy as np11import pandas as pd12from scipy import stats131415def hac_variance(d: np.ndarray, bandwidth: int) -> float:16 d = np.asarray(d, float)17 d = d - d.mean()18 n = len(d)19 v = float(np.dot(d, d) / n)20 for j in range(1, bandwidth + 1):21 w = 1.0 - j / (bandwidth + 1.0)22 v += 2.0 * w * float(np.dot(d[j:], d[:-j]) / n)23 return max(v, 1e-300)242526def dm_test(loss_a: np.ndarray, loss_b: np.ndarray, h: int = 1, bandwidth: int | None = None) -> dict:27 """H0: E[loss_a - loss_b] = 0. Negative statistic favours a. Returns dict(stat, p, mean_diff, n, bandwidth)."""28 d = np.asarray(loss_a, float) - np.asarray(loss_b, float)29 d = d[np.isfinite(d)]30 n = len(d)31 if n < 30:32 return dict(stat=np.nan, p=np.nan, mean_diff=np.nan, n=n, bandwidth=0)33 bw = bandwidth if bandwidth is not None else (h - 1 + int(math.floor(4 * (n / 100.0) ** (2.0 / 9.0))))34 v = hac_variance(d, bw)35 stat = d.mean() / math.sqrt(v / n)36 corr = math.sqrt((n + 1 - 2 * h + h * (h - 1) / n) / n) # HLN (1997)37 stat_c = stat * corr38 p = 2 * stats.t.sf(abs(stat_c), df=n - 1)39 return dict(stat=float(stat_c), p=float(p), mean_diff=float(d.mean()), n=int(n), bandwidth=int(bw))404142def mcs(losses: pd.DataFrame, size: float = 0.10, reps: int = 2000, block: int = 22, seed: int = 46, method: str = "max") -> pd.DataFrame:43 """losses: DataFrame (origins x models), complete cases only. Returns DataFrame(model, pvalue, included)."""44 from arch.bootstrap import MCS45 L = losses.dropna()46 if L.shape[0] < 50 or L.shape[1] < 2:47 return pd.DataFrame({"model": L.columns, "pvalue": np.nan, "included": False})48 m = MCS(L, size=size, reps=reps, block_size=block, method=method, bootstrap="stationary", seed=seed)49 m.compute()50 pv = m.pvalues["Pvalue"]51 return pd.DataFrame({"model": pv.index, "pvalue": pv.values, "included": [c in set(m.included) for c in pv.index]})525354def spa(bench: np.ndarray, models: pd.DataFrame, reps: int = 2000, block: int = 22, seed: int = 46) -> dict:55 """Hansen (2005) SPA: H0 = no model in `models` beats the benchmark. Returns consistent/upper/lower p-values."""56 from arch.bootstrap import SPA57 M = models.copy()58 M["_b"] = bench59 M = M.dropna()60 s = SPA(M["_b"].values, M.drop(columns="_b"), reps=reps, block_size=block, bootstrap="stationary", seed=seed)61 s.compute()62 pv = s.pvalues63 return {"lower": float(pv["lower"]), "consistent": float(pv["consistent"]), "upper": float(pv["upper"]), "n": int(len(M))}646566def block_bootstrap_mean_ci(D: pd.DataFrame, reps: int = 2000, block: int = 22, seed: int = 46, alpha: float = 0.10) -> dict:67 """D: DataFrame (dates x assets) of loss differentials (model - benchmark). Stationary bootstrap over DATE blocks,68 assets resampled jointly. Returns the pooled mean, the (alpha/2, 1-alpha/2) percentile CI and a two-sided p-value for H0: mean = 0."""69 X = D.values70 n = X.shape[0]71 rng = np.random.default_rng(seed)72 stat = float(np.nanmean(X))73 boots = np.empty(reps)74 p = 1.0 / block75 for b in range(reps):76 idx = np.empty(n, dtype=int)77 i = 078 while i < n:79 start = rng.integers(0, n)80 length = rng.geometric(p)81 for k in range(length):82 if i >= n:83 break84 idx[i] = (start + k) % n85 i += 186 boots[b] = np.nanmean(X[idx])87 lo, hi = np.percentile(boots, [100 * alpha / 2, 100 * (1 - alpha / 2)])88 centered = boots - boots.mean()89 pval = float(np.mean(np.abs(centered) >= abs(stat)))90 return {"mean": stat, "ci_low": float(lo), "ci_high": float(hi), "p": pval, "n_dates": int(n), "n_assets": int(X.shape[1])}91