1"""P7 inference: pairwise contrasts between saved models (DM-HAC + day-cluster bootstrap, Holm/BH), subgroup2analysis of a contrast, and stress-day calibration (H3). Works from per_origin_*.npz + origins_<tag>.csv.3Usage: pairwise.py --exp main --tag test --pairs mrpt_small:mrpt_small_1m_only,mrpt_small:har_rv,... [--subgroups mrpt_small:har_rv]"""4from __future__ import annotations56import argparse7import json8import sys9from pathlib import Path1011import numpy as np12import pandas as pd1314ROOT = Path(__file__).resolve().parents[1]15sys.path.insert(0, str(ROOT / "src"))16from chronos_x import pipeline as P # noqa: E40217from chronos_x.eval.statistics import holm, benjamini_hochberg, cluster_bootstrap_mean # noqa: E40218from chronos_x.eval.subgroups import add_subgroups # noqa: E40219from chronos_x.eval.calibration import pit_histogram # noqa: E402202122def load(exp: str, name: str, tag: str) -> dict:23 return np.load(ROOT / "results" / "metrics" / exp / f"per_origin_{name}_{tag}.npz")242526def main() -> None:27 ap = argparse.ArgumentParser()28 ap.add_argument("--exp", default="main"); ap.add_argument("--tag", default="test"); ap.add_argument("--pairs", required=True)29 ap.add_argument("--subgroups", default=None); ap.add_argument("--stride", type=int, default=15)30 ap.add_argument("--metrics", default="crps,wis,qlike"); ap.add_argument("--horizons", default="5,15,30,60")31 ap.add_argument("--universe", default="data/manifests/universe_main.json"); ap.add_argument("--out-suffix", default="")32 a = ap.parse_args()33 d = ROOT / "results" / "metrics" / a.exp34 tab = pd.read_csv(d / f"origins_{a.tag}.csv")35 man = json.loads((ROOT / a.universe).read_text())36 meta = add_subgroups(tab, man)37 H = [int(h) for h in a.horizons.split(",")]38 rows = []39 cache = {}40 def per(name):41 if name not in cache:42 cache[name] = load(a.exp, name, a.tag)43 return cache[name]44 for pair in a.pairs.split(","):45 m, ref = pair.split(":")46 for metric in a.metrics.split(","):47 for h in H:48 r = P.compare(per(m)[f"{metric}_{h}"], per(ref)[f"{metric}_{h}"], tab, h, a.stride)49 rows.append({"model": m, "reference": ref, "metric": metric, **r})50 df = pd.DataFrame(rows)51 for metric in a.metrics.split(","):52 sel = df["metric"] == metric53 df.loc[sel, "p_holm"] = holm(df.loc[sel, "dm_p"].tolist()); df.loc[sel, "p_bh"] = benjamini_hochberg(df.loc[sel, "dm_p"].tolist())54 out = d / f"pairwise_{a.tag}{a.out_suffix}.csv"; df.to_csv(out, index=False)55 print(df[df.metric == "crps"][["model", "reference", "H", "rel_diff", "lo", "hi", "dm_t", "dm_p", "p_holm"]].to_string(index=False, float_format=lambda x: f"{x:.4g}"))56 # subgroups for one contrast (relative CRPS difference with day-cluster CI inside each group)57 if a.subgroups:58 m, ref = a.subgroups.split(":")59 srows = []60 for col in ("asset_class", "adv_tercile", "vol_tercile", "group", "tod", "year", "stress_day"):61 for g, idx in meta.groupby(col, observed=True).indices.items():62 for h in H:63 pa_, pb_ = per(m)[f"crps_{h}"][idx], per(ref)[f"crps_{h}"][idx]64 cl = pd.factorize(meta["session_date"].values[idx])[0]65 cb = cluster_bootstrap_mean(pa_ - pb_, cl, n_boot=500)66 base = pb_.mean()67 srows.append({"col": col, "group": str(g), "H": h, "n": int(len(idx)), "rel": cb["mean"] / base, "lo": cb["lo"] / base, "hi": cb["hi"] / base,68 "crps_model": pa_.mean(), "crps_ref": base})69 sdf = pd.DataFrame(srows); sdf.to_csv(d / f"subgroups_{m}_vs_{ref}_{a.tag}.csv", index=False)70 print(sdf[sdf.H == 60].to_string(index=False, float_format=lambda x: f"{x:.4g}"))71 # stress-day calibration (H3): coverage and PIT dispersion on stress vs normal days for every model in the pairs72 names = sorted({x for pair in a.pairs.split(",") for x in pair.split(":")})73 crows = []74 for name in names:75 p = per(name)76 for h in H:77 for lab, idx in (("stress", np.flatnonzero(meta["stress_day"].values)), ("normal", np.flatnonzero(~meta["stress_day"].values))):78 u = p[f"pit_{h}"][idx]79 ph = pit_histogram(u)80 crows.append({"model": name, "H": h, "days": lab, "n": int(len(idx)), "cov_50": float(((u >= .25) & (u <= .75)).mean()),81 "cov_80": float(((u >= .1) & (u <= .9)).mean()), "cov_95": float(((u >= .025) & (u <= .975)).mean()),82 "pit_var": ph["var"], "pit_ks": ph["ks_stat"], "crps": float(p[f"crps_{h}"][idx].mean()),83 "crps_scaled": float(p[f"crps_scaled_{h}"][idx].mean())})84 cdf = pd.DataFrame(crows); cdf.to_csv(d / f"stress_calibration_{a.tag}{a.out_suffix}.csv", index=False)85 print(cdf[cdf.H == 60].to_string(index=False, float_format=lambda x: f"{x:.4g}"))868788if __name__ == "__main__":89 main()90