1"""Robustness (pre-specified): fidelity against GSW-based daily factors (1–30y zero curve, no bills).2Extension (post-hoc, disclosed): add the front 3-month SOFR future (SR3) as a seventh instrument from 2021."""3import sys, json4from pathlib import Path5import numpy as np, pandas as pd6ROOT = Path(__file__).resolve().parents[1]; sys.path.insert(0, str(ROOT / "src"))7from wp47.returns import ROOTS8CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"9FAC = ["level", "slope", "curvature"]10D = pd.read_parquet(CACHE / "returns_daily.parquet"); D.index = pd.to_datetime(D.index)11MIN_MIN = 10012def run(dF, roots_all, start_q, label, min_minutes=MIN_MIN, extra=None):13 X = D[roots_all].where(D[[f"n_{r}" for r in roots_all]].ge(min_minutes).values)14 if extra is not None: X = X.join(extra, how="left"); roots_all = roots_all + list(extra.columns)15 panel = dF.join(X, how="inner")16 quarters = pd.period_range(start_q, panel.index[-1].to_period("Q"), freq="Q")17 preds = []18 for q in quarters:19 train = panel.loc[q.start_time - pd.DateOffset(years=3): q.start_time - pd.Timedelta(days=1)]20 test = panel.loc[q.start_time:q.end_time]21 roots = [r for r in roots_all if train[r].notna().mean() >= 0.85]22 tr = train.dropna(subset=roots); te = test.dropna(subset=roots)23 if len(tr) < 500 or len(te) == 0: continue24 Xa = np.column_stack([np.ones(len(tr)), tr[roots].values])25 for f in FAC:26 b = np.linalg.lstsq(Xa, tr[f].values, rcond=None)[0]27 preds.append(pd.DataFrame({"date": te.index, "factor": f, "actual": te[f].values, "pred": te[roots].values @ b[1:]}))28 P = pd.concat(preds)29 def fid(s):30 return {"n": int(len(s)), "r2_oos": round(float(1 - ((s.actual - s.pred) ** 2).sum() / ((s.actual - s.actual.mean()) ** 2).sum()), 4),31 "slope": round(float(np.polyfit(s.pred, s.actual, 1)[0]), 4), "rmse_bp": round(float(np.sqrt(((s.actual - s.pred) ** 2).mean()) * 100), 2)}32 res = {"label": label}33 for per, (a, b) in {"full": ("2011-01-01", "2026-12-31"), "2011_2015": ("2011-01-01", "2015-12-15"), "2016_2019": ("2015-12-16", "2019-12-31"),34 "2020_2021": ("2020-01-01", "2021-12-31"), "2022_2026": ("2022-01-01", "2026-12-31")}.items():35 s = P[(P.date >= a) & (P.date <= b)]36 if len(s): res[per] = {f: fid(s[s.factor == f]) for f in FAC}37 return res, P38out = {}39FG = pd.read_csv(OUT / "daily_factors_gsw.csv", index_col=0, parse_dates=True).diff().dropna().loc["2008":]40out["gsw"], P = run(FG, ROOTS, "2011Q1", "GSW 1-30y zero curve factors")41P.to_parquet(CACHE / "mimicking_oos_daily_gsw.parquet")42FC = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).diff().dropna().loc["2008":]43# CMT factors re-estimated without the 1m/3m/6m bills (maturities >= 1y) to isolate the bill effect44from wp47.ns import read_fred, fit_panel45y = read_fred(ROOT / "data/external/fred").loc["2007-01-01":]46F1, _ = fit_panel(y[[c for c in y.columns if c not in ("DGS1MO", "DGS3MO", "DGS6MO")]])47F1.to_csv(OUT / "daily_factors_cmt_1y_plus.csv")48out["cmt_1y_plus"], _ = run(F1.diff().dropna().loc["2008":], ROOTS, "2011Q1", "CMT factors from maturities >= 1y")49# SR3 extension (daily returns from cached 1m SR3)50sr = pd.read_parquet(CACHE / "sr3_daily.parquet") if (CACHE / "sr3_daily.parquet").exists() else None51if sr is not None:52 out["cmt_with_sr3"], P3 = run(FC, ROOTS, "2021Q1", "CMT factors, six roots + SR3 (from 2021)", extra=sr[["SR3"]])53 P3.to_parquet(CACHE / "mimicking_oos_daily_sr3.parquet")54json.dump(out, open(OUT / "mimicking_fidelity_robustness.json", "w"), indent=1)55for k, v in out.items():56 print(k, {per: {f: vv[f]["r2_oos"] for f in FAC} for per, vv in v.items() if per != "label"})57