1#!/usr/bin/env python2"""Phase 3/4 evaluation: every metric and every paired contrast is computed from saved prediction files.3Outputs (artifacts/evaluation + results/): metrics_all.csv, contrasts.csv, horizon_ladder.csv, secondary_targets.csv,4firm_holdout.csv, walkforward_by_year.csv, reliability.csv, continuous.csv, power.json."""5from __future__ import annotations67import json8import sys9from pathlib import Path1011import numpy as np12import pandas as pd13import polars as pl1415sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "src"))16from wp39.inference import (binary_metrics, continuous_metrics, dm_gw_test, holm, paired_block_bootstrap, # noqa: E40217 pointwise_logloss, quarter_blocks, reliability_table)18from wp39.pipeline import Paths, get_logger, load_config # noqa: E40219from wp39.provenance import now_iso, write_json # noqa: E4022021CONTRASTS = [("B5", "B2"), ("B4", "B2"), ("B3", "B2"), ("B5", "B4"), ("B5", "B3"), ("B5all", "B2"), ("B5all", "B5"), ("B2", "B1"), ("B1", "B0")]22FAMILIES = {"freq_vs_daily": [("B5", "B2"), ("B4", "B2"), ("B3", "B2")], "fine_vs_coarse": [("B5", "B4"), ("B5", "B3")]}232425def load_preds(P: Paths) -> pl.DataFrame:26 files = sorted((P.evaluation / "preds").glob("*.parquet"))27 frames = []28 for f in files:29 d = pl.read_parquet(f)30 target = [c for c in d.columns if c.startswith("y_") or c == "S_opinc"][0]31 d = d.rename({target: "y"}).with_columns(pl.lit(target).alias("target"))32 if "fold" not in d.columns:33 d = d.with_columns(pl.lit(None, dtype=pl.Int64).alias("fold"))34 d = d.with_columns([pl.col("fold").cast(pl.Int64), pl.col("h").cast(pl.Int64), pl.col("y").cast(pl.Float64), pl.col("p").cast(pl.Float64)])35 frames.append(d.select(["entity_id", "symbol", "fiscal_period_end", "h", "cutoff", "A_q", "sic_division", "release_time",36 "y", "p", "model", "part", "tag", "target", "fold"]))37 return pl.concat(frames)383940def contrast(a: pl.DataFrame, b: pl.DataFrame, n_boot: int, seed: int) -> dict:41 key = ["entity_id", "fiscal_period_end", "cutoff"]42 j = a.select(key + ["y", "p", "A_q"]).join(b.select(key + [pl.col("p").alias("pb")]), on=key, how="inner")43 y = j["y"].to_numpy().astype(int); pa, pb = j["p"].to_numpy(), j["pb"].to_numpy()44 d_ll = pointwise_logloss(y, pa) - pointwise_logloss(y, pb)45 d_br = (pa - y) ** 2 - (pb - y) ** 246 blocks = quarter_blocks(j["A_q"].to_numpy())47 bll, bbr = paired_block_bootstrap(d_ll, blocks, n_boot, seed), paired_block_bootstrap(d_br, blocks, n_boot, seed)48 gw = dm_gw_test(d_ll, blocks)49 return {"n": int(len(y)), "d_logloss": bll["mean"], "ci_lo": bll["ci_lo"], "ci_hi": bll["ci_hi"], "p_boot": bll["p_boot"],50 "gw_t": gw["t"], "gw_p": gw["p"], "n_quarters": gw["T"], "d_brier": bbr["mean"], "brier_ci_lo": bbr["ci_lo"], "brier_ci_hi": bbr["ci_hi"]}515253def main() -> int:54 P = Paths.default()55 log = get_logger("07_eval", P.logs / "07_eval.log")56 cfg = load_config("research")57 hp, nb, seed = cfg["cutoffs"]["primary_horizon_td"], cfg["inference"]["n_bootstrap"], cfg["seed"]58 preds = load_preds(P)59 log.info("prediction rows: %d (%d files)", preds.height, len(list((P.evaluation / "preds").glob("*.parquet"))))6061 # 1) metrics for every binary prediction set62 rows = []63 grp = preds.filter(pl.col("target") != "S_opinc").with_columns(pl.col("cutoff").dt.year().alias("year"))64 for keys, g in grp.group_by(["h", "target", "model", "part", "tag"], maintain_order=True):65 h, target, model, part, tag = keys66 m = binary_metrics(g["y"].to_numpy().astype(int), g["p"].to_numpy())67 rows.append({"h": h, "target": target, "model": model, "part": part, "tag": tag, **m})68 metrics = pd.DataFrame(rows).sort_values(["target", "h", "tag", "part", "model"])69 metrics.to_csv(P.evaluation / "metrics_all.csv", index=False)7071 # 2) primary contrasts (h = hp, primary target, chronological test) for logistic and LightGBM72 out = []73 base = preds.filter((pl.col("h") == hp) & (pl.col("target") == "y_adverse_opinc") & (pl.col("part") == "test") & (pl.col("tag") == "chrono"))74 for prefix in ("", "L"):75 pv = {}76 for a, b in CONTRASTS:77 A, B = base.filter(pl.col("model") == prefix + a), base.filter(pl.col("model") == (prefix + b if b != "B0" else "B0"))78 if A.height == 0 or B.height == 0:79 continue80 c = contrast(A, B, nb, seed); c.update({"family": "logistic" if prefix == "" else "lightgbm", "a": prefix + a, "b": prefix + b})81 out.append(c); pv[(a, b)] = c["gw_p"]82 for fam, pairs in FAMILIES.items():83 adj = holm({f"{a}-{b}": pv[(a, b)] for a, b in pairs if (a, b) in pv})84 for r in out:85 if r["family"] == ("logistic" if prefix == "" else "lightgbm") and f"{r['a'][len(prefix):]}-{r['b'][len(prefix):]}" in adj:86 r["holm_family"], r["gw_p_holm"] = fam, adj[f"{r['a'][len(prefix):]}-{r['b'][len(prefix):]}"]87 contrasts = pd.DataFrame(out); contrasts.to_csv(P.evaluation / "contrasts.csv", index=False)88 log.info("primary contrasts:\n%s", contrasts[["family", "a", "b", "n", "d_logloss", "ci_lo", "ci_hi", "gw_p"]].to_string())8990 # 3) horizon ladder (B5 vs B2) logistic chrono + LightGBM fixed91 lad = []92 for h in sorted(preds["h"].unique().to_list()):93 for fam, tag, a, b in (("logistic", "chrono", "B5", "B2"), ("lightgbm", "fixed", "LB5", "LB2")):94 s = preds.filter((pl.col("h") == h) & (pl.col("target") == "y_adverse_opinc") & (pl.col("part") == "test") & (pl.col("tag") == tag))95 A, B = s.filter(pl.col("model") == a), s.filter(pl.col("model") == b)96 if A.height and B.height:97 c = contrast(A, B, nb, seed); c.update({"h": h, "family": fam})98 c["auc_a"] = binary_metrics(A["y"].to_numpy().astype(int), A["p"].to_numpy()).get("auc")99 c["auc_b"] = binary_metrics(B["y"].to_numpy().astype(int), B["p"].to_numpy()).get("auc")100 lad.append(c)101 ladder = pd.DataFrame(lad); ladder.to_csv(P.evaluation / "horizon_ladder.csv", index=False)102 h4 = {}103 for fam, g in ladder.groupby("family"):104 g = g.sort_values("h", ascending=False)105 from scipy.stats import spearmanr106 rho, p = spearmanr(-g["h"], -g["d_logloss"]) # gain (= -d) should rise as h falls107 h4[fam] = {"spearman_gain_vs_minus_h": float(rho), "p": float(p), "horizons": g["h"].tolist(), "d_logloss": g["d_logloss"].round(5).tolist()}108 write_json(P.evaluation / "h4_timing.json", h4)109110 # 4) secondary targets (B5 vs B2) logistic + lightgbm at hp111 sec = []112 for target in ("y_adverse_netinc", "y_roa_decline", "y_ocf_negative", "y_liquidity_deterioration", "y_abs_surprise", "y_adverse_opinc"):113 for fam, a, b in (("logistic", "B5", "B2"), ("lightgbm", "LB5", "LB2")):114 s = preds.filter((pl.col("h") == hp) & (pl.col("target") == target) & (pl.col("part") == "test") & (pl.col("tag") == "chrono"))115 A, B = s.filter(pl.col("model") == a), s.filter(pl.col("model") == b)116 if A.height and B.height:117 c = contrast(A, B, nb, seed); c.update({"target": target, "family": fam, "base_rate": float(A["y"].mean()),118 "auc_a": binary_metrics(A["y"].to_numpy().astype(int), A["p"].to_numpy()).get("auc"),119 "auc_b": binary_metrics(B["y"].to_numpy().astype(int), B["p"].to_numpy()).get("auc")})120 sec.append(c)121 pd.DataFrame(sec).to_csv(P.evaluation / "secondary_targets.csv", index=False)122123 # 5) firm hold-out124 fh = []125 for tag in ("firmhold_heldout", "firmhold_seen"):126 for fam, a, b in (("logistic", "B5", "B2"), ("lightgbm", "LB5", "LB2")):127 s = preds.filter((pl.col("h") == hp) & (pl.col("target") == "y_adverse_opinc") & (pl.col("part") == "test") & (pl.col("tag") == tag))128 A, B = s.filter(pl.col("model") == a), s.filter(pl.col("model") == b)129 if A.height and B.height:130 c = contrast(A, B, nb, seed); c.update({"tag": tag, "family": fam}); fh.append(c)131 pd.DataFrame(fh).to_csv(P.evaluation / "firm_holdout.csv", index=False)132133 # 6) walk-forward by year134 wf = []135 oof = preds.filter((pl.col("tag") == "walkforward") & (pl.col("h") == hp)).with_columns(pl.col("cutoff").dt.year().alias("year"))136 for keys, g in oof.group_by(["model", "year"], maintain_order=True):137 m = binary_metrics(g["y"].to_numpy().astype(int), g["p"].to_numpy()); wf.append({"model": keys[0], "year": keys[1], **m})138 wfd = pd.DataFrame(wf).sort_values(["model", "year"]); wfd.to_csv(P.evaluation / "walkforward_by_year.csv", index=False)139 wfc = []140 for fam, a, b in (("logistic", "B5", "B2"), ("lightgbm", "LB5", "LB2")):141 A, B = oof.filter(pl.col("model") == a), oof.filter(pl.col("model") == b)142 if A.height and B.height:143 for y in sorted(A["year"].unique().to_list()):144 c = contrast(A.filter(pl.col("year") == y), B.filter(pl.col("year") == y), 300, seed); c.update({"family": fam, "year": y}); wfc.append(c)145 c = contrast(A, B, nb, seed); c.update({"family": fam, "year": "all"}); wfc.append(c)146 pd.DataFrame(wfc).to_csv(P.evaluation / "walkforward_contrasts.csv", index=False)147148 # 7) reliability (test, chrono) for B2, B5, LB2, LB5149 rel = []150 for model in ("B2", "B5", "LB2", "LB5"):151 s = base.filter(pl.col("model") == model)152 if s.height:153 t = reliability_table(s["y"].to_numpy().astype(int), s["p"].to_numpy()); t["model"] = model; rel.append(t)154 if rel:155 pd.concat(rel).to_csv(P.evaluation / "reliability.csv", index=False)156157 # 8) continuous target158 cont = []159 fitlog = json.load(open(P.evaluation / "baselines_fit_log.json"))["fits"]160 for model in ("B2", "B5"):161 s = preds.filter((pl.col("target") == "S_opinc") & (pl.col("model") == model) & (pl.col("tag") == "chrono_cont"))162 if s.height:163 rec = next(r for r in fitlog if r.get("tag") == "chrono_cont" and r["model"] == model)164 lo, hi = rec["winsor"]165 y = np.clip(s["y"].to_numpy(), lo, hi)166 cont.append({"model": model, **continuous_metrics(y, s["p"].to_numpy(), rec["train_mean"])})167 if cont:168 A = preds.filter((pl.col("target") == "S_opinc") & (pl.col("model") == "B5") & (pl.col("tag") == "chrono_cont"))169 B = preds.filter((pl.col("target") == "S_opinc") & (pl.col("model") == "B2") & (pl.col("tag") == "chrono_cont"))170 key = ["entity_id", "fiscal_period_end", "cutoff"]171 j = A.join(B.select(key + [pl.col("p").alias("pb")]), on=key)172 lo, hi = next(r for r in fitlog if r.get("tag") == "chrono_cont" and r["model"] == "B5")["winsor"]173 y = np.clip(j["y"].to_numpy(), lo, hi)174 d = (j["p"].to_numpy() - y) ** 2 - (j["pb"].to_numpy() - y) ** 2175 bb = paired_block_bootstrap(d, quarter_blocks(j["A_q"].to_numpy()), nb, seed)176 cont.append({"model": "B5-B2 d_mse", **bb})177 pd.DataFrame(cont).to_csv(P.evaluation / "continuous.csv", index=False)178179 # 9) power (validation period, logistic B5 vs B2): MDE at 80 % power, alpha 5 %, quarter clustering180 v = preds.filter((pl.col("h") == hp) & (pl.col("target") == "y_adverse_opinc") & (pl.col("part") == "validation") & (pl.col("tag") == "chrono"))181 A, B = v.filter(pl.col("model") == "B5"), v.filter(pl.col("model") == "B2")182 power = {}183 if A.height and B.height:184 key = ["entity_id", "fiscal_period_end", "cutoff"]185 j = A.join(B.select(key + [pl.col("p").alias("pb")]), on=key)186 y = j["y"].to_numpy().astype(int)187 d = pointwise_logloss(y, j["p"].to_numpy()) - pointwise_logloss(y, j["pb"].to_numpy())188 bl = quarter_blocks(j["A_q"].to_numpy()); ub = np.unique(bl)189 dq = np.array([d[bl == b].mean() for b in ub]); se = dq.std(ddof=1) / np.sqrt(len(ub))190 n_test_q = preds.filter((pl.col("part") == "test") & (pl.col("tag") == "chrono") & (pl.col("h") == hp))["A_q"].to_numpy()191 T_test = len(np.unique(quarter_blocks(n_test_q))) if len(n_test_q) else 23192 se_test = dq.std(ddof=1) / np.sqrt(T_test)193 power = {"validation_quarters": int(len(ub)), "sd_quarterly_diff": float(dq.std(ddof=1)), "test_quarters_assumed": int(T_test),194 "mde_logloss_80pct_power": float((1.96 + 0.84) * se_test), "validation_mean_diff": float(d.mean())}195 write_json(P.evaluation / "power.json", power)196197 # copies to results/198 for f in ("metrics_all", "contrasts", "horizon_ladder", "secondary_targets", "firm_holdout", "walkforward_by_year",199 "walkforward_contrasts", "reliability", "continuous"):200 src = P.evaluation / f"{f}.csv"201 if src.exists():202 (P.results / f"{f}.csv").write_text(src.read_text())203 write_json(P.results / "evaluation_summary.json", {"run_at": now_iso(), "power": power, "h4": h4,204 "primary": contrasts.to_dict(orient="records") if len(contrasts) else []})205 print(contrasts[["family", "a", "b", "n", "d_logloss", "ci_lo", "ci_hi", "gw_p"]].to_string() if len(contrasts) else "no contrasts")206 return 0207208209if __name__ == "__main__":210 raise SystemExit(main())211