1"""Export small aggregated series for the paper's figures (results/figdata/). Runs on M3U96a. Never writes raw bars."""2import sys3from pathlib import Path4import numpy as np, pandas as pd5ROOT = Path(__file__).resolve().parents[1]; sys.path.insert(0, str(ROOT / "src"))6from wp47.ns import loadings, LAMBDA_DL7CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"; FIG = ROOT / "results/figdata"; FIG.mkdir(parents=True, exist_ok=True)8ET = "America/New_York"9Fm = pd.read_parquet(CACHE / "factors_1m.parquet").set_index("ts")10Fm.index = pd.to_datetime(Fm.index, utc=True).tz_convert(ET)11Fm = Fm.dropna(subset=["dL"])12Fm["yday"] = pd.to_datetime(Fm["yday"]).dt.normalize()13# 1. sample days: cumulative changes over the yield day + implied 2y/10y14L2, L10 = loadings(np.array([2.0]), LAMBDA_DL)[0], loadings(np.array([10.0]), LAMBDA_DL)[0]15DAYS = {"2013-06-19": "FOMC", "2022-06-15": "FOMC", "2022-11-10": "CPI", "2024-08-02": "EMP"}16parts = []17for d, typ in DAYS.items():18 s = Fm[Fm["yday"] == pd.Timestamp(d)][["dL", "dS", "dC"]].copy()19 if len(s) == 0:20 print("missing day", d); continue21 c = s.cumsum() * 10022 c.columns = ["cum_L_bp", "cum_S_bp", "cum_C_bp"]23 c["cum_y2_bp"] = (s["dL"] * L2[0] + s["dS"] * L2[1] + s["dC"] * L2[2]).cumsum() * 10024 c["cum_y10_bp"] = (s["dL"] * L10[0] + s["dS"] * L10[1] + s["dC"] * L10[2]).cumsum() * 10025 c["day"] = d; c["type"] = typ26 c.index.name = "ts"27 parts.append(c.reset_index())28pd.concat(parts).to_csv(FIG / "sample_days.csv", index=False)29# 4. fidelity by year for CMT (primary) and GSW target30rows = []31for lab, f in (("cmt", "mimicking_oos_daily.parquet"), ("gsw", "mimicking_oos_daily_gsw.parquet")):32 P = pd.read_parquet(CACHE / f)33 pcol = "pred_oos_nocst" if "pred_oos_nocst" in P.columns else "pred"34 P["year"] = pd.to_datetime(P["date"]).dt.year35 for (y, fac), g in P.groupby(["year", "factor"]):36 if len(g) < 30: continue37 r2 = 1 - np.sum((g.actual - g[pcol]) ** 2) / np.sum((g.actual - g.actual.mean()) ** 2)38 rows.append({"target": lab, "year": int(y), "factor": fac, "r2_oos": float(r2), "n": int(len(g))})39pd.DataFrame(rows).to_csv(FIG / "fidelity_by_year.csv", index=False)40# 7. ETF vs futures FOMC-window pairs41E = pd.read_parquet(CACHE / "factors_1m_etf.parquet")42E.index = pd.to_datetime(E.index, utc=True).tz_convert(ET) if E.index.tz is None or str(E.index.tz) != ET else E.index43J = E.join(Fm[["dL", "dS", "dC"]], how="inner", lsuffix="_etf", rsuffix="_fut")44ev = pd.read_csv(OUT / "events.csv", parse_dates=["date"]); ev["ts"] = pd.to_datetime(ev["ts"], utc=True).dt.tz_convert(ET)45rows = []46for _, e in ev[(ev.type == "FOMC") & (ev.scheduled) & (ev.ts >= J.index.min())].iterrows():47 s = J.loc[e.ts - pd.Timedelta(minutes=5): e.ts + pd.Timedelta(minutes=24)]48 if len(s) >= 24:49 w = s.sum() * 10050 rows.append({"date": e.date.strftime("%Y-%m-%d"), **{k: float(v) for k, v in w.items()}})51pd.DataFrame(rows).to_csv(FIG / "etf_fomc_pairs.csv", index=False)52print("exported", {f.name: f.stat().st_size for f in FIG.glob("*.csv")})5354# ------------------------------------------------------------------------------------------------55# 8. Event heatmaps: cumulative implied 2-year yield change (bp) from -5 min, columns -30..+120 minutes, rows = events56rng = np.random.default_rng(47)57L2v = loadings(np.array([2.0]), LAMBDA_DL)[0]58Fm["dy2"] = Fm["dL"] * L2v[0] + Fm["dS"] * L2v[1] + Fm["dC"] * L2v[2]59grid = np.arange(-30, 121)60for typ in ("FOMC", "CPI", "EMP"):61 rows = []62 for _, e in ev[(ev.type == typ) & (ev.scheduled) & (ev.date >= "2011-01-01")].iterrows():63 seg = Fm.loc[e.ts - pd.Timedelta(minutes=30): e.ts + pd.Timedelta(minutes=120), ["dy2", "dL", "dS", "dC"]]64 if len(seg) < 120: continue65 idx = pd.date_range(e.ts - pd.Timedelta(minutes=30), e.ts + pd.Timedelta(minutes=120), freq="1min")66 seg = seg.reindex(idx, fill_value=0.0)67 cum = seg["dy2"].cumsum().values * 10068 base = cum[25] # value at -5 min69 rows.append([e.date.strftime("%Y-%m-%d")] + list(np.round(cum - base, 3)))70 pd.DataFrame(rows, columns=["date"] + [str(m) for m in grid]).to_csv(FIG / f"event_heatmap_{typ}.csv", index=False)71# 9. Event paths with bootstrap bands (RMS across events; 300 resamples of events)72out = []73for typ in ("FOMC", "CPI", "EMP"):74 mats = {c: [] for c in ("dL", "dS", "dC", "dSR", "dy10")}75 L10v = loadings(np.array([10.0]), LAMBDA_DL)[0]76 for _, e in ev[(ev.type == typ) & (ev.scheduled) & (ev.date >= "2011-01-01")].iterrows():77 seg = Fm.loc[e.ts - pd.Timedelta(minutes=60): e.ts + pd.Timedelta(minutes=120), ["dL", "dS", "dC"]]78 if len(seg) < 150: continue79 idx = pd.date_range(e.ts - pd.Timedelta(minutes=60), e.ts + pd.Timedelta(minutes=120), freq="1min"); seg = seg.reindex(idx, fill_value=0.0)80 seg["dSR"] = seg["dL"] + seg["dS"]; seg["dy10"] = seg["dL"] * L10v[0] + seg["dS"] * L10v[1] + seg["dC"] * L10v[2]81 cum = seg.cumsum() * 10082 for c in mats: mats[c].append(cum[c].values)83 for c, v in mats.items():84 A = np.array(v); n = len(A)85 rms = np.sqrt((A ** 2).mean(0))86 bs = np.array([np.sqrt((A[rng.integers(0, n, n)] ** 2).mean(0)) for _ in range(300)])87 lo, hi = np.percentile(bs, 5, axis=0), np.percentile(bs, 95, axis=0)88 for i, m in enumerate(range(-60, 121)): out.append({"type": typ, "series": c, "minute": m, "rms": rms[i], "lo": lo[i], "hi": hi[i], "n": n})89pd.DataFrame(out).to_csv(FIG / "event_paths_boot.csv", index=False)90# 10. Fidelity by year with bootstrap 90% bands (resampling days within the year, 500 draws)91rows = []92for lab, f in (("cmt", "mimicking_oos_daily.parquet"), ("gsw", "mimicking_oos_daily_gsw.parquet")):93 P = pd.read_parquet(CACHE / f); pcol = "pred_oos_nocst" if "pred_oos_nocst" in P.columns else "pred"; P["year"] = pd.to_datetime(P["date"]).dt.year94 for (y, fac), g in P.groupby(["year", "factor"]):95 if len(g) < 30: continue96 a, p = g.actual.values, g[pcol].values97 def r2(a, p): return 1 - np.sum((a - p) ** 2) / np.sum((a - a.mean()) ** 2)98 bs = [r2(a[i], p[i]) for i in (rng.integers(0, len(a), len(a)) for _ in range(500))]99 rows.append({"target": lab, "year": int(y), "factor": fac, "r2_oos": float(r2(a, p)), "lo": float(np.percentile(bs, 5)), "hi": float(np.percentile(bs, 95)), "n": int(len(g))})100pd.DataFrame(rows).to_csv(FIG / "fidelity_by_year.csv", index=False)101print("exported", {f.name: f.stat().st_size for f in FIG.glob("*.csv")})102