1#!/usr/bin/env python32"""Stage 2 of the evaluation: per-instrument and pooled statistics from the wide loss tables.3usage: analyze.py --run main [--period test|validation] [--symbols ...] [--reps 2000]4Outputs (results/metrics/analysis_<run>_<period>/): per_symbol/<sym>.json (pairwise ratios, DM, coverage, MCS, SPA, regimes,5tail risk, failure accounting), pooled.csv/json, mcs_counts.csv, tailrisk.csv, compute.csv, and results/tables/*.csv copies.6Everything is computed on the requested period only (the test period is used only after the protocol freeze)."""7from __future__ import annotations89import argparse10import glob11import json12import sys13from pathlib import Path1415import numpy as np16import pandas as pd17import yaml1819ROOT = Path(__file__).resolve().parents[1]20sys.path.insert(0, str(ROOT / "src"))21from wp46.diagnostics import backtest_table, var_es_forecasts # noqa: E40222from wp46.inference import block_bootstrap_mean_ci, dm_test, mcs, spa # noqa: E40223from wp46.preprocessing import load_panel # noqa: E4022425PROT = yaml.safe_load(open(ROOT / "configs" / "forecast_protocol.yaml"))26REG = yaml.safe_load(open(ROOT / "configs" / "model_registry.yaml"))["models"]27INST = yaml.safe_load(open(ROOT / "configs" / "instruments.yaml"))28CLASS_OF = {s: c for c, syms in (INST.get("universe") or INST["candidates"]).items() for s in syms}29PRIMARY_PROXY, PRIMARY_LOSS = PROT["primary_proxy"], PROT["primary_loss"]30HORIZONS = tuple(PROT["horizons"])31MIN_COVERAGE_MCS = 0.90323334def spec_meta(spec_id: str) -> dict:35 mid, dist, rest = spec_id.split("__")36 m = REG.get(mid, {})37 return {"model_id": mid, "dist": dist, "family": m.get("family", "?"), "tier": m.get("tier", "?"), "variant": rest}383940def bench_id(specs) -> str:41 b = [s for s in specs if s.startswith("garch11__normal__")]42 if not b:43 raise RuntimeError("benchmark GARCH(1,1)-normal missing")44 return b[0]454647def load_symbol(run: str, sym: str):48 d = ROOT / "results" / "metrics" / f"losses_{run}"49 W = pd.read_parquet(d / f"{sym}.parquet"); W["origin"] = pd.to_datetime(W["origin"])50 T = pd.read_parquet(d / f"targets_{sym}.parquet"); T["origin"] = pd.to_datetime(T["origin"])51 S = pd.read_parquet(d / f"status_{sym}.parquet"); S["origin"] = pd.to_datetime(S["origin"])52 G = pd.read_parquet(d / f"fits_{sym}.parquet") if (d / f"fits_{sym}.parquet").exists() else pd.DataFrame()53 return W, T.set_index("origin"), S.set_index("origin"), G545556def block(W: pd.DataFrame, kind: str, proxy: str = "", loss: str = "", h: int = 1) -> pd.DataFrame:57 m = (W["kind"] == kind) & (W["proxy"] == proxy) & (W["loss"] == loss) & (W["h"] == h)58 return W[m].drop(columns=["kind", "proxy", "loss", "h", "symbol"]).set_index("origin").sort_index()596061def regimes_for(panel: pd.DataFrame, origins: pd.DatetimeIndex, val_start: str, val_end: str) -> pd.DataFrame:62 p = panel.set_index(pd.to_datetime(panel["session"]))63 rv = p["rv_5m_ss"].where(p["valid_rm"]).fillna(p["r"] ** 2)64 trail22 = rv.rolling(22, min_periods=15).mean()65 trail66 = rv.rolling(66, min_periods=40).mean()66 v = trail22[(trail22.index >= val_start) & (trail22.index <= val_end)]67 q = v.quantile([1 / 3, 2 / 3]).values if len(v) > 50 else trail22.quantile([1 / 3, 2 / 3]).values68 v66 = trail66[(trail66.index >= val_start) & (trail66.index <= val_end)]69 q66 = v66.quantile([1 / 3, 2 / 3]).values if len(v66) > 50 else trail66.quantile([1 / 3, 2 / 3]).values70 R = pd.DataFrame(index=origins)71 t22 = trail22.reindex(origins); t66 = trail66.reindex(origins)72 R["vol_tercile"] = np.where(t22 <= q[0], "low", np.where(t22 <= q[1], "mid", "high"))73 R["vol_tercile66"] = np.where(t66 <= q66[0], "low", np.where(t66 <= q66[1], "mid", "high"))74 R["prev_sign"] = np.where(p["r"].reindex(origins) < 0, "neg", "pos")75 ep = PROT["regimes"]["calendar_episodes"]76 R["episode"] = "none"77 for name, w in ep.items():78 R.loc[(origins >= w["start"]) & (origins <= w["end"]), "episode"] = name79 return R808182def analyze_symbol(run: str, sym: str, period: str, reps: int) -> dict:83 W, T, S, G = load_symbol(run, sym)84 panel = load_panel(sym)85 per = PROT["periods"][period]86 lo, hi = pd.Timestamp(per["start"]), pd.Timestamp(per["end"])87 out = {"symbol": sym, "asset_class": CLASS_OF.get(sym), "period": period, "horizons": {}, "regimes": {}, "secondary": {}}88 specs = None89 for h in HORIZONS:90 L = block(W, "loss", PRIMARY_PROXY, PRIMARY_LOSS, h)91 L = L[(L.index >= lo) & (L.index <= hi)]92 specs = list(L.columns)93 b = bench_id(specs)94 y = T[f"__target__{PRIMARY_PROXY}_{h}"].reindex(L.index)95 tv = np.isfinite(y.values)96 n_target = int(tv.sum())97 res = {"n_origins": int(len(L)), "n_valid_targets": n_target, "benchmark": b, "specs": {}}98 for s in specs:99 pair = np.isfinite(L[s].values) & np.isfinite(L[b].values)100 cov = float(np.isfinite(L[s].values[tv]).mean()) if n_target else np.nan101 if pair.sum() >= 30:102 ls, lb = L[s].values[pair], L[b].values[pair]103 dm = dm_test(ls, lb, h)104 res["specs"][s] = {"n_pair": int(pair.sum()), "coverage": cov, "mean": float(ls.mean()), "mean_bench": float(lb.mean()),105 "ratio": float(ls.mean() / lb.mean()), "dm_stat": dm["stat"], "dm_p": dm["p"],106 "own_mean": float(np.nanmean(L[s].values)), "own_n": int(np.isfinite(L[s].values).sum())}107 else:108 res["specs"][s] = {"n_pair": int(pair.sum()), "coverage": cov, "ratio": np.nan, "dm_stat": np.nan, "dm_p": np.nan}109 # MCS / SPA per tier on the common-success sample of well-covered specs110 for tier in ("A", "B"):111 cand = [s for s in specs if spec_meta(s)["tier"] == tier and res["specs"][s]["coverage"] >= MIN_COVERAGE_MCS]112 if tier == "B":113 cand = cand + ([b] if b not in cand else []) # GARCH(1,1) reference inside the realized tier (labelled)114 if len(cand) < 2:115 continue116 M = L[cand].dropna()117 res[f"common_sample_{tier}"] = {"n": int(len(M)), "specs": cand, "mean": {s: float(M[s].mean()) for s in cand},118 "ratio": {s: float(M[s].mean() / M[b].mean()) for s in cand} if b in M else {}}119 if len(M) >= 100:120 try:121 mc = mcs(M, size=PROT["inference"]["mcs"]["alpha"], reps=reps, block=PROT["inference"]["mcs"]["block_length"])122 res[f"mcs_{tier}"] = {r.model: {"p": float(r.pvalue), "in": bool(r.included)} for r in mc.itertuples()}123 except Exception as e: # noqa: BLE001124 res[f"mcs_{tier}"] = {"error": repr(e)[:200]}125 if b in M and len(cand) > 1:126 try:127 res[f"spa_{tier}"] = spa(M[b].values, M[[c for c in cand if c != b]], reps=reps, block=PROT["inference"]["spa"]["block_length"])128 except Exception as e: # noqa: BLE001129 res[f"spa_{tier}"] = {"error": repr(e)[:200]}130 # exploratory: MCS on one representative specification per architecture (Student-t where available, else the131 # single available law) to reduce the model count and sharpen the set (not pre-registered; labelled as such)132 rep = {}133 for sp in specs:134 mid, dist, _ = sp.split("__")135 if res["specs"][sp]["coverage"] < MIN_COVERAGE_MCS:136 continue137 pref = {"t": 0, "normal": 1, "empirical": 2, "ged": 3, "skewt": 4}[dist]138 if mid not in rep or pref < rep[mid][0]:139 rep[mid] = (pref, sp)140 for tier in ("A", "B"):141 cand = [v[1] for v in rep.values() if spec_meta(v[1])["tier"] == tier]142 if tier == "B" and b not in cand:143 cand.append(b)144 M = L[cand].dropna() if len(cand) >= 2 else None145 if M is not None and len(M) >= 100:146 try:147 mc = mcs(M, size=PROT["inference"]["mcs"]["alpha"], reps=reps, block=PROT["inference"]["mcs"]["block_length"])148 res[f"mcs_rep_{tier}"] = {r.model: {"p": float(r.pvalue), "in": bool(r.included)} for r in mc.itertuples()}149 res[f"mcs_rep_{tier}_n"] = int(len(M))150 except Exception as e: # noqa: BLE001151 res[f"mcs_rep_{tier}"] = {"error": repr(e)[:200]}152 out["horizons"][h] = res153 # regimes (pairwise ratios + DM within regime)154 R = regimes_for(panel, L.index, PROT["periods"]["validation"]["start"], PROT["periods"]["validation"]["end"])155 reg = {}156 for col in ("vol_tercile", "vol_tercile66", "prev_sign", "episode"):157 for lab in sorted(R[col].unique()):158 idx = R.index[R[col] == lab]159 Lr = L.loc[idx]160 if len(Lr) < 60:161 continue162 d = {}163 for s in specs:164 pair = np.isfinite(Lr[s].values) & np.isfinite(Lr[b].values)165 if pair.sum() >= 30:166 dm = dm_test(Lr[s].values[pair], Lr[b].values[pair], h)167 d[s] = {"n": int(pair.sum()), "ratio": float(Lr[s].values[pair].mean() / Lr[b].values[pair].mean()), "dm_p": dm["p"], "dm_stat": dm["stat"]}168 reg[f"{col}={lab}"] = d169 out["regimes"][h] = reg170 # secondary losses / proxies at every horizon (pairwise ratios only)171 b = bench_id(specs)172 for (kind, proxy, loss) in [("loss", PRIMARY_PROXY, ln) for ln in PROT["secondary_losses"]] + [("loss", px, PRIMARY_LOSS) for px in PROT["secondary_proxies"]]:173 for h in HORIZONS:174 L = block(W, kind, proxy, loss, h); L = L[(L.index >= lo) & (L.index <= hi)]175 if L.empty or b not in L:176 continue177 d = {}178 for s in L.columns:179 pair = np.isfinite(L[s].values) & np.isfinite(L[b].values)180 if pair.sum() >= 30:181 d[s] = float(L[s].values[pair].mean() / L[b].values[pair].mean()) if loss != "mse" or True else np.nan182 out["secondary"][f"{proxy}|{loss}|{h}"] = d183 # tail risk, h = 1184 F1 = block(W, "forecast", "", "", 1); F1 = F1[(F1.index >= lo) & (F1.index <= hi)]185 MU = block(W, "mu", "", "", 1).reindex(F1.index); D1 = block(W, "dp1", "", "", 1).reindex(F1.index); D2 = block(W, "dp2", "", "", 1).reindex(F1.index)186 p = panel.set_index(pd.to_datetime(panel["session"]))["r"]187 r_next = p.shift(-1).reindex(F1.index).values # realized return of the day after the origin188 tail = {}189 for s in F1.columns:190 meta = spec_meta(s)191 dist = meta["dist"]192 h1 = F1[s].values; mu = MU[s].values if s in MU else np.zeros(len(F1))193 ok = np.isfinite(h1) & np.isfinite(r_next)194 if ok.sum() < 100:195 continue196 emp = None197 if dist == "empirical" and len(G):198 g = G[G["spec_id"] == s].sort_values("origin")199 if "emp_q01" in g:200 gi = g.set_index(pd.to_datetime(g["origin"]))201 qq = {}202 for lvl, k in ((0.01, "01"), (0.05, "05")):203 qq[lvl] = {"q": gi[f"emp_q{k}"].reindex(F1.index, method="ffill").values, "es": gi[f"emp_es{k}"].reindex(F1.index, method="ffill").values}204 emp = qq205 try:206 fc = var_es_forecasts(np.nan_to_num(mu), np.where(ok, h1, np.nan), dist, np.nan_to_num(D1[s].values if s in D1 else np.full(len(F1), np.nan), nan=5.0),207 np.nan_to_num(D2[s].values if s in D2 else np.full(len(F1), np.nan)), emp_q=emp)208 tail[s] = {str(lvl): backtest_table(np.where(ok, r_next, np.nan), fc, lvl) for lvl in (0.01, 0.05)}209 except Exception as e: # noqa: BLE001210 tail[s] = {"error": repr(e)[:200]}211 out["tailrisk"] = tail212 # failure and compute accounting (refits in period)213 acc = {}214 if len(G):215 G["origin"] = pd.to_datetime(G["origin"])216 Gp = G[(G["origin"] >= lo) & (G["origin"] <= hi)]217 for s, g in Gp.groupby("spec_id"):218 acc[s] = {"n_refits": int(len(g)), "failed": float((g["status"] == "failed").mean()), "inadmissible": float((g["status"] == "inadmissible").mean()),219 "fit_time_mean": float(g["time_s"].mean()), "fit_time_median": float(g["time_s"].median()), "iter_mean": float(g["n_iter"].mean()),220 "restarts_mean": float(g["n_restarts"].mean())}221 Sp = S[(S.index >= lo) & (S.index <= hi)]222 for s in Sp.columns:223 acc.setdefault(s, {})224 vc = Sp[s].value_counts(normalize=True).to_dict()225 acc[s].update({"status_" + k: float(v) for k, v in vc.items()})226 out["accounting"] = acc227 return out228229230def pooled(per_sym: dict, run: str, period: str, reps: int, outdir: Path) -> None:231 """Pooled ratios (equal weight per instrument) with a joint block bootstrap over dates; counts of significant DM wins/losses232 (Holm across instruments per (spec, h)); MCS inclusion counts."""233 rows = []234 syms = sorted(per_sym)235 specs = sorted({s for sym in syms for s in per_sym[sym]["horizons"][HORIZONS[0]]["specs"]})236 # per-origin loss matrices are needed for the bootstrap: reload the primary blocks237 cache = {}238 for sym in syms:239 W, T, S, G = load_symbol(run, sym)240 per = PROT["periods"][period]241 for h in HORIZONS:242 L = block(W, "loss", PRIMARY_PROXY, PRIMARY_LOSS, h)243 cache[(sym, h)] = L[(L.index >= per["start"]) & (L.index <= per["end"])]244 for h in HORIZONS:245 b = bench_id(specs)246 for s in specs:247 for cls in ["all"] + sorted(set(CLASS_OF.values())):248 members = [sym for sym in syms if cls == "all" or CLASS_OF.get(sym) == cls]249 ratios, dms, n = [], [], 0250 D = {}251 for sym in members:252 r = per_sym[sym]["horizons"][h]["specs"].get(s, {})253 if np.isfinite(r.get("ratio", np.nan)):254 ratios.append(r["ratio"]); dms.append(r["dm_p"]); n += 1255 L = cache[(sym, h)]256 if s in L:257 pair = np.isfinite(L[s].values) & np.isfinite(L[b].values)258 D[sym] = pd.Series((L[s].values - L[b].values) / np.nanmean(L[b].values[pair]), index=L.index).where(pair)259 if n == 0:260 continue261 dmp = np.array([p for p in dms if np.isfinite(p)])262 # Holm across instruments263 m = len(dmp); order = np.argsort(dmp); adj = np.minimum.accumulate((dmp[order] * (m - np.arange(m)))[::-1])[::-1] if m else np.array([])264 holm = np.empty(m); holm[order] = np.minimum(adj, 1.0) if m else holm265 wins = sum(1 for sym, p_ in zip([x for x in members if np.isfinite(per_sym[x]["horizons"][h]["specs"].get(s, {}).get("ratio", np.nan))], holm)266 if p_ < 0.10 and per_sym[sym]["horizons"][h]["specs"][s]["ratio"] < 1)267 losses_ = sum(1 for sym, p_ in zip([x for x in members if np.isfinite(per_sym[x]["horizons"][h]["specs"].get(s, {}).get("ratio", np.nan))], holm)268 if p_ < 0.10 and per_sym[sym]["horizons"][h]["specs"][s]["ratio"] > 1)269 boot = {"mean": np.nan, "ci_low": np.nan, "ci_high": np.nan, "p": np.nan}270 if D and s != b:271 Dm = pd.DataFrame(D)272 Dm = Dm.dropna(how="all")273 if len(Dm) > 100:274 boot = block_bootstrap_mean_ci(Dm, reps=min(reps, 1000), block=22)275 mcs_in = sum(1 for sym in members for tier in ("A", "B") if per_sym[sym]["horizons"][h].get(f"mcs_{tier}", {}).get(s, {}).get("in"))276 rows.append({"h": h, "spec_id": s, **spec_meta(s), "asset_class": cls, "n_instruments": n, "mean_ratio": float(np.mean(ratios)),277 "median_ratio": float(np.median(ratios)), "share_ratio_below_1": float(np.mean(np.array(ratios) < 1)),278 "dm_wins_holm10": wins, "dm_losses_holm10": losses_, "boot_mean_rel_diff": boot["mean"], "boot_ci_low": boot["ci_low"],279 "boot_ci_high": boot["ci_high"], "boot_p": boot["p"], "mcs_inclusions": mcs_in})280 P = pd.DataFrame(rows)281 P.to_csv(outdir / "pooled.csv", index=False)282 return P283284285def main() -> None:286 ap = argparse.ArgumentParser()287 ap.add_argument("--run", default="main")288 ap.add_argument("--period", default="test", choices=["test", "validation"])289 ap.add_argument("--symbols", default="")290 ap.add_argument("--reps", type=int, default=2000)291 a = ap.parse_args()292 d = ROOT / "results" / "metrics" / f"losses_{a.run}"293 syms = a.symbols.split(",") if a.symbols else sorted(Path(f).stem for f in glob.glob(str(d / "*.parquet")) if not Path(f).stem.startswith(("targets_", "status_", "fits_")))294 outdir = ROOT / "results" / "metrics" / f"analysis_{a.run}_{a.period}"; (outdir / "per_symbol").mkdir(parents=True, exist_ok=True)295 per_sym = {}296 for sym in syms:297 res = analyze_symbol(a.run, sym, a.period, a.reps)298 (outdir / "per_symbol" / f"{sym}.json").write_text(json.dumps(res, indent=0, default=float))299 per_sym[sym] = res300 h1 = res["horizons"][1]301 best = min(((s, v["ratio"]) for s, v in h1["specs"].items() if np.isfinite(v.get("ratio", np.nan))), key=lambda x: x[1])302 print(f"{sym:7s} n={h1['n_origins']} valid={h1['n_valid_targets']} best@h1={best[0]} ratio={best[1]:.3f} mcsA={len([k for k,v in h1.get('mcs_A',{}).items() if isinstance(v,dict) and v.get('in')])}", flush=True)303 P = pooled(per_sym, a.run, a.period, a.reps, outdir)304 tabs = ROOT / "results" / "tables"; tabs.mkdir(exist_ok=True)305 P.to_csv(tabs / f"pooled_{a.run}_{a.period}.csv", index=False)306 print("pooled table", outdir / "pooled.csv", len(P))307308309if __name__ == "__main__":310 main()311