1"""H2 (and H4/H5) on Seshat: Y1 crisis severity (primary), Y1b severe binary, Y2 polity duration (Cox), Y3 succession.2This is the first module that reads outcomes together with P. Usage: python -m wp35.econometrics.h2_seshat3Everything follows research/preregistration.md §5.2–5.3, §5.5–5.7, §6.4"""5from __future__ import annotations6import json7from pathlib import Path8import numpy as np9import pandas as pd10import statsmodels.api as sm11import statsmodels.formula.api as smf12from lifelines import CoxPHFitter13from lifelines.statistics import proportional_hazard_test14from wp35.econometrics.common import fit_formula, coef_row, cv_compare, fold_win_share, hinge_fit, binary_metrics, oos_loglik_binomial1516ROOT = Path(__file__).resolve().parents[3]17PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical"18RES.mkdir(parents=True, exist_ok=True)19BIN = sm.families.Binomial()20FE = "C(macro_region) + C(era)"21XS = "X_warfare + X_prod_infra + epidemic_v"222324def load_crisis():25 cr = pd.read_parquet(PROC / "seshat_crisis_indices.parquet")26 y = pd.read_parquet(PROC / "seshat_crisis.parquet")[["crisis_id", "y1_severity", "y1b_severe", "y_reform", "epidemic_v", "y1_n_coded"]]27 d = cr.merge(y, on="crisis_id")28 meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "centralization", "homeNGA"]]29 d = d.merge(meta, on="polity_id", how="left", suffixes=("", "_m"))30 d["homeNGA"] = d["homeNGA"].fillna(d.get("homeNGA_m"))31 # persistence baseline: an earlier crisis record in the same NGA32 d = d.sort_values("crisis_start")33 d["prior_crisis_nga"] = 0.034 for i, r in d.iterrows():35 d.loc[i, "prior_crisis_nga"] = float(((d.homeNGA == r.homeNGA) & (d.crisis_start < r.crisis_start)).any())36 d["quasi"] = d.centralization.eq("quasi-polity")37 d["P_placebo"] = (d.D_pc1 - d.D_pc1.mean()) / d.D_pc1.std() - (d.X_prod_infra - d.X_prod_infra.mean()) / d.X_prod_infra.std()38 for c in ["P_main", "D_pc1", "K_mean", "K_irt", "complexity_pc1", "log_pop"]:39 d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std()40 return d414243def primary_models(d: pd.DataFrame):44 rows = []45 specs = {46 "A1 primary: Y1 ~ P + FE": (f"y1_severity ~ P_main_z + {FE}", "P_main_z"),47 "A2: Y1 ~ P + X + FE": (f"y1_severity ~ P_main_z + {XS} + {FE}", "P_main_z"),48 "A3: Y1 ~ P + D + K + X + FE (P given components)": (f"y1_severity ~ P_main_z + complexity_pc1_z + {XS} + {FE}", "P_main_z"),49 "B1 baseline D only": (f"y1_severity ~ D_pc1_z + {XS} + {FE}", "D_pc1_z"),50 "B2 baseline K only": (f"y1_severity ~ K_mean_z + {XS} + {FE}", "K_mean_z"),51 "B3 baseline complexity PC1": (f"y1_severity ~ complexity_pc1_z + {XS} + {FE}", "complexity_pc1_z"),52 "B4 baseline persistence": (f"y1_severity ~ prior_crisis_nga + {XS} + {FE}", "prior_crisis_nga"),53 "B5 D and K separately": (f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}", "D_pc1_z"),54 "Placebo: P built with production infrastructure instead of K": (f"y1_severity ~ P_placebo + {XS} + {FE}", "P_placebo"),55 }56 for label, (f, term) in specs.items():57 r = fit_formula(f, d, BIN, groups="seshat_region"); rows.append(coef_row(r, term, label))58 if "B5" in label:59 rows.append(coef_row(r, "K_mean_z", label))60 # binary secondary61 # binary secondary outcome: several region/era cells are all-severe (quasi-separation) -> region FE only, HC1 SE62 for label, f in {"A1b: Y1b ~ P + region FE (logit, HC1; era FE dropped: separation)": "y1b_severe ~ P_main_z + C(macro_region)",63 "A2b: Y1b ~ P + X + region FE (logit, HC1)": f"y1b_severe ~ P_main_z + {XS} + C(macro_region)",64 "A3b: Y1b ~ P (logit, HC1, no FE)": "y1b_severe ~ P_main_z"}.items():65 r = fit_formula(f, d, BIN, groups=None); rows.append(coef_row(r, "P_main_z", label))66 # exploratory: reforms67 r = fit_formula(f"y_reform ~ P_main_z + {XS} + {FE}", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "Exploratory: reform consequences ~ P + X + FE"))68 return pd.DataFrame(rows)697071def robustness(d: pd.DataFrame):72 rows = []73 variants = {"P_D_pc1_K_irt": "K = IRT latent", "P_log_pop_K_mean": "D = log population only", "P_D_pc1_K_pc1": "K = PC1 of blocks",74 "P_main_lo": "D at lower bounds", "P_main_hi": "D at upper bounds", "P_main_noinf": "inferred codes dropped",75 "P_within_era": "standardized within era", "P_within_region": "standardized within macro-region", "P_rank": "rank-based P"}76 for v, lab in variants.items():77 if v not in d.columns:78 # crisis-level file lacks some polity-level variants → merge from polity file79 continue80 dd = d.copy(); dd["_P"] = (dd[v] - dd[v].mean()) / dd[v].std()81 r = fit_formula(f"y1_severity ~ _P + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "_P", f"R: {lab}"))82 dd = d[~d.quasi]; r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: quasi-polities dropped"))83 dd = d[d.years_since_start >= 25]; r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: crises starting ≥25 years after polity start"))84 dd = d[d.y1_n_coded >= 12]; r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: ≥12 of 14 items coded"))85 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + C(macro_region)", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: no era FE"))86 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + C(era)", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: no region FE"))87 r = fit_formula(f"y1_severity ~ P_main_z + {XS}", d, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_main_z", "R: no FE"))88 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", d, BIN, groups="homeNGA"); rows.append(coef_row(r, "P_main_z", "R: SE clustered by NGA"))89 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {FE}", d, BIN, groups="polity_id"); rows.append(coef_row(r, "P_main_z", "R: SE clustered by polity"))90 return pd.DataFrame(rows)919293def out_of_sample(d: pd.DataFrame):94 formulas = {"A2 P": f"y1_severity ~ P_main_z + {XS} + {FE}", "B0 FE only": f"y1_severity ~ {XS} + {FE}", "B1 D only": f"y1_severity ~ D_pc1_z + {XS} + {FE}",95 "B2 K only": f"y1_severity ~ K_mean_z + {XS} + {FE}", "B3 complexity": f"y1_severity ~ complexity_pc1_z + {XS} + {FE}", "B4 persistence": f"y1_severity ~ prior_crisis_nga + {XS} + {FE}",96 "B5 D+K": f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}"}97 # era FE with leave-one-region-out can leave unseen levels; macro-region FE cannot be used when the region is held out98 f_loro = {k: v.replace("C(macro_region) + ", "") for k, v in formulas.items()}99 loro = cv_compare(d, f_loro, "y1_severity", BIN, fold_col="macro_region", fractional=True)100 kf = cv_compare(d, formulas, "y1_severity", BIN, n_repeats=10, k=5, fractional=True)101 loro["scheme"] = "leave-one-macro-region-out"; kf["scheme"] = "10x5-fold"102 cv = pd.concat([loro, kf]); cv.to_csv(RES / "h2_seshat_cv_folds.csv", index=False)103 wins = []104 for scheme, c in [("loro", loro), ("kfold", kf)]:105 for b in ["B0 FE only", "B1 D only", "B2 K only", "B3 complexity", "B4 persistence", "B5 D+K"]:106 w = fold_win_share(c, "A2 P", b); w["scheme"] = scheme; wins.append(w)107 wins = pd.DataFrame(wins); wins.to_csv(RES / "h2_seshat_cv_wins.csv", index=False)108 # binary Y1b out-of-sample metrics109 fb = {k: v.replace("y1_severity", "y1b_severe") for k, v in formulas.items()}110 kfb = cv_compare(d, fb, "y1b_severe", BIN, n_repeats=10, k=5)111 kfb.to_csv(RES / "h2_seshat_cv_folds_binary.csv", index=False)112 summ = kfb.groupby("model")[["oos_loglik", "auc", "brier"]].mean().reset_index(); summ.to_csv(RES / "h2_seshat_cv_binary_summary.csv", index=False)113 return wins, summ114115116def h4_nonlinear(d: pd.DataFrame):117 """Linear vs hinge vs spline, out of sample (10x5-fold); hinge knot chosen on the training fold only."""118 rng = np.random.default_rng(20261101)119 base = f"y1_severity ~ P_main_z + {XS} + {FE}"120 dd = d.dropna(subset=["y1_severity", "P_main_z", "X_warfare", "X_prod_infra", "epidemic_v", "macro_region", "era"]).reset_index(drop=True)121 rows = []122 for rep in range(10):123 perm = rng.permutation(len(dd)); parts = np.array_split(perm, 5)124 for i, te in enumerate(parts):125 m = np.zeros(len(dd), bool); m[te] = True; tr, ts = dd[~m], dd[m].copy()126 try:127 lin = smf.glm(base, tr, family=BIN).fit(); ll_lin = oos_loglik_binomial(ts.y1_severity, lin.predict(ts))128 h = hinge_fit(tr, "y1_severity", base, "P_main_z", BIN)129 ts["_hinge"] = np.clip(ts.P_main_z - h[1], 0, None); ll_h = oos_loglik_binomial(ts.y1_severity, h[2].predict(ts))130 sp = smf.glm(f"y1_severity ~ cr(P_main_z, df=3) + {XS} + {FE}", tr, family=BIN).fit(); ll_s = oos_loglik_binomial(ts.y1_severity, sp.predict(ts))131 rows.append({"fold": f"{rep}_{i}", "ll_linear": ll_lin, "ll_hinge": ll_h, "ll_spline": ll_s, "knot": h[1]})132 except Exception as e:133 rows.append({"fold": f"{rep}_{i}", "error": str(e)[:80]})134 r = pd.DataFrame(rows); r.to_csv(RES / "h4_seshat_folds.csv", index=False)135 ok = r.dropna(subset=["ll_linear"])136 summ = {"folds": int(len(ok)), "hinge_beats_linear_share": float((ok.ll_hinge > ok.ll_linear).mean()), "spline_beats_linear_share": float((ok.ll_spline > ok.ll_linear).mean()),137 "mean_ll_linear": float(ok.ll_linear.mean()), "mean_ll_hinge": float(ok.ll_hinge.mean()), "mean_ll_spline": float(ok.ll_spline.mean()), "median_knot": float(ok.knot.median())}138 # in-sample hinge for description139 h = hinge_fit(dd, "y1_severity", base, "P_main_z", BIN); summ["insample_knot"] = float(h[1]); summ["insample_hinge_coef"] = float(h[2].params["_hinge"]); summ["insample_P_coef"] = float(h[2].params["P_main_z"])140 return summ141142143def h5_blocks(d: pd.DataFrame):144 rows = []145 for col, name in [("macro_region", "region"), ("era", "era")]:146 for g in sorted(d[col].dropna().unique()):147 held = d[d[col] != g]; inb = d[d[col] == g]148 other = "C(era)" if col == "macro_region" else "C(macro_region)"149 try:150 try:151 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups="seshat_region"); cov = "cluster"152 except Exception:153 try:154 r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups=None); cov = "HC1 (cluster cov singular)"155 except Exception:156 r = fit_formula(f"y1_severity ~ P_main_z + {XS}", held, BIN, groups="seshat_region"); cov = "cluster, FE dropped (design singular)"157 row = coef_row(r, "P_main_z", f"H5 leave-{name}-out: {g}"); row["held_out_n"] = int(inb.y1_severity.notna().sum()); row["cov"] = cov158 if inb.P_main_z.notna().sum() >= 12:159 try:160 rb = fit_formula(f"y1_severity ~ P_main_z + {XS}", inb, BIN, groups="seshat_region"); row["within_block_coef"] = rb.params["P_main_z"]; row["within_block_se"] = rb.bse["P_main_z"]161 except Exception:162 try:163 rb = fit_formula("y1_severity ~ P_main_z", inb, BIN, groups=None); row["within_block_coef"] = rb.params["P_main_z"]; row["within_block_se"] = rb.bse["P_main_z"]; row["within_block_note"] = "P only, HC1"164 except Exception as e2:165 row["within_block_note"] = str(e2)[:60]166 rows.append(row)167 except Exception as e:168 rows.append({"model": f"H5 {name} {g}", "error": str(e)[:80]})169 out = pd.DataFrame(rows); out.to_csv(RES / "h5_seshat_blocks.csv", index=False)170 return out171172173def y2_cox_and_y3(dp: pd.DataFrame):174 rows = []175 d = dp.dropna(subset=["P_main", "duration", "macro_region"]).copy()176 d["event"] = (~d.censored).astype(int); d["P_z"] = (d.P_main - d.P_main.mean()) / d.P_main.std()177 for c in ["X_warfare", "X_prod_infra", "X_military", "D_pc1", "K_mean", "complexity_pc1"]:178 d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std()179 d["era_num"] = d.midpoint180 specs = {"Y2-1 Cox: P + strata(macro-region)": ["P_z"], "Y2-2 Cox: P + X + midpoint + strata": ["P_z", "X_warfare_z", "X_prod_infra_z", "midpoint"],181 "Y2-3 Cox: D + K + X + midpoint + strata": ["D_pc1_z", "K_mean_z", "X_warfare_z", "X_prod_infra_z", "midpoint"], "Y2-4 Cox: complexity PC1 + X + strata": ["complexity_pc1_z", "X_warfare_z", "X_prod_infra_z", "midpoint"]}182 for label, cols in specs.items():183 dd = d.dropna(subset=cols + ["seshat_region"])[cols + ["duration", "event", "macro_region", "seshat_region"]].copy()184 try:185 cph = CoxPHFitter(); cph.fit(dd, duration_col="duration", event_col="event", strata=["macro_region"], cluster_col="seshat_region", robust=True)186 for t in cols[:2] if "Y2-3" in label else cols[:1]:187 s = cph.summary.loc[t]188 rows.append({"model": label, "term": t, "coef": s["coef"], "se": s["se(coef)"], "ci_lo": s["coef lower 95%"], "ci_hi": s["coef upper 95%"], "p": s["p"], "n": len(dd), "events": int(dd.event.sum()), "hr": np.exp(s["coef"])})189 if "Y2-2" in label:190 ph = proportional_hazard_test(cph, dd, time_transform="rank"); rows[-1]["schoenfeld_p_P"] = float(ph.summary.loc["P_z", "p"])191 except Exception as e:192 rows.append({"model": label, "error": str(e)[:100]})193 # Y3 discontinuous succession194 meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "y3_discontinuous"]]195 d3 = dp.merge(meta, on="polity_id").dropna(subset=["y3_discontinuous", "P_main"]).copy(); d3["P_z"] = (d3.P_main - d3.P_main.mean()) / d3.P_main.std()196 for label, f in {"Y3-1 logit: discontinuous succession ~ P + FE": "y3_discontinuous ~ P_z + C(macro_region) + C(era)", "Y3-2 logit: + X": "y3_discontinuous ~ P_z + X_warfare + X_prod_infra + C(macro_region) + C(era)"}.items():197 try:198 r = fit_formula(f, d3, BIN, groups="seshat_region"); rows.append(coef_row(r, "P_z", label, {"events": int(r._data_used.y3_discontinuous.sum())}))199 except Exception as e:200 rows.append({"model": label, "error": str(e)[:100]})201 out = pd.DataFrame(rows); out.to_csv(RES / "h2_seshat_y2_y3.csv", index=False)202 return out203204205def main():206 d = load_crisis(); dp = pd.read_parquet(PROC / "seshat_indices.parquet")207 d.to_parquet(PROC / "seshat_h2_sample.parquet", index=False)208 prim = primary_models(d); prim.to_csv(RES / "h2_seshat_primary.csv", index=False)209 rob = robustness(d); rob.to_csv(RES / "h2_seshat_robustness.csv", index=False)210 wins, summ = out_of_sample(d)211 h4 = h4_nonlinear(d); (RES / "h4_seshat_summary.json").write_text(json.dumps(h4, indent=2))212 h5 = h5_blocks(d)213 y23 = y2_cox_and_y3(dp)214 pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70)215 cols = ["model", "term", "coef", "se", "p", "n"]216 print("=== PRIMARY (Y1 fractional logit, cluster region) ===\n", prim[cols].round(3).to_string(index=False))217 print("\n=== ROBUSTNESS ===\n", rob[cols].round(3).to_string(index=False))218 print("\n=== OOS wins of A2 over baselines ===\n", wins.round(3).to_string(index=False))219 print("\n=== Y1b OOS summary ===\n", summ.round(3).to_string(index=False))220 print("\n=== H4 ===", json.dumps(h4, indent=1))221 print("\n=== H5 ===\n", h5[[c for c in h5.columns if c in ("model", "coef", "se", "p", "n", "held_out_n", "within_block_coef", "within_block_se")]].round(3).to_string(index=False))222 print("\n=== Y2 Cox / Y3 ===\n", y23[[c for c in y23.columns if c in ("model", "term", "coef", "se", "p", "n", "events", "hr", "schoenfeld_p_P", "error")]].round(3).to_string(index=False))223224225if __name__ == "__main__":226 main()227