1#!/usr/bin/env python2"""Phase 5: negative controls, subgroup stress tests (H3), panel regressions (H3), pricing tests (H5)."""3from __future__ import annotations45import json6import sys7from pathlib import Path89import numpy as np10import pandas as pd11import polars as pl1213sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))14from wp39.inference import binary_metrics, paired_block_bootstrap, pointwise_logloss, quarter_blocks, two_way_cluster_ols # noqa: E40215from wp39.models import Preprocessor, columns_for, predict_proba # noqa: E40216from wp39.pipeline import Paths, get_logger, load_config # noqa: E40217from wp39.portfolio import car_window, download_french, factor_alpha, monthly_decile_portfolios, read_french_daily # noqa: E40218from wp39.provenance import now_iso, write_json # noqa: E40219from wp39.robustness import add_terciles, shuffle_labels_within_quarter, subgroup_masks # noqa: E40220from wp39.samples import feature_columns # noqa: E40221from wp39.splits import chronological # noqa: E4022223KEY = ["entity_id", "fiscal_period_end", "cutoff"]24FIXED = {"num_leaves": 15, "learning_rate": 0.03, "min_child_samples": 100, "feature_fraction": 0.7, "bagging_fraction": 0.8,25 "bagging_freq": 1, "lambda_l2": 1.0, "n_estimators": 600, "verbose": -1, "n_jobs": 12}262728def fixed_fit(train, test, cols, target, seed):29 import lightgbm as lgb30 pp = Preprocessor(cols).fit(train)31 cut = train["cutoff"].quantile(0.85)32 tr, va = train.filter(pl.col("cutoff") < cut), train.filter(pl.col("cutoff") >= cut)33 m = lgb.LGBMClassifier(random_state=seed, **FIXED)34 m.fit(pp.transform(tr), tr[target].to_numpy().astype(int), eval_set=[(pp.transform(va), va[target].to_numpy().astype(int))],35 callbacks=[lgb.early_stopping(50, verbose=False)])36 return predict_proba(m, pp.transform(test))373839def paired(A: pl.DataFrame, B: pl.DataFrame, y: str, nb: int, seed: int) -> dict:40 j = A.select(KEY + [y, "p", "A_q"]).join(B.select(KEY + [pl.col("p").alias("pb")]), on=KEY)41 yy = j[y].to_numpy().astype(int)42 d = pointwise_logloss(yy, j["p"].to_numpy()) - pointwise_logloss(yy, j["pb"].to_numpy())43 r = paired_block_bootstrap(d, quarter_blocks(j["A_q"].to_numpy()), nb, seed)44 r["auc_a"] = binary_metrics(yy, j["p"].to_numpy()).get("auc"); r["auc_b"] = binary_metrics(yy, j["pb"].to_numpy()).get("auc")45 r["base_rate"] = float(yy.mean())46 return r474849def main() -> int:50 P = Paths.default()51 log = get_logger("08_robust", P.logs / "08_robust.log")52 cfg = load_config("research"); hp, nb, seed = cfg["cutoffs"]["primary_horizon_td"], cfg["inference"]["n_bootstrap"], cfg["seed"]53 y = "y_adverse_opinc"54 df = pl.read_parquet(P.processed / "sample_event.parquet").filter((pl.col("h") == hp) & pl.col(y).is_not_null())55 fc = feature_columns(df)56 sp = chronological(df, cfg["splits"])57 trv = pl.concat([sp["train"], sp["validation"]])5859 # ---------------------------------------------------------------- placebo 1: labels shuffled within calendar quarter60 plc = []61 for s in range(3):62 shuf_tr = shuffle_labels_within_quarter(trv, y, seed + s); shuf_te = shuffle_labels_within_quarter(sp["test"], y, seed + 100 + s)63 A = shuf_te.with_columns(pl.Series("p", fixed_fit(shuf_tr, shuf_te, columns_for("B5", fc), y, seed)))64 B = shuf_te.with_columns(pl.Series("p", fixed_fit(shuf_tr, shuf_te, columns_for("B2", fc), y, seed)))65 r = paired(A, B, y, nb, seed); r.update({"placebo": "shuffled_labels", "draw": s}); plc.append(r)66 log.info("placebo shuffled labels draw %d: d=%.5f [%.5f, %.5f] auc5=%.3f", s, r["mean"], r["ci_lo"], r["ci_hi"], r["auc_a"])67 pd.DataFrame(plc).to_csv(P.evaluation / "placebo_shuffled.csv", index=False)6869 # ---------------------------------------------------------------- subgroups on the chronological LightGBM test predictions70 preds = {m: pl.read_parquet(P.evaluation / "preds" / f"event_h{hp}_{y}_{m}_test_chrono.parquet") for m in ("LB2", "LB5", "B2", "B5")}71 tr_t, te_t = add_terciles(sp["train"], [sp["test"]])72 te_t = te_t.with_columns([73 ((pl.col("cutoff") - pl.col("A_prev")).dt.total_days() >= 35).fill_null(True).alias("window_clean"),74 (pl.col("w1m_present_share") >= 0.98).alias("dense_bars"),75 ])76 extra = {"window_clean": pl.col("window_clean"), "window_straddles_prev_release": ~pl.col("window_clean"),77 "dense_bars": pl.col("dense_bars"), "sparse_bars": ~pl.col("dense_bars")}78 sub = []79 for fam, a, b in (("lightgbm", "LB5", "LB2"), ("logistic", "B5", "B2")):80 for name, mask in {**subgroup_masks(te_t), **extra}.items():81 ids = te_t.filter(mask).select(KEY)82 A, B = preds[a].join(ids, on=KEY), preds[b].join(ids, on=KEY)83 if A.height < 500:84 continue85 r = paired(A, B, y, nb, seed); r.update({"family": fam, "subgroup": name}); sub.append(r)86 subdf = pd.DataFrame(sub); subdf.to_csv(P.evaluation / "subgroups.csv", index=False)87 log.info("subgroups:\n%s", subdf[subdf.family == "lightgbm"][["subgroup", "n", "mean", "ci_lo", "ci_hi", "auc_a", "auc_b"]].to_string())8889 # ---------------------------------------------------------------- H3 panel regressions on walk-forward OOF signals90 oof5 = pl.read_parquet(P.evaluation / "preds" / f"event_h{hp}_{y}_LB5_oof_walkforward.parquet").select(KEY + [pl.col("p").alias("lfds5")])91 oof2 = pl.read_parquet(P.evaluation / "preds" / f"event_h{hp}_{y}_LB2_oof_walkforward.parquet").select(KEY + [pl.col("p").alias("lfds2")])92 pan = df.join(oof5, on=KEY).join(oof2, on=KEY)93 pan = pan.with_columns([94 (pl.col("cutoff").dt.year().cast(pl.Utf8) + "Q" + ((pl.col("cutoff").dt.month() - 1) // 3 + 1).cast(pl.Utf8)).alias("cq"),95 ])96 pan = pan.with_columns((pl.col("sic_division") + "|" + pl.col("cq")).alias("sector_q"))97 # ex-ante cross-sectional terciles within calendar quarter98 for name, col in (("size", "f_log_mcap"), ("illiq", "wd_amihud60"), ("lev", "f_leverage")):99 pan = pan.with_columns(((pl.col(col).rank("average").over("cq") - 1) / pl.len().over("cq") * 3).floor().clip(0, 2).cast(pl.Int8).alias(f"t_{name}"))100 pd_ = pan.to_pandas()101 lo, hi = np.quantile(pd_["S_opinc"].dropna(), [0.01, 0.99]); pd_["S_w"] = pd_["S_opinc"].clip(lo, hi)102 pd_["d_lfds"] = pd_["lfds5"] - pd_["lfds2"]103 for name in ("size", "illiq", "lev"):104 for k in (0, 2):105 pd_[f"lfds5_x_{name}{k}"] = pd_["lfds5"] * (pd_[f"t_{name}"] == k)106 pd_[f"dlfds_x_{name}{k}"] = pd_["d_lfds"] * (pd_[f"t_{name}"] == k)107 regs = {}108 controls = ["f_roa", "f_leverage", "f_log_mcap", "wd_ret20", "wd_vol20", "last_surprise", "mean_surprise4"]109 regs["S_on_lfds5"] = two_way_cluster_ols(pd_, "S_w", ["lfds5"] + controls, "entity_id", "sector_q", "entity_id", "cq")110 regs["S_on_lfds2_and_dlfds"] = two_way_cluster_ols(pd_, "S_w", ["lfds2", "d_lfds"] + controls, "entity_id", "sector_q", "entity_id", "cq")111 regs["y_on_lfds5"] = two_way_cluster_ols(pd_.assign(yb=pd_[y].astype(float)), "yb", ["lfds5"] + controls, "entity_id", "sector_q", "entity_id", "cq")112 for name in ("size", "illiq", "lev"):113 regs[f"dlfds_x_{name}"] = two_way_cluster_ols(pd_, "S_w", ["lfds2", "d_lfds", f"dlfds_x_{name}0", f"dlfds_x_{name}2"] + controls,114 "entity_id", "sector_q", "entity_id", "cq")115 allreg = pd.concat([r.assign(spec=k) for k, r in regs.items()]); allreg.to_csv(P.evaluation / "panel_regressions.csv", index=False)116 log.info("panel:\n%s", allreg[allreg.spec.isin(["S_on_lfds5", "S_on_lfds2_and_dlfds"])][["spec", "var", "coef", "se", "t", "p", "n"]].head(12).to_string())117118 # ---------------------------------------------------------------- H5a: announcement-window CAR on the ex-ante signal119 ext = P.cache / "external"120 ff = read_french_daily(download_french(load_config("data_contract")["external"]["ff5_daily"], ext / "ff5_daily.zip"))121 mom = read_french_daily(download_french(load_config("data_contract")["external"]["mom_daily"], ext / "mom_daily.zip"))122 ff = ff.join(mom.iloc[:, :1].rename(columns={mom.columns[0]: "Mom"}), how="left")123 daily = pl.concat([pl.read_parquet(p).select(["symbol", "session_date", "d_ret_cc", "d_close", "d_dollar_vol"]) for p in (P.cache / "features").glob("*.parquet")])124 syms = set(pan["symbol"].unique().to_list())125 daily = daily.filter(pl.col("symbol").is_in(list(syms)))126 car = car_window(daily, ff, pan.select(["entity_id", "symbol", "cutoff", "A_q", "fiscal_period_end"]))127 pc = pd_.merge(car.to_pandas(), on=["entity_id", "symbol", "cutoff", "A_q", "fiscal_period_end"], how="inner")128 lo, hi = np.quantile(pc["car"].dropna(), [0.01, 0.99]); pc["car_w"] = pc["car"].clip(lo, hi)129 h5 = {"car_on_lfds5": two_way_cluster_ols(pc, "car_w", ["lfds5"] + controls, "entity_id", "sector_q", "entity_id", "cq"),130 "car_on_lfds2_dlfds": two_way_cluster_ols(pc, "car_w", ["lfds2", "d_lfds"] + controls, "entity_id", "sector_q", "entity_id", "cq")}131 pd.concat([r.assign(spec=k) for k, r in h5.items()]).to_csv(P.evaluation / "h5_car_regressions.csv", index=False)132 log.info("H5 CAR:\n%s", h5["car_on_lfds2_dlfds"][["var", "coef", "se", "t", "p", "n"]].head(3).to_string())133134 # ---------------------------------------------------------------- H5b: monthly decile sorts (calendar design, OOF scores)135 port_out = {}136 for m in ("LB5", "LB2"):137 f = P.evaluation / "preds" / f"calendar_h0_{y}_{m}_oof_walkforward.parquet"138 if not f.exists():139 continue140 sig = pl.read_parquet(f).select(["symbol", "cutoff", pl.col("p").alias("score")]).to_pandas()141 res = monthly_decile_portfolios(sig, daily, ff, cost_bps=20.0)142 port = res["portfolios"]; port.to_csv(P.evaluation / f"portfolios_{m}.csv", index=False)143 port_out[m] = {"gross": factor_alpha(port, ff, 1, 10, net=False), "net20bps": factor_alpha(port, ff, 1, 10, net=True),144 "decile_means": port.groupby("decile")["ret"].mean().round(5).to_dict(),145 "avg_turnover": float(port["turnover"].mean())}146 log.info("portfolio %s: %s", m, port_out[m]["gross"])147 write_json(P.evaluation / "h5_portfolios.json", port_out)148149 for f in ("placebo_shuffled", "subgroups", "panel_regressions", "h5_car_regressions"):150 (P.results / f"{f}.csv").write_text((P.evaluation / f"{f}.csv").read_text())151 write_json(P.results / "robustness_summary.json", {"run_at": now_iso(), "placebo": plc, "portfolios": port_out})152 return 0153154155if __name__ == "__main__":156 raise SystemExit(main())157