1"""r3k_17_event_study_v2.py — Stage D (FIRC-0002, package L5, cluster F: ECON-08, EDIT-06, ADVE-02, AI_C-08, THEO-05, AI_C-09, DATA-07).2ChatGPT (2022-11-30) event study rebuilt (v2 of r3k_05_event_study.py):3 window anchored at the launch: day 0 = 2022-12-01 (first trading day after the 30 Nov 2022 release); the pre-launch block4 [-20, -1] (1-29 Nov 2022, inside the author's window) is reported separately; horizons [-1,+1], [0,+5], [0,+20], 60, 125, 250, 395 days;5 terciles on rank(method='first') (removes the pd.qcut tie fragility found at L0); sub-samples full / largest-1,000 / remaining6 (corrected cap, and the author's v1 cap for the before/after) / S&P 500 / non-S&P, plus terciles formed WITHIN size group;7 benchmarks (a) IWV market model, est. 2022-01-03..2022-10-31 (author) · (b) beta-only AR = r - b m · (c) market-adjusted r - m ·8 (d) market model estimated over calendar 2021 · (e) FF5+MOM daily calendar-time alpha of the H-L portfolio over the event window9 (HAC t) · (f) size x momentum matched (5 x 5 cells from the 2021-cohort panel; AR = r - mean of the cell excluding the firm);10 inference for every cell: H-L daily-portfolio CDA t (estimation-window SD of the H-L portfolio), event-window-SD t, and the11 BMP (Boehmer-Musumeci-Poulsen 1991) standardized-CAR t for the H-L difference (independence-based, reported for reference), KP12 factor stated once; placebo launches 30 Nov 2015..2021 (same design, exposure = the preceding cohort) with the 2022 rank.13Inputs: data/processed_r3k_v2/{prices_daily_v2, benchmarks_daily_v2, ff_daily, analysis_panel_v2}.parquet (+ staged14prices_daily_event.parquet only to reproduce the author's baseline). Writes results/r3k_v2/results_event.json, t_event_cda.csv,15event_car_v2.parquet, event_car_paths_v2.parquet. Run: OMP_NUM_THREADS=4 nice -n 10 .venv-review/bin/python -I scripts/r3k_17_event_study_v2.py16"""17import os, sys, json, time, warnings, numpy as np, pandas as pd, statsmodels.api as sm, statsmodels.formula.api as smf18sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))); import r3k_inference as I19from scipy import stats20warnings.filterwarnings("ignore"); T0 = time.time()21V = "data/processed_r3k_v2"; T = "results/r3k_v2"22panel = pd.read_parquet(f"{V}/analysis_panel_v2.parquet")23dp = pd.read_parquet(f"{V}/prices_daily_v2.parquet"); dp["date"] = pd.to_datetime(dp["date"]); dp = dp.sort_values(["ticker", "date"]); dp["ret"] = dp.groupby("ticker")["adj_close"].pct_change()24dp.loc[(dp.ret < -0.99) | (dp.adj_close <= 0) | ~np.isfinite(dp.ret), "ret"] = np.nan25bm = pd.read_parquet(f"{V}/benchmarks_daily_v2.parquet"); bm["date"] = pd.to_datetime(bm["date"]); bm = bm.sort_values(["ticker", "date"]); bm["ret"] = bm.groupby("ticker")["adj_close"].pct_change()26mkt_all = bm[bm.ticker == "IWV"][["date", "ret"]].rename(columns={"ret": "mkt"}).set_index("date")["mkt"]27ffd = pd.read_parquet(f"{V}/ff_daily.parquet"); ffd["date"] = pd.to_datetime(ffd["date"]); ffd = ffd.set_index("date")28RES = {"script": "r3k_17_event_study_v2.py"}29HORIZONS = [("car_m1_p1", -1, 1), ("car_0_5", 0, 5), ("car_0_20", 0, 20), ("car_60d", 0, 59), ("car_6m", 0, 124), ("car_1y", 0, 249), ("car_18m", 0, 394)]3031def exposure(cohort_year):32 e = panel[panel.fiscal_year == cohort_year].copy()33 e["ai_pre"] = (e.ai_adoption_wmean.transform(lambda g: (g - g.mean()) / g.std()) + e.ai_infrastructure_wmean.transform(lambda g: (g - g.mean()) / g.std())) / 234 e = e[["ticker", "ai_pre", "log_mktcap", "mktcap", "mktcap_v1", "mom_12m", "sector", "size_group", "size_group_v1", "size_group_pct42", "in_sp500"]].dropna(subset=["ai_pre"]).drop_duplicates("ticker")35 e["tercile"] = pd.qcut(e.ai_pre.rank(method="first"), 3, labels=["Low", "Mid", "High"]); e["tercile_within_size"] = e.groupby("size_group", dropna=False)["ai_pre"].transform(lambda s: pd.qcut(s.rank(method="first"), 3, labels=["Low", "Mid", "High"]) if len(s) >= 30 else pd.Series([None] * len(s), index=s.index))36 e["size5"] = pd.qcut(e.log_mktcap.rank(method="first"), 5, labels=False); e["mom5"] = pd.qcut(e.mom_12m.rank(method="first"), 5, labels=False); return e3738def build(ev_year, est0, est1, day0, exp, label):39 """Returns per-benchmark dicts of {ticker: AR series} for estimation and event windows (event = day0 .. +395 trading days, plus pre-block)."""40 trading = mkt_all.loc[(mkt_all.index >= est0)].index; pre_start_idx = trading.get_indexer([day0])[0]41 if pre_start_idx < 0: pre_start_idx = trading.searchsorted(day0)42 ev_dates = trading[pre_start_idx - 20: pre_start_idx + 395]; event_pos = {d: i - 20 for i, d in enumerate(ev_dates)} # position -20..+394, day 0 = first day >= day043 sub = dp[dp.ticker.isin(exp.ticker) & (dp.date >= est0) & (dp.date <= ev_dates[-1])].merge(mkt_all.rename("mkt"), left_on="date", right_index=True, how="left").dropna(subset=["ret", "mkt"])44 out = {b: {"est": {}, "ev": {}} for b in ["mm", "beta", "madj", "sizemom"]}; beta = {}; alpha = {}; sd_est = {}45 for t, g in sub.groupby("ticker"):46 est = g[(g.date >= est0) & (g.date <= est1)]47 if len(est) < 60: continue48 m = sm.OLS(est["ret"], sm.add_constant(est["mkt"])).fit(); a, b = m.params["const"], m.params["mkt"]; beta[t] = b; alpha[t] = a49 ev = g[g.date.isin(ev_dates)].set_index("date")50 for key, f_est, f_ev in [("mm", lambda e: e["ret"] - (a + b * e["mkt"]), lambda e: e["ret"] - (a + b * e["mkt"])), ("beta", lambda e: e["ret"] - b * e["mkt"], lambda e: e["ret"] - b * e["mkt"]), ("madj", lambda e: e["ret"] - e["mkt"], lambda e: e["ret"] - e["mkt"])]:51 out[key]["est"][t] = f_est(est.set_index("date")); out[key]["ev"][t] = f_ev(ev)52 sd_est[t] = float(out["mm"]["est"][t].std(ddof=2))53 # size x momentum matched benchmark: raw return minus the mean raw return of the firm's (size5, mom5) cell (excluding itself)54 cells = exp.set_index("ticker")[["size5", "mom5"]]; rets = pd.DataFrame({t: g.set_index("date")["ret"] for t, g in sub.groupby("ticker") if t in beta}); cell_of = {t: (cells.loc[t, "size5"], cells.loc[t, "mom5"]) for t in rets.columns if t in cells.index and pd.notna(cells.loc[t, "size5"]) and pd.notna(cells.loc[t, "mom5"])}55 groups = {}56 for t, c in cell_of.items(): groups.setdefault(c, []).append(t)57 for c, tks in groups.items():58 if len(tks) < 3: continue59 block = rets[tks]; s = block.sum(1); n = block.notna().sum(1)60 for t in tks:61 other = (s - block[t].fillna(0)) / (n - block[t].notna().astype(int)).replace(0, np.nan); ar = block[t] - other; out["sizemom"]["est"][t] = ar.loc[est0:est1].dropna(); out["sizemom"]["ev"][t] = ar.loc[ev_dates[0]:].dropna()62 return out, beta, alpha, sd_est, ev_dates, event_pos6364def hl_stats(ARs, hi, lo, event_pos, sd_est=None):65 """H-L daily portfolio statistics from AR matrices (dict ticker -> series). CDA t uses the estimation-window SD of the H-L portfolio."""66 AE = pd.DataFrame({t: ARs["est"][t] for t in hi + lo if t in ARs["est"]}); AV = pd.DataFrame({t: ARs["ev"][t] for t in hi + lo if t in ARs["ev"]})67 hi_ = [t for t in hi if t in AE.columns and t in AV.columns]; lo_ = [t for t in lo if t in AE.columns and t in AV.columns]68 if len(hi_) < 10 or len(lo_) < 10: return None, None69 hl_est = AE[hi_].mean(1) - AE[lo_].mean(1); hl_ev = (AV[hi_].mean(1) - AV[lo_].mean(1)); hl_ev.index = [event_pos.get(d, None) for d in hl_ev.index]; hl_ev = hl_ev[[i is not None for i in hl_ev.index]]70 s_est = hl_est.std(ddof=1); c = AE[hi_ + lo_].corr().values; n = c.shape[0]; rho = float((c.sum() - n) / (n * (n - 1)))71 out = {"n_high": len(hi_), "n_low": len(lo_), "rho_bar": round(rho, 4), "kp_factor": round(float(np.sqrt((1 - rho) / (1 + (n - 1) * rho))), 3), "horizons": {}}72 for lab, a, b in [("pre_m20_m1", -20, -1)] + HORIZONS:73 seg = hl_ev[(hl_ev.index >= a) & (hl_ev.index <= b)]; L = len(seg)74 if L == 0: continue75 car = seg.sum(); row = {"days": int(L), "HL_CAR_pct": round(car * 100, 2), "t_cda": round(car / (s_est * np.sqrt(L)), 2), "t_event_sd": round(car / (hl_ev[hl_ev.index >= 0].std(ddof=1) * np.sqrt(L)), 2) if (hl_ev.index >= 0).sum() > 5 else None, "High_pct": round(AV[hi_].loc[[d for d in AV.index if a <= event_pos.get(d, -99) <= b]].mean(1).sum() * 100, 2), "Low_pct": round(AV[lo_].loc[[d for d in AV.index if a <= event_pos.get(d, -99) <= b]].mean(1).sum() * 100, 2)}76 if sd_est is not None:77 dd = [d for d in AV.index if a <= event_pos.get(d, -99) <= b]; sc_h = np.array([AV[t].loc[dd].sum() / (sd_est[t] * np.sqrt(L)) for t in hi_ if t in sd_est and sd_est[t] > 0]); sc_l = np.array([AV[t].loc[dd].sum() / (sd_est[t] * np.sqrt(L)) for t in lo_ if t in sd_est and sd_est[t] > 0])78 row["t_bmp_HL"] = round(float((sc_h.mean() - sc_l.mean()) / np.sqrt(sc_h.var(ddof=1) / len(sc_h) + sc_l.var(ddof=1) / len(sc_l))), 2); row["t_bmp_HL_kp"] = round(row["t_bmp_HL"] * out["kp_factor"], 2)79 out["horizons"][lab] = row80 return out, hl_ev8182# ---------------- 2022 launch ----------------83exp = exposure(2021); est0, est1, day0 = pd.Timestamp("2022-01-03"), pd.Timestamp("2022-10-31"), pd.Timestamp("2022-12-01")84AR, beta, alpha, sd_est, ev_dates, event_pos = build(2022, est0, est1, day0, exp, "2022")85AR21, beta21, _, sd21, _, _ = build(2022, pd.Timestamp("2021-01-04"), pd.Timestamp("2021-12-31"), day0, exp, "2022_est2021")86RES["n_firms"] = len(beta); RES["window"] = {"estimation": "2022-01-03..2022-10-31 (209 days)", "day0": str(ev_dates[20].date()), "pre_block": f"{ev_dates[0].date()}..{ev_dates[19].date()}", "last_day": str(ev_dates[-1].date()), "n_event_days": int(len(ev_dates) - 20)}87SUBS = {"full": exp.ticker, "large": exp[exp.size_group == "large"].ticker, "small": exp[exp.size_group == "small"].ticker, "large_v1": exp[exp.size_group_v1 == "large"].ticker, "small_v1": exp[exp.size_group_v1 == "small"].ticker, "large_pct42": exp[exp.size_group_pct42 == "large"].ticker, "small_pct42": exp[exp.size_group_pct42 == "small"].ticker, "sp500": exp[exp.in_sp500 == True].ticker, "non_sp500": exp[exp.in_sp500 == False].ticker}88terc = exp.set_index("ticker")["tercile"]; terc_ws = exp.set_index("ticker")["tercile_within_size"]89RES["by_benchmark"] = {}; HLSERIES = {}90for bname, ARs, sde in [("mm_2022", AR["mm"], sd_est), ("beta_only", AR["beta"], None), ("market_adjusted", AR["madj"], None), ("mm_est2021", AR21["mm"], sd21), ("size_mom_matched", AR["sizemom"], None)]:91 RES["by_benchmark"][bname] = {}92 for sname, tks in SUBS.items():93 tks = set(tks); hi = [t for t in terc.index if terc[t] == "High" and t in tks]; lo = [t for t in terc.index if terc[t] == "Low" and t in tks]94 r, hl = hl_stats(ARs, hi, lo, event_pos, sde); RES["by_benchmark"][bname][sname] = r95 if r is not None and bname == "mm_2022": HLSERIES[sname] = hl96 for sname in ["large", "small"]: # terciles formed within size group97 tks = set(SUBS[sname]); hi = [t for t in terc_ws.index if terc_ws[t] == "High" and t in tks]; lo = [t for t in terc_ws.index if terc_ws[t] == "Low" and t in tks]98 r, _ = hl_stats(ARs, hi, lo, event_pos, sde); RES["by_benchmark"][bname][sname + "_within"] = r99 print(bname, {s: (v["horizons"]["car_1y"]["HL_CAR_pct"], v["horizons"]["car_1y"]["t_cda"]) for s, v in RES["by_benchmark"][bname].items() if v and "car_1y" in v["horizons"]}, flush=True)100# (e) FF5+MOM calendar-time alpha of the H-L daily portfolio (raw returns) over each horizon, HAC t101rets = pd.DataFrame({t: g.set_index("date")["ret"] for t, g in dp[dp.ticker.isin(beta) & dp.date.isin(ev_dates)].groupby("ticker")})102RES["ff6_calendar_time"] = {}103for sname, tks in SUBS.items():104 tks = set(tks); hi = [t for t in terc.index if terc[t] == "High" and t in tks and t in rets.columns]; lo = [t for t in terc.index if terc[t] == "Low" and t in tks and t in rets.columns]105 if len(hi) < 10 or len(lo) < 10: continue106 hl = (rets[hi].mean(1) - rets[lo].mean(1)).rename("hl").to_frame().join(ffd, how="inner"); hl["pos"] = [event_pos[d] for d in hl.index]; blk = {}107 for lab, a, b in HORIZONS:108 s = hl[(hl.pos >= a) & (hl.pos <= b)].dropna()109 if len(s) < 15: continue110 r = sm.OLS(s["hl"], sm.add_constant(s[["MktRF", "SMB", "HML", "RMW", "CMA", "MOM"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 5}); r1 = sm.OLS(s["hl"], sm.add_constant(s[["MktRF"]])).fit(cov_type="HAC", cov_kwds={"maxlags": 5})111 blk[lab] = {"days": int(len(s)), "ff6_alpha_cum_pct": round(float(r.params["const"] * len(s)) * 100, 2), "ff6_t": round(float(r.tvalues["const"]), 2), "capm_alpha_cum_pct": round(float(r1.params["const"] * len(s)) * 100, 2), "capm_t": round(float(r1.tvalues["const"]), 2), "betas": {k: round(float(r.params[k]), 3) for k in ["MktRF", "SMB", "HML", "RMW", "CMA", "MOM"]}}112 RES["ff6_calendar_time"][sname] = blk113print("FF6 calendar-time 1y:", {s: (v.get("car_1y", {}).get("ff6_alpha_cum_pct"), v.get("car_1y", {}).get("ff6_t")) for s, v in RES["ff6_calendar_time"].items()}, flush=True)114# firm-level CARs (baseline) + cross-section + Welch for the before/after115car_rows = []116for t in beta:117 s = AR["mm"]["ev"][t].copy(); s.index = [event_pos[d] for d in s.index]118 row = {"ticker": t, "beta": beta[t], "alpha_est": alpha[t]}119 for lab, a, b in [("pre_m20_m1", -20, -1)] + HORIZONS: row[lab] = float(s[(s.index >= a) & (s.index <= b)].sum())120 car_rows.append(row)121car = pd.DataFrame(car_rows).merge(exp, on="ticker"); car.to_parquet(f"{V}/event_car_v2.parquet", index=False)122RES["welch_reference"] = {}123for sname, tks in SUBS.items():124 c = car[car.ticker.isin(set(tks))]; hi, lo = c[c.tercile == "High"]["car_1y"].dropna(), c[c.tercile == "Low"]["car_1y"].dropna()125 if len(hi) > 10 and len(lo) > 10: tt = stats.ttest_ind(hi, lo, equal_var=False); RES["welch_reference"][sname] = {"n": int(len(c)), "H_minus_L_1y_pct": round((hi.mean() - lo.mean()) * 100, 2), "t_welch": round(float(tt.statistic), 2)}126reg = car.dropna(subset=["car_1y", "ai_pre", "log_mktcap", "beta"]).copy()127m1 = sm.OLS(reg["car_1y"], sm.add_constant(reg[["ai_pre", "log_mktcap", "beta"]])).fit(cov_type="HC1"); m2 = smf.ols("car_1y ~ ai_pre + log_mktcap + beta + C(sector)", data=reg).fit(cov_type="HC1")128RES["xsec_1y"] = {"ai_pre_pct": round(m1.params["ai_pre"] * 100, 3), "t": round(m1.tvalues["ai_pre"], 2), "n": int(m1.nobs), "r2": round(m1.rsquared, 3)}; RES["xsec_1y_sector"] = {"ai_pre_pct": round(m2.params["ai_pre"] * 100, 3), "t": round(m2.tvalues["ai_pre"], 2)}129RES["tercile_characteristics"] = car.groupby("tercile", observed=True).agg(beta=("beta", "mean"), log_mktcap=("log_mktcap", "mean"), share_large=("size_group", lambda s: (s == "large").mean()), n=("ticker", "size")).round(3).to_dict("index")130# CAR paths (baseline) for the figure131AV = pd.DataFrame({t: AR["mm"]["ev"][t] for t in beta}); paths = {tc: AV[[t for t in AV.columns if terc.get(t) == tc]].mean(1).cumsum() for tc in ["Low", "Mid", "High"]}; pdf = pd.DataFrame(paths); pdf["High_minus_Low"] = pdf["High"] - pdf["Low"]; pdf["pos"] = [event_pos[d] for d in pdf.index]; pdf.to_parquet(f"{V}/event_car_paths_v2.parquet")132# CDA table (baseline, full sample) in the author's layout + all sub-samples133cda_rows = []134for sname in SUBS:135 r = RES["by_benchmark"]["mm_2022"][sname]136 if r is None: continue137 for lab, h in r["horizons"].items(): cda_rows.append({"sample": sname, "horizon": lab, "days": h["days"], "HL_CAR_pct": h["HL_CAR_pct"], "t_cda": h["t_cda"], "t_event_sd": h["t_event_sd"], "t_bmp_HL": h.get("t_bmp_HL"), "kp_factor": r["kp_factor"], "rho_bar": r["rho_bar"], "n_high": r["n_high"], "n_low": r["n_low"]})138pd.DataFrame(cda_rows).to_csv(f"{T}/t_event_cda.csv", index=False)139# multiple testing over the event family: sub-samples x horizons (baseline), Holm / BH on normal p-values of the CDA t140fam = [(s, lab, h["t_cda"]) for s in ["full", "large", "small", "sp500", "non_sp500"] for lab, h in RES["by_benchmark"]["mm_2022"][s]["horizons"].items() if lab != "pre_m20_m1"]141p = 2 * (1 - stats.norm.cdf(np.abs([x[2] for x in fam]))); RES["mht_event"] = {"members": [f"{s}:{l}" for s, l, _ in fam], "t_cda": [x[2] for x in fam], "p": [round(float(x), 4) for x in p], "p_holm": [round(float(x), 4) for x in I.holm(p)], "p_bh": [round(float(x), 4) for x in I.bh(p)], "family_size": len(fam)}142143# ---------------- placebo launches: 30 Nov 2015 .. 2021 ----------------144RES["placebo_launches"] = {}145for Y in range(2015, 2022):146 try:147 e = exposure(Y - 1); e = e[e.ticker.isin(dp.ticker.unique())]148 ARp, bp, ap, sdp, evd, evp = build(Y, pd.Timestamp(f"{Y}-01-03"), pd.Timestamp(f"{Y}-10-31"), pd.Timestamp(f"{Y}-12-01"), e, str(Y)); tp = e.set_index("ticker")["tercile"]149 blk = {"n_firms": len(bp)}150 for sname, tks in [("full", e.ticker), ("large", e[e.size_group == "large"].ticker), ("small", e[e.size_group == "small"].ticker)]:151 tks = set(tks); hi = [t for t in tp.index if tp[t] == "High" and t in tks]; lo = [t for t in tp.index if tp[t] == "Low" and t in tks]; r, _ = hl_stats(ARp["mm"], hi, lo, evp, None)152 blk[sname] = {h: {"HL_CAR_pct": r["horizons"][h]["HL_CAR_pct"], "t_cda": r["horizons"][h]["t_cda"]} for h in ["car_0_20", "car_6m", "car_1y", "car_18m"] if r and h in r["horizons"]} if r else None153 RES["placebo_launches"][str(Y)] = blk; print("placebo", Y, blk["full"].get("car_1y") if blk["full"] else None, flush=True)154 except Exception as ex:155 RES["placebo_launches"][str(Y)] = {"error": repr(ex)[:200]}; print("placebo", Y, "failed", repr(ex)[:200], flush=True)156obs = RES["by_benchmark"]["mm_2022"]["full"]["horizons"]["car_1y"]; vals = {y: v["full"]["car_1y"] for y, v in RES["placebo_launches"].items() if v.get("full") and "car_1y" in v["full"]}157RES["placebo_summary"] = {"HL_1y_by_launch_pct": {y: v["HL_CAR_pct"] for y, v in vals.items()} | {"2022": obs["HL_CAR_pct"]}, "t_cda_by_launch": {y: v["t_cda"] for y, v in vals.items()} | {"2022": obs["t_cda"]},158 "rank_2022_by_HL": int(1 + sum(v["HL_CAR_pct"] > obs["HL_CAR_pct"] for v in vals.values())), "rank_2022_by_t": int(1 + sum(v["t_cda"] > obs["t_cda"] for v in vals.values())), "n_launches": len(vals) + 1,159 "large_HL_1y_by_launch_pct": {y: v["large"]["car_1y"]["HL_CAR_pct"] for y, v in RES["placebo_launches"].items() if v.get("large") and "car_1y" in v["large"]} | {"2022": RES["by_benchmark"]["mm_2022"]["large"]["horizons"]["car_1y"]["HL_CAR_pct"]}}160RES["runtime_s"] = round(time.time() - T0); json.dump(RES, open(f"{T}/results_event.json", "w"), indent=2, default=str); print("saved results_event.json in", RES["runtime_s"], "s; placebo", RES["placebo_summary"], flush=True)161