"""Event study (H2): factor changes in announcement windows, variance shares, event/placebo ratios, cumulative response paths, pre-announcement drift, by event type and regime. Also yield-space rotation: short rate y(0)=L+S, long rate L, implied Δy(τ).""" import sys, json from pathlib import Path import numpy as np, pandas as pd ROOT = Path(__file__).resolve().parents[1]; sys.path.insert(0, str(ROOT / "src")) from wp47.ns import loadings, LAMBDA_DL CACHE = ROOT / "data/cache"; OUT = ROOT / "results/metrics" rng = np.random.default_rng(47) Fm = pd.read_parquet(CACHE / "factors_1m.parquet").set_index("ts") Fm.index = pd.to_datetime(Fm.index, utc=True).tz_convert("America/New_York") Fm = Fm.dropna(subset=["dL"]) Fm["dSR"] = Fm["dL"] + Fm["dS"] # short-rate change (yield at tau -> 0) L = loadings(np.array([0.25, 2.0, 5.0, 10.0, 30.0]), LAMBDA_DL) for i, t in enumerate(["y3m", "y2", "y5", "y10", "y30"]): Fm[f"d{t}"] = Fm["dL"] * L[i, 0] + Fm["dS"] * L[i, 1] + Fm["dC"] * L[i, 2] COLS = ["dL", "dS", "dC", "dSR", "dy3m", "dy2", "dy5", "dy10", "dy30"] daily = pd.read_parquet(CACHE / "factors_daily_realized.parquet") ev = pd.read_csv(OUT / "events.csv", parse_dates=["date"]); ev["ts"] = pd.to_datetime(ev["ts"], utc=True).dt.tz_convert("America/New_York") ev = ev[(ev.ts >= Fm.index.min()) & (ev.ts <= Fm.index.max())].reset_index(drop=True) PRE, POST = 5, 25 def window_sum(ts, pre=PRE, post=POST, cols=COLS): s = Fm.loc[ts - pd.Timedelta(minutes=pre): ts + pd.Timedelta(minutes=post) - pd.Timedelta(minutes=1), cols] if len(s) < (pre + post) * 0.8: return None return s.sum() def regime(d): d = pd.Timestamp(d) if d <= pd.Timestamp("2015-12-15"): return "2011_2015" if d <= pd.Timestamp("2019-12-31"): return "2016_2019" if d <= pd.Timestamp("2021-12-31"): return "2020_2021" return "2022_2026" rows = [] for _, e in ev.iterrows(): ws = window_sum(e.ts) if ws is None: continue yday = Fm.loc[e.ts - pd.Timedelta(minutes=PRE): e.ts, "yday"] yd = pd.Timestamp(yday.iloc[-1]) if len(yday) else None rec = {"date": e.date, "type": e.type, "scheduled": e.scheduled, "ts": e.ts, "regime": regime(e.date), "yday": yd} rec.update({f"w_{c}": float(ws[c]) for c in COLS}) for pre_, post_, tag in ((5, 15, "w15"), (5, 90, "w90")): ws2 = window_sum(e.ts, pre_, post_, cols=["dL", "dS", "dC", "dSR", "dy2", "dy10"]) if ws2 is not None: rec.update({f"{tag}_{c}": float(ws2[c]) for c in ws2.index}) if yd is not None and yd in daily.index: for c, k in (("dL", "L"), ("dS", "S"), ("dC", "C")): rv = daily.loc[yd, f"RV_{k}_1m"] rec[f"share_{c}"] = float(ws[c] ** 2 / rv) if rv > 0 else np.nan # share of window RV (sum of squares in window) in the day's RV sq = (Fm.loc[e.ts - pd.Timedelta(minutes=PRE): e.ts + pd.Timedelta(minutes=POST) - pd.Timedelta(minutes=1), c] ** 2).sum() rec[f"rvshare_{c}"] = float(sq / rv) if rv > 0 else np.nan # pre-announcement drift: cumulative change from 24h before to -5 min 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() rec.update({f"pre24h_{c}": float(pre[c]) for c in pre.index}) post = Fm.loc[e.ts + pd.Timedelta(minutes=POST): e.ts + pd.Timedelta(hours=24), ["dL", "dS", "dC", "dSR", "dy10"]].sum() rec.update({f"post24h_{c}": float(post[c]) for c in post.index}) rows.append(rec) E = pd.DataFrame(rows); E.to_csv(OUT / "event_windows.csv", index=False) # placebo: same clock time on non-event business days (sample 3 per event, same regime) evdays = set(pd.to_datetime(ev.date).dt.normalize()) alldays = sorted(set(pd.to_datetime(Fm["yday"]).dt.normalize()) - evdays) alldays = pd.DatetimeIndex(alldays) prow = [] for _, e in E.iterrows(): cands = alldays[(alldays >= pd.Timestamp(e.date) - pd.Timedelta(days=200)) & (alldays <= pd.Timestamp(e.date) + pd.Timedelta(days=200))] cands = [c for c in cands if c.weekday() == pd.Timestamp(e.date).weekday()] or list(cands) for c in rng.choice(cands, size=min(3, len(cands)), replace=False): ts = pd.Timestamp(c).tz_localize("America/New_York") + (e.ts - e.ts.normalize()) ws = window_sum(ts, cols=["dL", "dS", "dC", "dSR", "dy2", "dy10"]) if ws is None: continue prow.append({"type": e.type, "regime": e.regime, **{f"w_{k}": float(v) for k, v in ws.items()}}) P = pd.DataFrame(prow); P.to_csv(OUT / "event_windows_placebo.csv", index=False) def boot_ratio(a, b, n=2000): """ratio of SDs with bootstrap CI (iid over events / placebo windows).""" a, b = np.asarray(a), np.asarray(b); a = a[~np.isnan(a)]; b = b[~np.isnan(b)] r = a.std() / b.std() bs = [rng.choice(a, len(a)).std() / rng.choice(b, len(b)).std() for _ in range(n)] return float(r), [float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))] def boot_mean(x, n=2000): x = np.asarray(x, float); x = x[~np.isnan(x)] bs = [rng.choice(x, len(x)).mean() for _ in range(n)] return float(x.mean()), [float(np.percentile(bs, 2.5)), float(np.percentile(bs, 97.5))] summary = {"window_minutes": [-PRE, POST], "benchmark_share_30min": 30 / 1380} for typ in ["FOMC", "CPI", "EMP"]: for reg in ["all", "2011_2015", "2016_2019", "2020_2021", "2022_2026"]: e = E[(E.type == typ) & ((E.regime == reg) if reg != "all" else True) & (E.scheduled)] p = P[(P.type == typ) & ((P.regime == reg) if reg != "all" else True)] if len(e) < 8: continue rec = {"n_events": int(len(e)), "n_placebo": int(len(p))} for c in ["dL", "dS", "dC", "dSR", "dy2", "dy10"]: rec[f"sd_{c}_bp"] = round(float(e[f"w_{c}"].std() * 100), 2) rec[f"mean_abs_{c}_bp"] = round(float(e[f"w_{c}"].abs().mean() * 100), 2) if len(p) > 8: 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)]] for c in ["dL", "dS", "dC"]: m, ci = boot_mean(e[f"rvshare_{c}"]); rec[f"rvshare_{c}"] = [round(m, 4), [round(ci[0], 4), round(ci[1], 4)]] for tag in ("w15", "w90"): for c in ["dL", "dS", "dC", "dSR", "dy2", "dy10"]: if f"{tag}_{c}" in e: rec[f"sd_{tag}_{c}_bp"] = round(float(e[f"{tag}_{c}"].std() * 100), 2) for c in ["dL", "dS", "dC"]: rec[f"rvshare_{c}_median"] = round(float(e[f"rvshare_{c}"].median()), 4) # 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^2 V = (e[["w_dL", "w_dS", "w_dC"]] * 100).var() rec["var_mix_LSC"] = (V / V.sum()).round(3).to_dict() # yield-space mix: variance of implied 3m, 2y, 10y, 30y window moves (bp^2) rec["var_dy_bp2"] = {t: round(float((e[f"w_dy{t}"] * 100).var()), 2) for t in ["3m", "2", "5", "10", "30"]} # pre-drift for c in ["dL", "dS", "dSR", "dy10"]: 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)]] summary[f"{typ}_{reg}"] = rec # cumulative response path: mean squared cumulative change at k minutes after (and before) the release, normalised by its value at +25 paths = {} grid = list(range(-60, 121)) for typ in ["FOMC", "CPI", "EMP"]: e = E[(E.type == typ) & (E.scheduled)] acc = {c: [] for c in ["dL", "dS", "dC", "dSR", "dy10"]} for _, r in e.iterrows(): seg = Fm.loc[r.ts - pd.Timedelta(minutes=60): r.ts + pd.Timedelta(minutes=120), ["dL", "dS", "dC", "dSR", "dy10"]] if len(seg) < 150: continue seg = seg.reindex(pd.date_range(r.ts - pd.Timedelta(minutes=60), r.ts + pd.Timedelta(minutes=120), freq="1min"), fill_value=0.0) cum = seg.cumsum() for c in acc: acc[c].append(cum[c].values) 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)} paths[typ]["grid_minutes"] = grid json.dump(summary, open(OUT / "event_summary.json", "w"), indent=1) json.dump(paths, open(OUT / "event_paths.json", "w")) for k in ["FOMC_all", "CPI_all", "EMP_all", "FOMC_2011_2015", "FOMC_2022_2026", "CPI_2011_2015", "CPI_2022_2026"]: if k in summary: 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")})