1"""Factor-mimicking portfolios: rolling OLS of daily NS factor changes on daily futures returns; strictly out-of-sample weights.23Weights for quarter q are estimated on the trailing 3 years ending the last day before q (>= 500 days with all roots), then applied4to q. Also: expanding-window and full-sample (in-sample) weights for comparison. Outputs weights, OOS fitted daily changes, fidelity.5"""6import sys, json7from pathlib import Path8import numpy as np, pandas as pd9ROOT = Path(__file__).resolve().parents[1]; sys.path.insert(0, str(ROOT / "src"))10from wp47.returns import ROOTS11CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"12FAC = ["level", "slope", "curvature"]13F = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).dropna()14dF = F.diff().dropna().loc["2008-01-01":]15D = pd.read_parquet(CACHE / "returns_daily.parquet")16D.index = pd.to_datetime(D.index)17MIN_MIN = 100 # a root counts as observed on a yield day if >= 100 minute returns (roll days of ZT/UB are sparse)18def active_roots(day):19 return [r for r in ROOTS if D.loc[day, f"n_{r}"] >= MIN_MIN] if day in D.index else []20X = D[ROOTS].where(D[[f"n_{r}" for r in ROOTS]].ge(MIN_MIN).values)21panel = dF.join(X, how="inner")22# root sets: TN enters 2016; UB enters 2010. We use the set of roots available for >= 95% of the window days.23def ols(Y, Xm):24 Xa = np.column_stack([np.ones(len(Xm)), Xm])25 b = np.linalg.lstsq(Xa, Y, rcond=None)[0]26 return b27quarters = pd.period_range("2011Q1", panel.index[-1].to_period("Q"), freq="Q")28rows, preds = [], []29for q in quarters:30 q_start, q_end = q.start_time, q.end_time31 train = panel.loc[q_start - pd.DateOffset(years=3): q_start - pd.Timedelta(days=1)]32 test = panel.loc[q_start:q_end]33 if len(test) == 0: continue34 roots = [r for r in ROOTS if train[r].notna().mean() >= 0.85]35 tr = train.dropna(subset=roots)36 if len(tr) < 500: continue37 te = test.dropna(subset=roots)38 W = {}39 for f in FAC:40 b = ols(tr[f].values, tr[roots].values)41 W[f] = dict(zip(["const"] + roots, b.round(6)))42 preds.append(pd.DataFrame({"date": te.index, "factor": f, "actual": te[f].values,43 "pred_oos": te[roots].values @ b[1:] + b[0], "pred_oos_nocst": te[roots].values @ b[1:], "quarter": str(q)}))44 rows.append({"quarter": str(q), "n_train": len(tr), "roots": roots, "weights": W,45 "r2_train": {f: float(1 - np.var(tr[f] - (np.column_stack([np.ones(len(tr)), tr[roots].values]) @ ols(tr[f].values, tr[roots].values))) / np.var(tr[f])) for f in FAC}})46P = pd.concat(preds, ignore_index=True)47P.to_parquet(CACHE / "mimicking_oos_daily.parquet")48json.dump(rows, open(OUT / "mimicking_weights_quarterly.json", "w"), indent=1)49def fidelity(sub):50 out = {}51 for f in FAC:52 s = sub[sub.factor == f]53 if len(s) < 30: continue54 b = np.polyfit(s.pred_oos_nocst, s.actual, 1)55 r2 = 1 - np.sum((s.actual - s.pred_oos_nocst) ** 2) / np.sum((s.actual - s.actual.mean()) ** 2)56 out[f] = {"n": int(len(s)), "r2_oos": round(float(r2), 4), "corr": round(float(np.corrcoef(s.pred_oos_nocst, s.actual)[0, 1]), 4),57 "slope": round(float(b[0]), 4), "rmse_bp": round(float(np.sqrt(np.mean((s.actual - s.pred_oos_nocst) ** 2)) * 100), 2),58 "sd_actual_bp": round(float(s.actual.std() * 100), 2)}59 return out60periods = {"full_2011_2026": ("2011-01-01", "2026-12-31"), "2011_2015_ZLB": ("2011-01-01", "2015-12-15"), "2016_2019": ("2015-12-16", "2019-12-31"),61 "2020_2021": ("2020-01-01", "2021-12-31"), "2022_2026": ("2022-01-01", "2026-12-31")}62fid = {k: fidelity(P[(P.date >= a) & (P.date <= b)]) for k, (a, b) in periods.items()}63# yearly R264fid["by_year"] = {str(y): fidelity(P[P.date.dt.year == y]) for y in sorted(P.date.dt.year.unique())}65# full-sample in-sample weights (upper bound) on 2016–2026 with six roots66full = panel.loc["2016-01-01":].dropna(subset=ROOTS)67fid["insample_2016_2026_six_roots"] = {}68for f in FAC:69 b = ols(full[f].values, full[ROOTS].values); fitv = np.column_stack([np.ones(len(full)), full[ROOTS].values]) @ b70 fid["insample_2016_2026_six_roots"][f] = {"n": len(full), "r2": round(float(1 - np.var(full[f] - fitv) / np.var(full[f])), 4), "weights": dict(zip(["const"] + ROOTS, b.round(4)))}71json.dump(fid, open(OUT / "mimicking_fidelity.json", "w"), indent=1)72print(json.dumps({k: fid[k] for k in periods}, indent=1)); print(json.dumps(fid["insample_2016_2026_six_roots"], indent=1))73