TeX 53.6%
Python 46.2%
1"""Stage F.1 — idiosyncratic shock detection and propagation measurement.23Per month m, on the LAGGED top-500 universe (RTH, flags excluded):45 Residualization (all inputs lagged, CLAUDE.md §8):6 beta_i and residual scale MAD_i estimated on month m-1; minute-of-day7 seasonal factor f_tau = cross-sectional median |residual| profile of m-1,8 normalized to mean 1. In month m:9 e_it = r_it - beta_i * r_mkt(-i),t (leave-one-out EW market)10 z_it = e_it / (1.4826 * MAD_i * f_tau)1112 Shock (baseline |z| >= 6, sign recorded; 4 and 5 stored for robustness):13 isolation: |z_mkt,t| < 1 AND share of |z_jt| > 4 at t <= 2%.1415 Outcomes for each event at h in {1,2,5,10,15,30,60,120} minutes, same session16 only (censored otherwise):17 breadth(h) = share of eligible j != i with |z_{j,t+h}| > 218 intensity(h) = mean |z_{j,t+h}|19 signed(h) = mean sign(e_i) * z_{j,t+h}20 pre-shock falsification: same outcomes at t-5.2122 Matched placebo per event: same ticker, same minute-of-day, random other day23 in month with valid quiet z (|z|<2) and calm market — same outcome code path.2425Output: outputs/estimates/events/events_YYYY-MM.parquet26Usage: python scripts/11_... [--workers N] [--month YYYY-MM]27"""2829from __future__ import annotations3031import argparse32from concurrent.futures import ProcessPoolExecutor, as_completed33from pathlib import Path3435import numpy as np36import polars as pl3738ROOT = Path(__file__).resolve().parents[1]39MINUTE = ROOT / "data" / "processed" / "minute"40UNIV = ROOT / "data" / "metadata" / "universe_monthly.parquet"41OUT = ROOT / "outputs" / "estimates" / "events"4243HORIZONS = [1, 2, 5, 10, 15, 30, 60, 120]44Z_SHOCK = 6.045Z_RECEIVER = 2.046Z_MKT_CALM = 1.047SIMUL_CAP = 0.0248PRE_H = 549SEED = 20260713505152def month_pair(month: str) -> str:53 y, m = int(month[:4]), int(month[5:7])54 m -= 155 if m == 0:56 y, m = y - 1, 1257 return f"{y:04d}-{m:02d}"585960def load_matrix(month: str, tickers: list[str]):61 """(returns matrix [rows x N], day-block slices, mods, dates) on the RTH grid."""62 y, mm = month.split("-")63 df = (64 pl.scan_parquet(MINUTE / f"year={y}" / f"month={mm}" / "data.parquet")65 .filter((pl.col("session") == "rth") & ~pl.col("flag_ohlc")66 & ~pl.col("flag_split_window") & pl.col("ticker").is_in(tickers))67 .select("ticker", "date", "mod", "close")68 .collect()69 )70 if df.height == 0:71 return None72 grid = df.pivot(on="ticker", index=["date", "mod"], values="close").sort("date", "mod")73 cols = [c for c in grid.columns if c not in ("date", "mod")]74 px = grid.select(cols).to_numpy()75 dates = grid["date"].to_numpy()76 mods = grid["mod"].to_numpy()77 logpx = np.log(px)78 ret = np.vstack([np.full((1, logpx.shape[1]), np.nan),79 logpx[1:] - logpx[:-1]])80 newday = np.empty(len(dates), dtype=bool)81 newday[0] = True82 newday[1:] = dates[1:] != dates[:-1]83 ret[newday] = np.nan # never bridge overnight84 # day block bounds85 starts = np.flatnonzero(newday)86 ends = np.append(starts[1:], len(dates))87 return ret, cols, dates, mods, list(zip(starts, ends))888990def residualize(ret: np.ndarray, beta: np.ndarray):91 """z-scores with LOO EW market; returns (e, mkt_raw)."""92 valid = ~np.isnan(ret)93 s = np.nansum(ret, axis=1)94 n = valid.sum(axis=1).astype(float)95 with np.errstate(invalid="ignore", divide="ignore"):96 mkt = s / n97 loo = (s[:, None] - np.where(valid, ret, 0.0)) / np.maximum(n[:, None] - valid, 1.0)98 e = ret - beta[None, :] * loo99 return e, mkt, n100101102def month_params(prev_month: str, tickers: list[str]):103 """beta, MAD, seasonal factor from the PRIOR month."""104 loaded = load_matrix(prev_month, tickers)105 if loaded is None:106 return None107 ret, cols, dates, mods, _ = loaded108 valid = ~np.isnan(ret)109 s = np.nansum(ret, axis=1)110 n = valid.sum(axis=1).astype(float)111 with np.errstate(invalid="ignore", divide="ignore"):112 mkt = s / n113 ok = np.isfinite(mkt) & (n >= 50)114 r, m = ret[ok], mkt[ok]115 mv = np.nanvar(m)116 beta = np.array([117 (np.nanmean((r[:, j] - np.nanmean(r[:, j])) * (m - m.mean())) / mv118 if np.isfinite(r[:, j]).sum() > 300 else 1.0)119 for j in range(r.shape[1])])120 beta = np.clip(np.nan_to_num(beta, nan=1.0), -3, 5)121 e = r - beta[None, :] * m[:, None]122 med = np.nanmedian(e, axis=0)123 mad = np.nanmedian(np.abs(e - med[None, :]), axis=0)124 mad = np.where((mad > 0) & np.isfinite(mad), mad, np.nan)125 # seasonal factor: cross-sectional median |e| per minute-of-day, normalized126 mods_ok = mods[ok]127 absmed = np.abs(e / (1.4826 * mad[None, :]))128 ftau = {}129 for mod in np.unique(mods_ok):130 ftau[int(mod)] = np.nanmedian(absmed[mods_ok == mod])131 mean_f = np.nanmean(list(ftau.values()))132 ftau = {k: (v / mean_f if np.isfinite(v) and v > 0 else 1.0)133 for k, v in ftau.items()}134 mkt_scale = 1.4826 * np.nanmedian(np.abs(m - np.nanmedian(m)))135 return dict(zip(cols, beta)), dict(zip(cols, mad)), ftau, mkt_scale136137138def outcomes_at(z: np.ndarray, row: int, j: int, d0: int, d1: int):139 zr = z[row]140 mask = np.isfinite(zr)141 mask[j] = False142 n = int(mask.sum())143 if n < 50 or not (d0 <= row < d1):144 return None145 az = np.abs(zr[mask])146 return (float((az > Z_RECEIVER).mean()), float(az.mean()),147 float(zr[mask].mean()), n)148149150def process_month(month: str, force: bool = False) -> str:151 dest = OUT / f"events_{month}.parquet"152 if dest.exists() and not force:153 return f"{month} skip"154 prev = month_pair(month)155 uni = (pl.read_parquet(UNIV)156 .filter((pl.col("month") == month) & pl.col("top500"))["ticker"].to_list())157 if not uni:158 dest.parent.mkdir(parents=True, exist_ok=True)159 pl.DataFrame([]).write_parquet(dest)160 return f"{month} no-universe"161 params = month_params(prev, uni)162 loaded = load_matrix(month, uni)163 if params is None or loaded is None:164 dest.parent.mkdir(parents=True, exist_ok=True)165 pl.DataFrame([]).write_parquet(dest)166 return f"{month} no-data"167 beta_d, mad_d, ftau, mkt_scale = params168 ret, cols, dates, mods, blocks = loaded169 beta = np.array([beta_d.get(c, 1.0) for c in cols])170 mad = np.array([mad_d.get(c, np.nan) for c in cols])171 e, mkt, n_t = residualize(ret, beta)172 f = np.array([ftau.get(int(m), 1.0) for m in mods])173 with np.errstate(invalid="ignore", divide="ignore"):174 z = e / (1.4826 * mad[None, :] * f[:, None])175 z_mkt = mkt / (mkt_scale * f)176177 absz = np.abs(z)178 simul = np.nanmean(absz > 4.0, axis=1)179 row_ok = (np.abs(z_mkt) < Z_MKT_CALM) & (simul <= SIMUL_CAP) & (n_t >= 50)180181 # day block id per row182 block_of = np.empty(len(mods), dtype=int)183 bounds = []184 for b, (s0, s1) in enumerate(blocks):185 block_of[s0:s1] = b186 bounds.append((s0, s1))187188 ev_rows, ev_cols = np.where((absz >= Z_SHOCK) & row_ok[:, None])189 rng = np.random.default_rng(SEED + hash(month) % 100000)190 records = []191192 def build_record(row, j, kind):193 d0, d1 = bounds[block_of[row]]194 rec = {195 "month": month, "kind": kind, "ticker": cols[j],196 "date": str(dates[row]), "mod": int(mods[row]),197 "z": float(z[row, j]), "sign": int(np.sign(z[row, j])),198 "mkt_z": float(z_mkt[row]),199 "mkt_vol30": float(np.nanstd(mkt[max(row - 30, d0):row]) * 1e4),200 }201 pre = outcomes_at(z, row - PRE_H, j, d0, d1) if row - PRE_H >= d0 else None202 rec["pre_breadth"], rec["pre_intensity"] = (203 (pre[0], pre[1]) if pre else (None, None))204 got = False205 for h in HORIZONS:206 o = outcomes_at(z, row + h, j, d0, d1) if row + h < d1 else None207 rec[f"breadth_{h}"] = o[0] if o else None208 rec[f"intensity_{h}"] = o[1] if o else None209 rec[f"signed_{h}"] = (o[2] * rec["sign"]) if o else None210 rec[f"nelig_{h}"] = o[3] if o else None211 got = got or (o is not None)212 return rec if got else None213214 for row, j in zip(ev_rows, ev_cols):215 rec = build_record(row, j, "real")216 if rec is None:217 continue218 records.append(rec)219 # matched placebo: same ticker & minute-of-day, other calm day220 cand = np.flatnonzero(221 (mods == rec["mod"]) & row_ok & np.isfinite(z[:, j])222 & (np.abs(z[:, j]) < 2.0) & (block_of != block_of[row]))223 if len(cand):224 prow = int(rng.choice(cand))225 prec = build_record(prow, j, "placebo")226 if prec is not None:227 prec["z"] = rec["z"] # carry source magnitude for matching228 prec["sign"] = rec["sign"]229 records.append(prec)230231 dest.parent.mkdir(parents=True, exist_ok=True)232 pl.DataFrame(records).write_parquet(dest)233 nre = sum(1 for r in records if r["kind"] == "real")234 return f"{month} events={nre} placebos={len(records)-nre}"235236237def main() -> None:238 ap = argparse.ArgumentParser()239 ap.add_argument("--workers", type=int, default=4)240 ap.add_argument("--month", default=None)241 ap.add_argument("--force", action="store_true")242 args = ap.parse_args()243 months = sorted(p.stem.replace("daily_", "") for p in244 (ROOT / "data" / "metadata" / "daily").glob("daily_*.parquet"))245 months = [m for m in months if m > "1992-01"] # first month has no lag246 if args.month:247 print(process_month(args.month, True))248 return249 with ProcessPoolExecutor(max_workers=args.workers) as ex:250 futs = {ex.submit(process_month, m, args.force): m for m in months}251 for i, fut in enumerate(as_completed(futs), 1):252 print(f"[{i}/{len(months)}] {fut.result()}", flush=True)253254255if __name__ == "__main__":256 main()257