1"""r3k_05_event_study.py — ChatGPT (2022-11-30) market-model event study on the Russell 3000 panel, with2dependence-robust inference (Brown-Warner CDA; Kolari-Pynnonen deflation) and size sub-samples."""3import json, warnings, numpy as np, pandas as pd, statsmodels.api as sm, statsmodels.formula.api as smf4from scipy import stats5warnings.filterwarnings("ignore")6D = "data/processed_r3k"; T = "results/r3k"7panel = pd.read_parquet(f"{D}/analysis_panel.parquet")8dp = pd.read_parquet(f"{D}/prices_daily_event.parquet"); dp["date"] = pd.to_datetime(dp["date"]); dp = dp.sort_values(["ticker", "date"])9dp["ret"] = dp.groupby("ticker")["adj_close"].pct_change()10bm = pd.read_parquet(f"{D}/benchmarks_daily_event.parquet"); bm["date"] = pd.to_datetime(bm["date"]); bm = bm.sort_values(["ticker", "date"]); bm["ret"] = bm.groupby("ticker")["adj_close"].pct_change()11mkt = bm[bm.ticker == "IWV"][["date", "ret"]].rename(columns={"ret": "mkt"}) # Russell 3000 ETF as the market12t0 = pd.Timestamp("2022-11-30"); est0, est1 = pd.Timestamp("2022-01-03"), pd.Timestamp("2022-10-31"); ev0, ev1 = pd.Timestamp("2022-11-01"), pd.Timestamp("2024-06-30")13exp = panel[panel.fiscal_year == 2021].copy()14exp["ai_pre"] = (exp.ai_adoption_wmean.transform(lambda g: (g - g.mean()) / g.std()) + exp.ai_infrastructure_wmean.transform(lambda g: (g - g.mean()) / g.std())) / 215exp = exp[["ticker", "ai_pre", "log_mktcap", "sector", "size_group", "in_sp500"]].dropna(subset=["ai_pre"]).drop_duplicates("ticker"); exp["tercile"] = pd.qcut(exp.ai_pre, 3, labels=["Low", "Mid", "High"])16car_rows, daily_car, beta, ar_est, ar_ev = [], {}, {}, {}, {}17for t, g in dp[dp.ticker.isin(exp.ticker)].groupby("ticker"):18 g = g.merge(mkt, on="date", how="left").dropna(subset=["ret", "mkt"]); est = g[(g.date >= est0) & (g.date <= est1)]19 if len(est) < 60: continue20 m = sm.OLS(est["ret"], sm.add_constant(est["mkt"])).fit(); a, b = m.params["const"], m.params["mkt"]; beta[t] = b21 ar_est[t] = (est["ret"] - (a + b * est["mkt"])).set_axis(est["date"]); ev = g[(g.date >= ev0) & (g.date <= ev1)].copy(); ev["ar"] = ev["ret"] - (a + b * ev["mkt"]); ar_ev[t] = ev.set_index("date")["ar"]22 daily_car[t] = ev.set_index("date")["ar"].cumsum()23 for hz, lab in [(60, "car_60d"), (125, "car_6m"), (250, "car_1y"), (395, "car_18m")]: car_rows.append({"ticker": t, "hz": lab, "car": ev["ar"].iloc[:hz].sum()})24car = pd.DataFrame(car_rows).pivot(index="ticker", columns="hz", values="car").reset_index().merge(exp, on="ticker"); car["beta"] = car.ticker.map(beta); car.to_parquet(f"{D}/event_car.parquet", index=False)25allc = pd.DataFrame(daily_car).dropna(how="all"); paths = {tc: allc[[t for t in exp[exp.tercile == tc].ticker if t in allc.columns]].mean(axis=1) for tc in ["Low", "Mid", "High"]}26pdf = pd.DataFrame(paths); pdf["High_minus_Low"] = pdf["High"] - pdf["Low"]; pdf.to_parquet(f"{D}/event_car_paths.parquet")27RES = {"n_firms": int(len(car))}28for hz in ["car_60d", "car_6m", "car_1y", "car_18m"]:29 hi, lo = car[car.tercile == "High"][hz].dropna(), car[car.tercile == "Low"][hz].dropna(); tt = stats.ttest_ind(hi, lo, equal_var=False)30 RES[hz] = {"High": round(hi.mean() * 100, 2), "Mid": round(car[car.tercile == "Mid"][hz].mean() * 100, 2), "Low": round(lo.mean() * 100, 2), "H_minus_L": round((hi.mean() - lo.mean()) * 100, 2), "t_welch": round(tt.statistic, 2), "p": round(tt.pvalue, 4)}31reg = car.dropna(subset=["car_1y", "ai_pre", "log_mktcap", "beta"]).copy()32m1 = 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")33RES["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)}34RES["xsec_1y_sector"] = {"ai_pre_pct": round(m2.params["ai_pre"] * 100, 3), "t": round(m2.tvalues["ai_pre"], 2)}35# size sub-samples36RES["by_size"] = {}37for gname, msk in [("large", car.size_group == "large"), ("small", car.size_group == "small"), ("sp500", car.in_sp500 == True), ("non_sp500", car.in_sp500 == False)]:38 c = car[msk]; hi, lo = c[c.tercile == "High"]["car_1y"].dropna(), c[c.tercile == "Low"]["car_1y"].dropna()39 if len(hi) > 10 and len(lo) > 10:40 tt = stats.ttest_ind(hi, lo, equal_var=False); RES["by_size"][gname] = {"n": int(len(c)), "H_minus_L_1y": round((hi.mean() - lo.mean()) * 100, 2), "t_welch": round(tt.statistic, 2)}41# dependence-robust inference (CDA + KP)42AE, AV = pd.DataFrame(ar_est), pd.DataFrame(ar_ev); terc = car.drop_duplicates("ticker").set_index("ticker")["tercile"]43hi_t = [t for t in AE.columns if terc.get(t) == "High"]; lo_t = [t for t in AE.columns if terc.get(t) == "Low"]44hl_est = AE[hi_t].mean(axis=1) - AE[lo_t].mean(axis=1); hl_ev = (AV[hi_t].mean(axis=1) - AV[lo_t].mean(axis=1)).dropna(); s_est = hl_est.std(ddof=1)45def avg_corr(M): c = M.corr().values; n = c.shape[0]; return float((c.sum() - n) / (n * (n - 1)))46rho = avg_corr(AE[hi_t + lo_t]); n = len(hi_t) + len(lo_t); adj = np.sqrt((1 - rho) / (1 + (n - 1) * rho))47cda = []48for hz, lab in [(60, "car_60d"), (125, "car_6m"), (250, "car_1y"), (395, "car_18m")]:49 sub = hl_ev.iloc[:hz]; c = sub.sum(); cda.append({"horizon": lab, "days": len(sub), "HL_CAR_pct": round(c * 100, 2), "t_cda": round(c / (s_est * np.sqrt(len(sub))), 2), "t_event_sd": round(c / (hl_ev.std(ddof=1) * np.sqrt(len(sub))), 2), "kp_factor": round(adj, 3), "rho_bar": round(rho, 4)})50RES["cda"] = cda; RES["n_high"] = len(hi_t); RES["n_low"] = len(lo_t)51pd.DataFrame(cda).to_csv(f"{T}/t_event_cda.csv", index=False)52json.dump(RES, open(f"{T}/results_event.json", "w"), indent=2, default=str); print(json.dumps({k: RES[k] for k in ["car_60d", "car_1y", "xsec_1y", "by_size"]}, indent=0)); print(cda)53