1"""Robustness computations quoted in Section 'Robustness': (a) expanding-window weights; (b) yield-day boundary at 15:00 and 16:00 ET;2(c) event-window RV shares with 5-minute denominators; (d) FOMC statement minute = minute of maximum volume (12:00–15:00 ET) vs the rule."""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"]10F = pd.read_csv(OUT / "daily_factors_cmt.csv", index_col=0, parse_dates=True).dropna(); dF = F.diff().dropna().loc["2008":]11R = pd.read_parquet(CACHE / "returns_1m.parquet").set_index("ts"); R.index = pd.to_datetime(R.index, utc=True).tz_convert("America/New_York")12def daily_from_boundary(hour, minute):13 b = pd.DatetimeIndex(F.index).tz_localize("America/New_York") + pd.Timedelta(hours=hour, minutes=minute)14 pos = np.searchsorted(b.values, (R.index + pd.Timedelta(minutes=1)).values, side="left"); pos = np.clip(pos, 0, len(b) - 1)15 yd = pd.Series(F.index.values[pos], index=R.index)16 g = R[ROOTS].groupby(yd.values); D = g.sum(min_count=1); N = g.count(); N.columns = [f"n_{c}" for c in ROOTS]17 return D.join(N)18def rolling(D, expanding=False):19 X = D[ROOTS].where(D[[f"n_{r}" for r in ROOTS]].ge(100).values); panel = dF.join(X, how="inner"); preds = []20 for q in pd.period_range("2011Q1", panel.index[-1].to_period("Q"), freq="Q"):21 start = panel.index[0] if expanding else q.start_time - pd.DateOffset(years=3)22 train = panel.loc[start: q.start_time - pd.Timedelta(days=1)]; test = panel.loc[q.start_time:q.end_time]23 roots = [r for r in ROOTS if train[r].notna().mean() >= 0.85]24 if expanding: # with an expanding window, presence is judged on the last 3 years25 roots = [r for r in ROOTS if panel.loc[q.start_time - pd.DateOffset(years=3): q.start_time - pd.Timedelta(days=1), r].notna().mean() >= 0.85]26 tr = train.dropna(subset=roots); te = test.dropna(subset=roots)27 if len(tr) < 500 or len(te) == 0: continue28 Xa = np.column_stack([np.ones(len(tr)), tr[roots].values])29 for f in FAC:30 b = np.linalg.lstsq(Xa, tr[f].values, rcond=None)[0]31 preds.append(pd.DataFrame({"date": te.index, "factor": f, "actual": te[f].values, "pred": te[roots].values @ b[1:]}))32 P = pd.concat(preds)33 return {f: round(float(1 - ((P[P.factor == f].actual - P[P.factor == f].pred) ** 2).sum() / ((P[P.factor == f].actual - P[P.factor == f].actual.mean()) ** 2).sum()), 4) for f in FAC}34res = {}35D1530 = pd.read_parquet(CACHE / "returns_daily.parquet"); D1530.index = pd.to_datetime(D1530.index)36res["rolling_3y_1530"] = rolling(D1530)37res["expanding_1530"] = rolling(D1530, expanding=True)38for h, mnt, lab in ((15, 0, "1500"), (16, 0, "1600")):39 res[f"rolling_3y_{lab}"] = rolling(daily_from_boundary(h, mnt))40# (c) RV shares with 5-minute denominators41E = pd.read_csv(OUT / "event_windows.csv", parse_dates=["date", "yday"]); daily = pd.read_parquet(CACHE / "factors_daily_realized.parquet"); daily.index = pd.to_datetime(daily.index)42E = E[E.scheduled]; E = E.join(daily[["RV_L_5m", "RV_S_5m", "RV_C_5m"]], on="yday")43res["rvshare_5m_denominator"] = {}44for typ in ("FOMC", "CPI", "EMP"):45 e = E[E.type == typ]46 res["rvshare_5m_denominator"][typ] = {c: round(float((e[f"w_{c}"] ** 2 / e[f"RV_{k}_5m"]).clip(upper=1).mean()), 4) for c, k in (("dL", "L"), ("dS", "S"), ("dC", "C"))}47 res["rvshare_5m_denominator"][typ + "_1m_single_sum"] = {c: round(float(e[f"share_{c}"].clip(upper=1).mean()), 4) for c in ("dL", "dS", "dC")}48# (d) FOMC release minute vs volume49ev = pd.read_csv(OUT / "events.csv", parse_dates=["date"]); ev["ts"] = pd.to_datetime(ev["ts"], utc=True).dt.tz_convert("America/New_York")50fo = ev[(ev.type == "FOMC") & (ev.scheduled) & (ev.date >= "2011-01-01")]51# activity measure: sum of absolute 1-minute returns of ZF, ZN, US (volume peaks at the 15:00 settlement, so the search window ends at 14:45)52V = R[["ZF", "ZN", "US"]].abs().sum(axis=1)53diffs = []54for _, e in fo.iterrows():55 day = e.ts.normalize(); w = V.loc[day + pd.Timedelta(hours=12): day + pd.Timedelta(hours=14, minutes=45)]56 if len(w) == 0: continue57 diffs.append((w.idxmax() - e.ts).total_seconds() / 60)58diffs = np.array(diffs)59res["fomc_time_check"] = {"n": int(len(diffs)), "share_within_0_2_min": round(float(np.mean((diffs >= 0) & (diffs <= 2))), 3), "share_within_0_5_min": round(float(np.mean((diffs >= 0) & (diffs <= 5))), 3),60 "median_offset_min": float(np.median(diffs)), "offsets_outside_0_5": sorted([float(d) for d in diffs if not (0 <= d <= 5)])}61json.dump(res, open(OUT / "robustness_checks.json", "w"), indent=1); print(json.dumps(res, indent=1))62