Repository Public Research preview
1"""Apply out-of-sample quarterly mimicking weights to one-minute futures returns -> minute-level NS factor changes (2011–2026).2Outputs: data/cache/factors_1m.parquet (ts, yday, dL, dS, dC, n_active), realized measures per yield day at 1/5/30-minute sampling,3intraday periodicity by minute-of-day, realized correlations by period (H3, H4 inputs)."""4import sys, json5from pathlib import Path6import numpy as np, pandas as pd7ROOT = Path(__file__).resolve().parents[1]; sys.path.insert(0, str(ROOT / "src"))8from wp47.returns import ROOTS9CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"10FAC = ["level", "slope", "curvature"]; SHORT = {"level": "dL", "slope": "dS", "curvature": "dC"}11W = json.load(open(OUT / "mimicking_weights_quarterly.json"))12R = pd.read_parquet(CACHE / "returns_1m.parquet").set_index("ts")13R.index = pd.to_datetime(R.index, utc=True).tz_convert("America/New_York")14R["yday"] = pd.to_datetime(R["yday"])15# active flag per root: a root is 'active' at minute m if it has a bar within the previous 15 minutes (else NaN, not zero)16act = {}17for r in ROOTS:18 has = R[r].notna()19 last = pd.Series(np.where(has, np.arange(len(R)), np.nan), index=R.index).ffill()20 act[r] = (np.arange(len(R)) - last.values) <= 1521A = pd.DataFrame(act, index=R.index)22Rz = R[ROOTS].fillna(0.0).where(A) # zero return when active but no trade; NaN when inactive23parts = []24for w in W:25 q = pd.Period(w["quarter"], freq="Q"); roots = w["roots"]26 sl = Rz.loc[q.start_time.tz_localize("America/New_York"): q.end_time.tz_localize("America/New_York")]27 if len(sl) == 0: continue28 ok = sl[roots].notna().all(axis=1)29 X = sl[roots].values30 out = pd.DataFrame(index=sl.index)31 for f in FAC:32 wv = np.array([w["weights"][f][r] for r in roots])33 out[SHORT[f]] = np.where(ok, X @ wv, np.nan)34 out["yday"] = R.loc[sl.index, "yday"].values35 out["n_roots"] = len(roots)36 parts.append(out)37Fm = pd.concat(parts)38Fm.index.name = "ts"39Fm.reset_index().to_parquet(CACHE / "factors_1m.parquet", index=False)40print("minutes with factors:", int(Fm[["dL"]].notna().sum().iloc[0]), "of", len(Fm), "first", Fm.index.min(), "last", Fm.index.max())41# realized measures per yield day42def rv_block(df, k):43 g = df.groupby("yday")44 if k == 1:45 s = g[["dL", "dS", "dC"]].agg(lambda x: np.nansum(x ** 2))46 cov = g.apply(lambda x: pd.Series({"cLS": np.nansum(x.dL * x.dS), "cLC": np.nansum(x.dL * x.dC), "cSC": np.nansum(x.dS * x.dC)}))47 else:48 # k-minute sampling: sum within consecutive k-minute bins of the yield day49 tmp = df.copy(); tmp["bin"] = tmp.groupby("yday").cumcount() // k50 agg = tmp.groupby(["yday", "bin"])[["dL", "dS", "dC"]].sum(min_count=1)51 s = agg.groupby("yday").agg(lambda x: np.nansum(x ** 2))52 cov = agg.groupby("yday").apply(lambda x: pd.Series({"cLS": np.nansum(x.dL * x.dS), "cLC": np.nansum(x.dL * x.dC), "cSC": np.nansum(x.dS * x.dC)}))53 s.columns = [f"RV_{c[1]}" for c in s.columns]54 return s.join(cov)55daily = None56for k in (1, 5, 30):57 b = rv_block(Fm.dropna(subset=["dL"]), k); b.columns = [f"{c}_{k}m" for c in b.columns]58 daily = b if daily is None else daily.join(b)59n_min = Fm.dropna(subset=["dL"]).groupby("yday").size().rename("n_minutes")60cum = Fm.groupby("yday")[["dL", "dS", "dC"]].sum(min_count=1); cum.columns = ["cum_dL", "cum_dS", "cum_dC"]61daily = daily.join(n_min).join(cum)62daily.to_parquet(CACHE / "factors_daily_realized.parquet"); daily.to_csv(OUT / "factors_daily_realized.csv")63# realized correlations by period and sampling64def corr(d, k):65 return {"LS": float(d[f"cLS_{k}m"].sum() / np.sqrt(d[f"RV_L_{k}m"].sum() * d[f"RV_S_{k}m"].sum())),66 "LC": float(d[f"cLC_{k}m"].sum() / np.sqrt(d[f"RV_L_{k}m"].sum() * d[f"RV_C_{k}m"].sum())),67 "SC": float(d[f"cSC_{k}m"].sum() / np.sqrt(d[f"RV_S_{k}m"].sum() * d[f"RV_C_{k}m"].sum()))}68periods = {"full": ("2011-01-01", "2026-12-31"), "2011_2015": ("2011-01-01", "2015-12-15"), "2016_2019": ("2015-12-16", "2019-12-31"), "2020_2021": ("2020-01-01", "2021-12-31"), "2022_2026": ("2022-01-01", "2026-12-31")}69H4 = {per: {f"{k}m": corr(daily.loc[a:b], k) for k in (1, 5, 30)} for per, (a, b) in periods.items()}70# daily (FRED) correlations for comparison71dF = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).diff().dropna()72H4["daily_fred"] = {per: {"LS": float(dF.loc[a:b].corr().loc["level", "slope"]), "LC": float(dF.loc[a:b].corr().loc["level", "curvature"]), "SC": float(dF.loc[a:b].corr().loc["slope", "curvature"])} for per, (a, b) in periods.items()}73# volatility signature: average RV by sampling74H4["rv_mean_bp2"] = {f"{k}m": {c: float(daily[f"RV_{c}_{k}m"].mean() * 1e4) for c in "LSC"} for k in (1, 5, 30)}75json.dump(H4, open(OUT / "realized_correlations.json", "w"), indent=1)76# fidelity of the CUMULATED INTRADAY series itself (all yield days with factors), overall and by coverage77dFc = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).diff().dropna()78test = daily[["cum_dL", "cum_dS", "cum_dC", "n_minutes"]].join(dFc, how="inner")79def r2(a, p): return float(1 - ((a - p) ** 2).sum() / ((a - a.mean()) ** 2).sum())80fidc = {}81for lab, sel in (("all_days", test.n_minutes > 0), ("ge600_minutes", test.n_minutes >= 600), ("lt600_minutes", test.n_minutes < 600)):82 t = test[sel]83 fidc[lab] = {"n": int(len(t)), **{f: {"r2": round(r2(t[f], t[c]), 4), "slope": round(float(np.polyfit(t[c], t[f], 1)[0]), 4)} for f, c in (("level", "cum_dL"), ("slope", "cum_dS"), ("curvature", "cum_dC"))}}84for per, (a, b) in periods.items():85 t = test.loc[a:b]; fidc[f"all_days_{per}"] = {"n": int(len(t)), **{f: {"r2": round(r2(t[f], t[c]), 4)} for f, c in (("level", "cum_dL"), ("slope", "cum_dS"), ("curvature", "cum_dC"))}}86fidc["coverage"] = {"median_minutes": float(test.n_minutes.median()), "share_days_ge600": float((test.n_minutes >= 600).mean()), "share_days_ge1200": float((test.n_minutes >= 1200).mean())}87json.dump(fidc, open(OUT / "fidelity_cumulated_intraday.json", "w"), indent=1); print("cumulated-intraday fidelity:", {k: v for k, v in fidc.items() if k in ("all_days", "ge600_minutes", "lt600_minutes", "coverage")})88# intraday periodicity: mean squared factor change by minute-of-day (ET), by period; also announcement-day vs non-announcement89Fm["tod"] = Fm.index.hour * 60 + Fm.index.minute90ev = pd.read_csv(OUT / "events.csv", parse_dates=["date"])91evdays = set(pd.to_datetime(ev.date).dt.normalize())92Fm["evday"] = pd.to_datetime(Fm["yday"]).dt.normalize().isin(evdays).values93def tod_profile(d):94 g = d.groupby("tod")95 p = g[["dL", "dS", "dC"]].agg(lambda x: np.nanmean(x ** 2)) * 1e4 # bp^2 per minute96 p["n"] = g.size(); return p97prof = {"all": tod_profile(Fm), "non_event_days": tod_profile(Fm[~Fm.evday]), "event_days": tod_profile(Fm[Fm.evday])}98for per, (a, b) in periods.items():99 if per != "full": prof[per] = tod_profile(Fm.loc[a:b])100P = pd.concat(prof, axis=1); P.to_csv(OUT / "intraday_periodicity.csv")101print(json.dumps({k: v for k, v in H4.items() if k in ("full", "2011_2015", "2022_2026", "daily_fred")}, indent=1))102