1"""Event study (H2): factor changes in announcement windows, variance shares, event/placebo ratios, cumulative response paths,2pre-announcement drift, by event type and regime. Also yield-space rotation: short rate y(0)=L+S, long rate L, implied Δy(τ)."""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.ns import loadings, LAMBDA_DL8CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics"9rng = np.random.default_rng(47)10Fm = pd.read_parquet(CACHE / "factors_1m.parquet").set_index("ts")11Fm.index = pd.to_datetime(Fm.index, utc=True).tz_convert("America/New_York")12Fm = Fm.dropna(subset=["dL"])13Fm["dSR"] = Fm["dL"] + Fm["dS"] # short-rate change (yield at tau -> 0)14L = loadings(np.array([0.25, 2.0, 5.0, 10.0, 30.0]), LAMBDA_DL)15for i, t in enumerate(["y3m", "y2", "y5", "y10", "y30"]):16 Fm[f"d{t}"] = Fm["dL"] * L[i, 0] + Fm["dS"] * L[i, 1] + Fm["dC"] * L[i, 2]17COLS = ["dL", "dS", "dC", "dSR", "dy3m", "dy2", "dy5", "dy10", "dy30"]18daily = pd.read_parquet(CACHE / "factors_daily_realized.parquet")19ev = pd.read_csv(OUT / "events.csv", parse_dates=["date"]); ev["ts"] = pd.to_datetime(ev["ts"], utc=True).dt.tz_convert("America/New_York")20ev = ev[(ev.ts >= Fm.index.min()) & (ev.ts <= Fm.index.max())].reset_index(drop=True)21PRE, POST = 5, 2522def window_sum(ts, pre=PRE, post=POST, cols=COLS):23 s = Fm.loc[ts - pd.Timedelta(minutes=pre): ts + pd.Timedelta(minutes=post) - pd.Timedelta(minutes=1), cols]24 if len(s) < (pre + post) * 0.8: return None25 return s.sum()26def regime(d):27 d = pd.Timestamp(d)28 if d <= pd.Timestamp("2015-12-15"): return "2011_2015"29 if d <= pd.Timestamp("2019-12-31"): return "2016_2019"30 if d <= pd.Timestamp("2021-12-31"): return "2020_2021"31 return "2022_2026"32rows = []33for _, e in ev.iterrows():34 ws = window_sum(e.ts)35 if ws is None: continue36 yday = Fm.loc[e.ts - pd.Timedelta(minutes=PRE): e.ts, "yday"]37 yd = pd.Timestamp(yday.iloc[-1]) if len(yday) else None38 rec = {"date": e.date, "type": e.type, "scheduled": e.scheduled, "ts": e.ts, "regime": regime(e.date), "yday": yd}39 rec.update({f"w_{c}": float(ws[c]) for c in COLS})40 for pre_, post_, tag in ((5, 15, "w15"), (5, 90, "w90")):41 ws2 = window_sum(e.ts, pre_, post_, cols=["dL", "dS", "dC", "dSR", "dy2", "dy10"])42 if ws2 is not None: rec.update({f"{tag}_{c}": float(ws2[c]) for c in ws2.index})43 if yd is not None and yd in daily.index:44 for c, k in (("dL", "L"), ("dS", "S"), ("dC", "C")):45 rv = daily.loc[yd, f"RV_{k}_1m"]46 rec[f"share_{c}"] = float(ws[c] ** 2 / rv) if rv > 0 else np.nan47 # share of window RV (sum of squares in window) in the day's RV48 sq = (Fm.loc[e.ts - pd.Timedelta(minutes=PRE): e.ts + pd.Timedelta(minutes=POST) - pd.Timedelta(minutes=1), c] ** 2).sum()49 rec[f"rvshare_{c}"] = float(sq / rv) if rv > 0 else np.nan50 # pre-announcement drift: cumulative change from 24h before to -5 min51 pre = Fm.loc[e.ts - pd.Timedelta(hours=24): e.ts - pd.Timedelta(minutes=PRE) - pd.Timedelta(minutes=1), ["dL", "dS", "dC", "dSR", "dy10"]].sum()52 rec.update({f"pre24h_{c}": float(pre[c]) for c in pre.index})53 post = Fm.loc[e.ts + pd.Timedelta(minutes=POST): e.ts + pd.Timedelta(hours=24), ["dL", "dS", "dC", "dSR", "dy10"]].sum()54 rec.update({f"post24h_{c}": float(post[c]) for c in post.index})55 rows.append(rec)56E = pd.DataFrame(rows); E.to_csv(OUT / "event_windows.csv", index=False)57# placebo: same clock time on non-event business days (sample 3 per event, same regime)58evdays = set(pd.to_datetime(ev.date).dt.normalize())59alldays = sorted(set(pd.to_datetime(Fm["yday"]).dt.normalize()) - evdays)60alldays = pd.DatetimeIndex(alldays)61prow = []62for _, e in E.iterrows():63 cands = alldays[(alldays >= pd.Timestamp(e.date) - pd.Timedelta(days=200)) & (alldays <= pd.Timestamp(e.date) + pd.Timedelta(days=200))]64 cands = [c for c in cands if c.weekday() == pd.Timestamp(e.date).weekday()] or list(cands)65 for c in rng.choice(cands, size=min(3, len(cands)), replace=False):66 ts = pd.Timestamp(c).tz_localize("America/New_York") + (e.ts - e.ts.normalize())67 ws = window_sum(ts, cols=["dL", "dS", "dC", "dSR", "dy2", "dy10"])68 if ws is None: continue69 prow.append({"type": e.type, "regime": e.regime, **{f"w_{k}": float(v) for k, v in ws.items()}})70P = pd.DataFrame(prow); P.to_csv(OUT / "event_windows_placebo.csv", index=False)71def boot_ratio(a, b, n=2000):72 """ratio of SDs with bootstrap CI (iid over events / placebo windows)."""73 a, b = np.asarray(a), np.asarray(b); a = a[~np.isnan(a)]; b = b[~np.isnan(b)]74 r = a.std() / b.std()75 bs = [rng.choice(a, len(a)).std() / rng.choice(b, len(b)).std() for _ in range(n)]76 return float(r), [float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))]77def boot_mean(x, n=2000):78 x = np.asarray(x, float); x = x[~np.isnan(x)]79 bs = [rng.choice(x, len(x)).mean() for _ in range(n)]80 return float(x.mean()), [float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))]81summary = {"window_minutes": [-PRE, POST], "benchmark_share_30min": 30 / 1380}82for typ in ["FOMC", "CPI", "EMP"]:83 for reg in ["all", "2011_2015", "2016_2019", "2020_2021", "2022_2026"]:84 e = E[(E.type == typ) & ((E.regime == reg) if reg != "all" else True) & (E.scheduled)]85 p = P[(P.type == typ) & ((P.regime == reg) if reg != "all" else True)]86 if len(e) < 8: continue87 rec = {"n_events": int(len(e)), "n_placebo": int(len(p))}88 for c in ["dL", "dS", "dC", "dSR", "dy2", "dy10"]:89 rec[f"sd_{c}_bp"] = round(float(e[f"w_{c}"].std() * 100), 2)90 rec[f"mean_abs_{c}_bp"] = round(float(e[f"w_{c}"].abs().mean() * 100), 2)91 if len(p) > 8:92 r, ci = boot_ratio(e[f"w_{c}"], p[f"w_{c}"]); rec[f"ratio_sd_vs_placebo_{c}"] = [round(r, 2), [round(ci[0], 2), round(ci[1], 2)]]93 for c in ["dL", "dS", "dC"]:94 m, ci = boot_mean(e[f"rvshare_{c}"]); rec[f"rvshare_{c}"] = [round(m, 4), [round(ci[0], 4), round(ci[1], 4)]]95 for tag in ("w15", "w90"):96 for c in ["dL", "dS", "dC", "dSR", "dy2", "dy10"]:97 if f"{tag}_{c}" in e: rec[f"sd_{tag}_{c}_bp"] = round(float(e[f"{tag}_{c}"].std() * 100), 2)98 for c in ["dL", "dS", "dC"]:99 rec[f"rvshare_{c}_median"] = round(float(e[f"rvshare_{c}"].median()), 4)100 # factor mix within the window: variance decomposition of the 2y/10y/30y implied moves? simpler: share of trace of window covariance of (dL,dS,dC) in bp^2101 V = (e[["w_dL", "w_dS", "w_dC"]] * 100).var()102 rec["var_mix_LSC"] = (V / V.sum()).round(3).to_dict()103 # yield-space mix: variance of implied 3m, 2y, 10y, 30y window moves (bp^2)104 rec["var_dy_bp2"] = {t: round(float((e[f"w_dy{t}"] * 100).var()), 2) for t in ["3m", "2", "5", "10", "30"]}105 # pre-drift106 for c in ["dL", "dS", "dSR", "dy10"]:107 m, ci = boot_mean(e[f"pre24h_{c}"] * 100); rec[f"pre24h_mean_{c}_bp"] = [round(m, 2), [round(ci[0], 2), round(ci[1], 2)]]108 summary[f"{typ}_{reg}"] = rec109# cumulative response path: mean squared cumulative change at k minutes after (and before) the release, normalised by its value at +25110paths = {}111grid = list(range(-60, 121))112for typ in ["FOMC", "CPI", "EMP"]:113 e = E[(E.type == typ) & (E.scheduled)]114 acc = {c: [] for c in ["dL", "dS", "dC", "dSR", "dy10"]}115 for _, r in e.iterrows():116 seg = Fm.loc[r.ts - pd.Timedelta(minutes=60): r.ts + pd.Timedelta(minutes=120), ["dL", "dS", "dC", "dSR", "dy10"]]117 if len(seg) < 150: continue118 seg = seg.reindex(pd.date_range(r.ts - pd.Timedelta(minutes=60), r.ts + pd.Timedelta(minutes=120), freq="1min"), fill_value=0.0)119 cum = seg.cumsum()120 for c in acc: acc[c].append(cum[c].values)121 paths[typ] = {c: {"rms_bp": (np.sqrt(np.mean(np.array(v) ** 2, axis=0)) * 100).round(3).tolist(), "n": len(v)} for c, v in acc.items() if len(v)}122 paths[typ]["grid_minutes"] = grid123json.dump(summary, open(OUT / "event_summary.json", "w"), indent=1)124json.dump(paths, open(OUT / "event_paths.json", "w"))125for k in ["FOMC_all", "CPI_all", "EMP_all", "FOMC_2011_2015", "FOMC_2022_2026", "CPI_2011_2015", "CPI_2022_2026"]:126 if k in summary:127 s = summary[k]; print(k, {kk: s[kk] for kk in s if kk.startswith(("n_", "sd_", "rvshare_dL", "rvshare_dS", "rvshare_dC", "var_mix", "ratio_sd_vs_placebo_dSR", "ratio_sd_vs_placebo_dy10")) and not kk.endswith("median")})128