"""H2 (and H4/H5) on Seshat: Y1 crisis severity (primary), Y1b severe binary, Y2 polity duration (Cox), Y3 succession. This is the first module that reads outcomes together with P. Usage: python -m wp35.econometrics.h2_seshat Everything follows research/preregistration.md §5.2–5.3, §5.5–5.7, §6. """ from __future__ import annotations import json from pathlib import Path import numpy as np import pandas as pd import statsmodels.api as sm import statsmodels.formula.api as smf from lifelines import CoxPHFitter from lifelines.statistics import proportional_hazard_test from wp35.econometrics.common import fit_formula, coef_row, cv_compare, fold_win_share, hinge_fit, binary_metrics, oos_loglik_binomial ROOT = Path(__file__).resolve().parents[3] PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical" RES.mkdir(parents=True, exist_ok=True) BIN = sm.families.Binomial() FE = "C(macro_region) + C(era)" XS = "X_warfare + X_prod_infra + epidemic_v" def load_crisis(): cr = pd.read_parquet(PROC / "seshat_crisis_indices.parquet") y = pd.read_parquet(PROC / "seshat_crisis.parquet")[["crisis_id", "y1_severity", "y1b_severe", "y_reform", "epidemic_v", "y1_n_coded"]] d = cr.merge(y, on="crisis_id") meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "centralization", "homeNGA"]] d = d.merge(meta, on="polity_id", how="left", suffixes=("", "_m")) d["homeNGA"] = d["homeNGA"].fillna(d.get("homeNGA_m")) # persistence baseline: an earlier crisis record in the same NGA d = d.sort_values("crisis_start") d["prior_crisis_nga"] = 0.0 for i, r in d.iterrows(): d.loc[i, "prior_crisis_nga"] = float(((d.homeNGA == r.homeNGA) & (d.crisis_start < r.crisis_start)).any()) d["quasi"] = d.centralization.eq("quasi-polity") 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() for c in ["P_main", "D_pc1", "K_mean", "K_irt", "complexity_pc1", "log_pop"]: d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std() return d def primary_models(d: pd.DataFrame): rows = [] specs = { "A1 primary: Y1 ~ P + FE": (f"y1_severity ~ P_main_z + {FE}", "P_main_z"), "A2: Y1 ~ P + X + FE": (f"y1_severity ~ P_main_z + {XS} + {FE}", "P_main_z"), "A3: Y1 ~ P + D + K + X + FE (P given components)": (f"y1_severity ~ P_main_z + complexity_pc1_z + {XS} + {FE}", "P_main_z"), "B1 baseline D only": (f"y1_severity ~ D_pc1_z + {XS} + {FE}", "D_pc1_z"), "B2 baseline K only": (f"y1_severity ~ K_mean_z + {XS} + {FE}", "K_mean_z"), "B3 baseline complexity PC1": (f"y1_severity ~ complexity_pc1_z + {XS} + {FE}", "complexity_pc1_z"), "B4 baseline persistence": (f"y1_severity ~ prior_crisis_nga + {XS} + {FE}", "prior_crisis_nga"), "B5 D and K separately": (f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}", "D_pc1_z"), "Placebo: P built with production infrastructure instead of K": (f"y1_severity ~ P_placebo + {XS} + {FE}", "P_placebo"), } for label, (f, term) in specs.items(): r = fit_formula(f, d, BIN, groups="seshat_region"); rows.append(coef_row(r, term, label)) if "B5" in label: rows.append(coef_row(r, "K_mean_z", label)) # binary secondary # binary secondary outcome: several region/era cells are all-severe (quasi-separation) -> region FE only, HC1 SE for label, f in {"A1b: Y1b ~ P + region FE (logit, HC1; era FE dropped: separation)": "y1b_severe ~ P_main_z + C(macro_region)", "A2b: Y1b ~ P + X + region FE (logit, HC1)": f"y1b_severe ~ P_main_z + {XS} + C(macro_region)", "A3b: Y1b ~ P (logit, HC1, no FE)": "y1b_severe ~ P_main_z"}.items(): r = fit_formula(f, d, BIN, groups=None); rows.append(coef_row(r, "P_main_z", label)) # exploratory: reforms 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")) return pd.DataFrame(rows) def robustness(d: pd.DataFrame): rows = [] 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", "P_main_lo": "D at lower bounds", "P_main_hi": "D at upper bounds", "P_main_noinf": "inferred codes dropped", "P_within_era": "standardized within era", "P_within_region": "standardized within macro-region", "P_rank": "rank-based P"} for v, lab in variants.items(): if v not in d.columns: # crisis-level file lacks some polity-level variants → merge from polity file continue dd = d.copy(); dd["_P"] = (dd[v] - dd[v].mean()) / dd[v].std() r = fit_formula(f"y1_severity ~ _P + {XS} + {FE}", dd, BIN, groups="seshat_region"); rows.append(coef_row(r, "_P", f"R: {lab}")) 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")) 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")) 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")) 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")) 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")) 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")) 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")) 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")) return pd.DataFrame(rows) def out_of_sample(d: pd.DataFrame): 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}", "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}", "B5 D+K": f"y1_severity ~ D_pc1_z + K_mean_z + {XS} + {FE}"} # era FE with leave-one-region-out can leave unseen levels; macro-region FE cannot be used when the region is held out f_loro = {k: v.replace("C(macro_region) + ", "") for k, v in formulas.items()} loro = cv_compare(d, f_loro, "y1_severity", BIN, fold_col="macro_region", fractional=True) kf = cv_compare(d, formulas, "y1_severity", BIN, n_repeats=10, k=5, fractional=True) loro["scheme"] = "leave-one-macro-region-out"; kf["scheme"] = "10x5-fold" cv = pd.concat([loro, kf]); cv.to_csv(RES / "h2_seshat_cv_folds.csv", index=False) wins = [] for scheme, c in [("loro", loro), ("kfold", kf)]: for b in ["B0 FE only", "B1 D only", "B2 K only", "B3 complexity", "B4 persistence", "B5 D+K"]: w = fold_win_share(c, "A2 P", b); w["scheme"] = scheme; wins.append(w) wins = pd.DataFrame(wins); wins.to_csv(RES / "h2_seshat_cv_wins.csv", index=False) # binary Y1b out-of-sample metrics fb = {k: v.replace("y1_severity", "y1b_severe") for k, v in formulas.items()} kfb = cv_compare(d, fb, "y1b_severe", BIN, n_repeats=10, k=5) kfb.to_csv(RES / "h2_seshat_cv_folds_binary.csv", index=False) summ = kfb.groupby("model")[["oos_loglik", "auc", "brier"]].mean().reset_index(); summ.to_csv(RES / "h2_seshat_cv_binary_summary.csv", index=False) return wins, summ def h4_nonlinear(d: pd.DataFrame): """Linear vs hinge vs spline, out of sample (10x5-fold); hinge knot chosen on the training fold only.""" rng = np.random.default_rng(20261101) base = f"y1_severity ~ P_main_z + {XS} + {FE}" dd = d.dropna(subset=["y1_severity", "P_main_z", "X_warfare", "X_prod_infra", "epidemic_v", "macro_region", "era"]).reset_index(drop=True) rows = [] for rep in range(10): perm = rng.permutation(len(dd)); parts = np.array_split(perm, 5) for i, te in enumerate(parts): m = np.zeros(len(dd), bool); m[te] = True; tr, ts = dd[~m], dd[m].copy() try: lin = smf.glm(base, tr, family=BIN).fit(); ll_lin = oos_loglik_binomial(ts.y1_severity, lin.predict(ts)) h = hinge_fit(tr, "y1_severity", base, "P_main_z", BIN) ts["_hinge"] = np.clip(ts.P_main_z - h[1], 0, None); ll_h = oos_loglik_binomial(ts.y1_severity, h[2].predict(ts)) 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)) rows.append({"fold": f"{rep}_{i}", "ll_linear": ll_lin, "ll_hinge": ll_h, "ll_spline": ll_s, "knot": h[1]}) except Exception as e: rows.append({"fold": f"{rep}_{i}", "error": str(e)[:80]}) r = pd.DataFrame(rows); r.to_csv(RES / "h4_seshat_folds.csv", index=False) ok = r.dropna(subset=["ll_linear"]) 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()), "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())} # in-sample hinge for description 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"]) return summ def h5_blocks(d: pd.DataFrame): rows = [] for col, name in [("macro_region", "region"), ("era", "era")]: for g in sorted(d[col].dropna().unique()): held = d[d[col] != g]; inb = d[d[col] == g] other = "C(era)" if col == "macro_region" else "C(macro_region)" try: try: r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups="seshat_region"); cov = "cluster" except Exception: try: r = fit_formula(f"y1_severity ~ P_main_z + {XS} + {other}", held, BIN, groups=None); cov = "HC1 (cluster cov singular)" except Exception: r = fit_formula(f"y1_severity ~ P_main_z + {XS}", held, BIN, groups="seshat_region"); cov = "cluster, FE dropped (design singular)" 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"] = cov if inb.P_main_z.notna().sum() >= 12: try: 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"] except Exception: try: 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" except Exception as e2: row["within_block_note"] = str(e2)[:60] rows.append(row) except Exception as e: rows.append({"model": f"H5 {name} {g}", "error": str(e)[:80]}) out = pd.DataFrame(rows); out.to_csv(RES / "h5_seshat_blocks.csv", index=False) return out def y2_cox_and_y3(dp: pd.DataFrame): rows = [] d = dp.dropna(subset=["P_main", "duration", "macro_region"]).copy() d["event"] = (~d.censored).astype(int); d["P_z"] = (d.P_main - d.P_main.mean()) / d.P_main.std() for c in ["X_warfare", "X_prod_infra", "X_military", "D_pc1", "K_mean", "complexity_pc1"]: d[c + "_z"] = (d[c] - d[c].mean()) / d[c].std() d["era_num"] = d.midpoint 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"], "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"]} for label, cols in specs.items(): dd = d.dropna(subset=cols + ["seshat_region"])[cols + ["duration", "event", "macro_region", "seshat_region"]].copy() try: cph = CoxPHFitter(); cph.fit(dd, duration_col="duration", event_col="event", strata=["macro_region"], cluster_col="seshat_region", robust=True) for t in cols[:2] if "Y2-3" in label else cols[:1]: s = cph.summary.loc[t] 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"])}) if "Y2-2" in label: ph = proportional_hazard_test(cph, dd, time_transform="rank"); rows[-1]["schoenfeld_p_P"] = float(ph.summary.loc["P_z", "p"]) except Exception as e: rows.append({"model": label, "error": str(e)[:100]}) # Y3 discontinuous succession meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet")[["polity_id", "y3_discontinuous"]] 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() 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(): try: 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())})) except Exception as e: rows.append({"model": label, "error": str(e)[:100]}) out = pd.DataFrame(rows); out.to_csv(RES / "h2_seshat_y2_y3.csv", index=False) return out def main(): d = load_crisis(); dp = pd.read_parquet(PROC / "seshat_indices.parquet") d.to_parquet(PROC / "seshat_h2_sample.parquet", index=False) prim = primary_models(d); prim.to_csv(RES / "h2_seshat_primary.csv", index=False) rob = robustness(d); rob.to_csv(RES / "h2_seshat_robustness.csv", index=False) wins, summ = out_of_sample(d) h4 = h4_nonlinear(d); (RES / "h4_seshat_summary.json").write_text(json.dumps(h4, indent=2)) h5 = h5_blocks(d) y23 = y2_cox_and_y3(dp) pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70) cols = ["model", "term", "coef", "se", "p", "n"] print("=== PRIMARY (Y1 fractional logit, cluster region) ===\n", prim[cols].round(3).to_string(index=False)) print("\n=== ROBUSTNESS ===\n", rob[cols].round(3).to_string(index=False)) print("\n=== OOS wins of A2 over baselines ===\n", wins.round(3).to_string(index=False)) print("\n=== Y1b OOS summary ===\n", summ.round(3).to_string(index=False)) print("\n=== H4 ===", json.dumps(h4, indent=1)) 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)) 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)) if __name__ == "__main__": main()