Repository Public Research preview
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 out["horizons"][h] = res131 # regimes (pairwise ratios + DM within regime)132 R = regimes_for(panel, L.index, PROT["periods"]["validation"]["start"], PROT["periods"]["validation"]["end"])133 reg = {}134 for col in ("vol_tercile", "vol_tercile66", "prev_sign", "episode"):135 for lab in sorted(R[col].unique()):136 idx = R.index[R[col] == lab]137 Lr = L.loc[idx]138 if len(Lr) < 60:139 continue140 d = {}141 for s in specs:142 pair = np.isfinite(Lr[s].values) & np.isfinite(Lr[b].values)143 if pair.sum() >= 30:144 dm = dm_test(Lr[s].values[pair], Lr[b].values[pair], h)145 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"]}146 reg[f"{col}={lab}"] = d147 out["regimes"][h] = reg148 # secondary losses / proxies at every horizon (pairwise ratios only)149 b = bench_id(specs)150 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"]]:151 for h in HORIZONS:152 L = block(W, kind, proxy, loss, h); L = L[(L.index >= lo) & (L.index <= hi)]153 if L.empty or b not in L:154 continue155 d = {}156 for s in L.columns:157 pair = np.isfinite(L[s].values) & np.isfinite(L[b].values)158 if pair.sum() >= 30:159 d[s] = float(L[s].values[pair].mean() / L[b].values[pair].mean()) if loss != "mse" or True else np.nan160 out["secondary"][f"{proxy}|{loss}|{h}"] = d161 # tail risk, h = 1162 F1 = block(W, "forecast", "", "", 1); F1 = F1[(F1.index >= lo) & (F1.index <= hi)]163 MU = block(W, "mu", "", "", 1).reindex(F1.index); D1 = block(W, "dp1", "", "", 1).reindex(F1.index); D2 = block(W, "dp2", "", "", 1).reindex(F1.index)164 p = panel.set_index(pd.to_datetime(panel["session"]))["r"]165 r_next = p.shift(-1).reindex(F1.index).values # realized return of the day after the origin166 tail = {}167 for s in F1.columns:168 meta = spec_meta(s)169 dist = meta["dist"]170 h1 = F1[s].values; mu = MU[s].values if s in MU else np.zeros(len(F1))171 ok = np.isfinite(h1) & np.isfinite(r_next)172 if ok.sum() < 100:173 continue174 emp = None175 if dist == "empirical" and len(G):176 g = G[G["spec_id"] == s].sort_values("origin")177 if "emp_q01" in g:178 gi = g.set_index(pd.to_datetime(g["origin"]))179 qq = {}180 for lvl, k in ((0.01, "01"), (0.05, "05")):181 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}182 emp = qq183 try:184 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),185 np.nan_to_num(D2[s].values if s in D2 else np.full(len(F1), np.nan)), emp_q=emp)186 tail[s] = {str(lvl): backtest_table(np.where(ok, r_next, np.nan), fc, lvl) for lvl in (0.01, 0.05)}187 except Exception as e: # noqa: BLE001188 tail[s] = {"error": repr(e)[:200]}189 out["tailrisk"] = tail190 # failure and compute accounting (refits in period)191 acc = {}192 if len(G):193 G["origin"] = pd.to_datetime(G["origin"])194 Gp = G[(G["origin"] >= lo) & (G["origin"] <= hi)]195 for s, g in Gp.groupby("spec_id"):196 acc[s] = {"n_refits": int(len(g)), "failed": float((g["status"] == "failed").mean()), "inadmissible": float((g["status"] == "inadmissible").mean()),197 "fit_time_mean": float(g["time_s"].mean()), "fit_time_median": float(g["time_s"].median()), "iter_mean": float(g["n_iter"].mean()),198 "restarts_mean": float(g["n_restarts"].mean())}199 Sp = S[(S.index >= lo) & (S.index <= hi)]200 for s in Sp.columns:201 acc.setdefault(s, {})202 vc = Sp[s].value_counts(normalize=True).to_dict()203 acc[s].update({"status_" + k: float(v) for k, v in vc.items()})204 out["accounting"] = acc205 return out206207208def pooled(per_sym: dict, run: str, period: str, reps: int, outdir: Path) -> None:209 """Pooled ratios (equal weight per instrument) with a joint block bootstrap over dates; counts of significant DM wins/losses210 (Holm across instruments per (spec, h)); MCS inclusion counts."""211 rows = []212 syms = sorted(per_sym)213 specs = sorted({s for sym in syms for s in per_sym[sym]["horizons"][HORIZONS[0]]["specs"]})214 # per-origin loss matrices are needed for the bootstrap: reload the primary blocks215 cache = {}216 for sym in syms:217 W, T, S, G = load_symbol(run, sym)218 per = PROT["periods"][period]219 for h in HORIZONS:220 L = block(W, "loss", PRIMARY_PROXY, PRIMARY_LOSS, h)221 cache[(sym, h)] = L[(L.index >= per["start"]) & (L.index <= per["end"])]222 for h in HORIZONS:223 b = bench_id(specs)224 for s in specs:225 for cls in ["all"] + sorted(set(CLASS_OF.values())):226 members = [sym for sym in syms if cls == "all" or CLASS_OF.get(sym) == cls]227 ratios, dms, n = [], [], 0228 D = {}229 for sym in members:230 r = per_sym[sym]["horizons"][h]["specs"].get(s, {})231 if np.isfinite(r.get("ratio", np.nan)):232 ratios.append(r["ratio"]); dms.append(r["dm_p"]); n += 1233 L = cache[(sym, h)]234 if s in L:235 pair = np.isfinite(L[s].values) & np.isfinite(L[b].values)236 D[sym] = pd.Series((L[s].values - L[b].values) / np.nanmean(L[b].values[pair]), index=L.index).where(pair)237 if n == 0:238 continue239 dmp = np.array([p for p in dms if np.isfinite(p)])240 # Holm across instruments241 m = len(dmp); order = np.argsort(dmp); adj = np.minimum.accumulate((dmp[order] * (m - np.arange(m)))[::-1])[::-1] if m else np.array([])242 holm = np.empty(m); holm[order] = np.minimum(adj, 1.0) if m else holm243 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)244 if p_ < 0.10 and per_sym[sym]["horizons"][h]["specs"][s]["ratio"] < 1)245 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)246 if p_ < 0.10 and per_sym[sym]["horizons"][h]["specs"][s]["ratio"] > 1)247 boot = {"mean": np.nan, "ci_low": np.nan, "ci_high": np.nan, "p": np.nan}248 if D and s != b:249 Dm = pd.DataFrame(D)250 Dm = Dm.dropna(how="all")251 if len(Dm) > 100:252 boot = block_bootstrap_mean_ci(Dm, reps=min(reps, 1000), block=22)253 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"))254 rows.append({"h": h, "spec_id": s, **spec_meta(s), "asset_class": cls, "n_instruments": n, "mean_ratio": float(np.mean(ratios)),255 "median_ratio": float(np.median(ratios)), "share_ratio_below_1": float(np.mean(np.array(ratios) < 1)),256 "dm_wins_holm10": wins, "dm_losses_holm10": losses_, "boot_mean_rel_diff": boot["mean"], "boot_ci_low": boot["ci_low"],257 "boot_ci_high": boot["ci_high"], "boot_p": boot["p"], "mcs_inclusions": mcs_in})258 P = pd.DataFrame(rows)259 P.to_csv(outdir / "pooled.csv", index=False)260 return P261262263def main() -> None:264 ap = argparse.ArgumentParser()265 ap.add_argument("--run", default="main")266 ap.add_argument("--period", default="test", choices=["test", "validation"])267 ap.add_argument("--symbols", default="")268 ap.add_argument("--reps", type=int, default=2000)269 a = ap.parse_args()270 d = ROOT / "results" / "metrics" / f"losses_{a.run}"271 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_")))272 outdir = ROOT / "results" / "metrics" / f"analysis_{a.run}_{a.period}"; (outdir / "per_symbol").mkdir(parents=True, exist_ok=True)273 per_sym = {}274 for sym in syms:275 res = analyze_symbol(a.run, sym, a.period, a.reps)276 (outdir / "per_symbol" / f"{sym}.json").write_text(json.dumps(res, indent=0, default=float))277 per_sym[sym] = res278 h1 = res["horizons"][1]279 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])280 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)281 P = pooled(per_sym, a.run, a.period, a.reps, outdir)282 tabs = ROOT / "results" / "tables"; tabs.mkdir(exist_ok=True)283 P.to_csv(tabs / f"pooled_{a.run}_{a.period}.csv", index=False)284 print("pooled table", outdir / "pooled.csv", len(P))285286287if __name__ == "__main__":288 main()289