1"""H2 (and H4/H5) on the modern panel (V-Dem v16 + Maddison 2023), pre-registration §5.4.2Primary outcome Y456 = any instability event in (t, t+5]; pooled logit with region and year FE, country-clustered SE.3Usage: python -m wp35.econometrics.h2_modern4"""5from __future__ import annotations6import json7from pathlib import Path8import numpy as np9import pandas as pd10import statsmodels.api as sm11import statsmodels.formula.api as smf12from wp35.econometrics.common import fit_formula, coef_row, binary_metrics, oos_loglik_binomial, hinge_fit, fold_win_share1314ROOT = Path(__file__).resolve().parents[3]15PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical"16BIN = sm.families.Binomial()17FE = "C(e_regionpol_6C) + C(year)"18XS = "log_gdppc + v2x_polyarchy + t_since + t_since2 + t_since3"192021def window_outcome(df: pd.DataFrame, col: str, h: int, forward: bool = True) -> pd.Series:22 """any(event in (t, t+h]) (forward) or any(event in [t-h, t-1]) (backward); NaN if no coded year in the window."""23 out = pd.Series(np.nan, index=df.index)24 for cid, g in df.groupby("country_text_id"):25 s = g.set_index("year")[col]26 yrs = g.year.values27 vals = []28 for t in yrs:29 rng = range(t + 1, t + h + 1) if forward else range(t - h, t)30 w = s.reindex(rng)31 vals.append(np.nan if w.notna().sum() == 0 else float((w == 1).any()))32 out.loc[g.index] = vals33 return out343536def build_sample(h: int = 5):37 p = pd.read_parquet(PROC / "modern_indices.parquet").sort_values(["country_text_id", "year"]).reset_index(drop=True)38 p["any_event"] = p[["y4_civilwar_onset", "y5_coup", "y6_irregular_regime_end"]].max(axis=1, skipna=True)39 p.loc[p[["y4_civilwar_onset", "y5_coup", "y6_irregular_regime_end"]].isna().all(axis=1), "any_event"] = np.nan40 for col, name in [("any_event", "y456"), ("y4_civilwar_onset", "y4"), ("y5_coup", "y5"), ("y6_irregular_regime_end", "y6")]:41 p[f"{name}_f{h}"] = window_outcome(p, col, h, forward=True)42 p["y456_b5"] = window_outcome(p, "any_event", 5, forward=False) # persistence / negative control43 # time since last event (years), cubic polynomial (Carter & Signorino 2010)44 p["t_since"] = np.nan45 for cid, g in p.groupby("country_text_id"):46 last = None; vals = []47 for y, e in zip(g.year, g.any_event):48 vals.append(np.nan if last is None else y - last)49 if e == 1:50 last = y51 p.loc[g.index, "t_since"] = vals52 p["t_since"] = p.t_since.fillna(p.t_since.max()) / 100.0 # no prior event observed → treated as long spell (max observed), flagged in the paper53 p["t_since2"] = p.t_since ** 2; p["t_since3"] = p.t_since ** 354 for c in ["P_modern", "D_modern", "K_modern", "P_modern_pc1", "P_modern_within_year", "P_modern_area"]:55 p[c + "_z"] = (p[c] - p[c].mean()) / p[c].std()56 p["P_placebo_z"] = ((p.D_modern - p.D_modern.mean()) / p.D_modern.std() - (p.v2x_polyarchy - p.v2x_polyarchy.mean()) / p.v2x_polyarchy.std())57 p["period"] = pd.cut(p.year, [1700, 1899, 1945, 1989, 2100], labels=["pre-1900", "1900-1945", "1946-1989", "1990+"])58 return p596061def primary(d: pd.DataFrame, y: str = "y456_f5"):62 rows = []63 specs = {64 "C1 primary: Y456 ~ P + X + region FE + year FE": (f"{y} ~ P_modern_z + {XS} + {FE}", "P_modern_z"),65 "C2: Y456 ~ P + FE only": (f"{y} ~ P_modern_z + {FE}", "P_modern_z"),66 "C3: Y456 ~ P + K + X + FE (P given K)": (f"{y} ~ P_modern_z + K_modern_z + {XS} + {FE}", "P_modern_z"),67 "B1 baseline population only": (f"{y} ~ D_modern_z + {XS} + {FE}", "D_modern_z"),68 "B2 baseline capacity only": (f"{y} ~ K_modern_z + {XS} + {FE}", "K_modern_z"),69 "B3 baseline development (gdppc) only": (f"{y} ~ log_gdppc + v2x_polyarchy + t_since + t_since2 + t_since3 + {FE}", "log_gdppc"),70 "B4 baseline persistence": (f"{y} ~ y456_b5 + {XS} + {FE}", "y456_b5"),71 "B5 D and K separately": (f"{y} ~ D_modern_z + K_modern_z + {XS} + {FE}", "D_modern_z"),72 "Placebo: P with polyarchy in place of K (polyarchy not separately controlled)": (f"{y} ~ P_placebo_z + log_gdppc + t_since + t_since2 + t_since3 + {FE}", "P_placebo_z"),73 "Placebo comparison: P_main with the same controls as the placebo": (f"{y} ~ P_modern_z + log_gdppc + t_since + t_since2 + t_since3 + {FE}", "P_modern_z"),74 "Negative control: past events (t-5..t-1) ~ P_t + gdppc + polyarchy + FE (no duration terms)": (f"y456_b5 ~ P_modern_z + log_gdppc + v2x_polyarchy + {FE}", "P_modern_z"),75 "Forward comparison for the negative control: future events ~ P_t + gdppc + polyarchy + FE": (f"{y} ~ P_modern_z + log_gdppc + v2x_polyarchy + {FE}", "P_modern_z"),76 }77 for label, (f, term) in specs.items():78 try:79 r = fit_formula(f, d, BIN, groups="country_text_id"); rows.append(coef_row(r, term, label, {"events": int(r._data_used[f.split(" ~")[0]].sum())}))80 if "B5" in label:81 rows.append(coef_row(r, "K_modern_z", label))82 except Exception as e:83 rows.append({"model": label, "error": str(e)[:100]})84 # country fixed effects: conditional logit via statsmodels (drops never-event countries)85 try:86 from statsmodels.discrete.conditional_models import ConditionalLogit87 dd = d.dropna(subset=[y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since"]).copy()88 dd["decade"] = (dd.year // 10) * 1089 X = pd.concat([dd[["P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "t_since2", "t_since3"]], pd.get_dummies(dd.decade, prefix="dec", drop_first=True, dtype=float)], axis=1)90 cl = ConditionalLogit(dd[y].values, X.values, groups=pd.factorize(dd.country_text_id)[0]).fit(disp=0, maxiter=200)91 b, se = cl.params[0], cl.bse[0]92 rows.append({"model": "C4: conditional logit, country FE + decade FE", "term": "P_modern_z", "coef": b, "se": se, "ci_lo": b - 1.96 * se, "ci_hi": b + 1.96 * se, "p": cl.pvalues[0], "n": len(dd), "events": int(dd[y].sum())})93 except Exception as e:94 rows.append({"model": "C4 conditional logit", "error": str(e)[:120]})95 # LPM two-way FE96 try:97 dd = d.dropna(subset=[y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since"]).copy()98 for c in [y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "t_since2", "t_since3"]:99 dd[c + "_dm"] = dd[c] - dd.groupby("country_text_id")[c].transform("mean")100 r = fit_formula(f"{y}_dm ~ P_modern_z_dm + log_gdppc_dm + v2x_polyarchy_dm + t_since_dm + t_since2_dm + t_since3_dm + C(year)", dd, None, groups="country_text_id")101 rows.append(coef_row(r, "P_modern_z_dm", "C5: linear probability, country + year FE", {"events": int(dd[y].sum())}))102 except Exception as e:103 rows.append({"model": "C5 LPM", "error": str(e)[:100]})104 return pd.DataFrame(rows)105106107def oos(d: pd.DataFrame, y: str = "y456_f5"):108 formulas = {"C1 P": f"{y} ~ P_modern_z + {XS} + {FE}", "B0 FE+X": f"{y} ~ {XS} + {FE}", "B1 D": f"{y} ~ D_modern_z + {XS} + {FE}", "B2 K": f"{y} ~ K_modern_z + {XS} + {FE}",109 "B4 persistence": f"{y} ~ y456_b5 + {XS} + {FE}", "B5 D+K": f"{y} ~ D_modern_z + K_modern_z + {XS} + {FE}"}110 rows = []111 need = [y, "P_modern_z", "D_modern_z", "K_modern_z", "y456_b5", "log_gdppc", "v2x_polyarchy", "t_since", "e_regionpol_6C", "year"]112 dd = d.dropna(subset=need).copy()113 # forward chaining (year FE replaced by decade FE so that test years have a fitted level: pre-registered deviation logged as A-001)114 fc = {k: v.replace("C(year)", "C(decade)") for k, v in formulas.items()}115 dd["decade"] = (dd.year // 10) * 10116 splits = [("train<=1945,test1946-1975", dd.year <= 1945, dd.year.between(1946, 1975)), ("train<=1975,test1976-2005", dd.year <= 1975, dd.year.between(1976, 2005)), ("train<=2005,test>=2006", dd.year <= 2005, dd.year >= 2006)]117 for name, f in fc.items():118 for sname, tr, te in splits:119 try:120 m = smf.glm(f, dd[tr], family=BIN).fit()121 te_d = dd[te].copy(); te_d["decade"] = np.minimum(te_d.decade, dd[tr].decade.max()) # unseen decades mapped to the last training decade122 p = m.predict(te_d); r = binary_metrics(te_d[y], p); r.update({"model": name, "fold": sname, "scheme": "forward"}); rows.append(r)123 except Exception as e:124 rows.append({"model": name, "fold": sname, "scheme": "forward", "error": str(e)[:80]})125 # leave-one-region-out (region FE removed)126 lr = {k: v.replace("C(e_regionpol_6C) + ", "") for k, v in formulas.items()}127 for name, f in lr.items():128 for reg in sorted(dd.e_regionpol_6C.unique()):129 try:130 m = smf.glm(f, dd[dd.e_regionpol_6C != reg], family=BIN).fit(); te_d = dd[dd.e_regionpol_6C == reg]131 p = m.predict(te_d); r = binary_metrics(te_d[y], p); r.update({"model": name, "fold": f"loro_{int(reg)}", "scheme": "loro"}); rows.append(r)132 except Exception as e:133 rows.append({"model": name, "fold": f"loro_{int(reg)}", "scheme": "loro", "error": str(e)[:80]})134 cv = pd.DataFrame(rows); cv.to_csv(RES / "h2_modern_cv_folds.csv", index=False)135 wins = []136 for scheme in ["forward", "loro"]:137 c = cv[cv.scheme == scheme]138 for b in ["B0 FE+X", "B1 D", "B2 K", "B4 persistence", "B5 D+K"]:139 for metric in ["oos_loglik", "auc"]:140 w = fold_win_share(c, "C1 P", b, metric); w.update({"scheme": scheme, "metric": metric}); wins.append(w)141 wins = pd.DataFrame(wins); wins.to_csv(RES / "h2_modern_cv_wins.csv", index=False)142 summ = cv.groupby(["scheme", "model"])[["oos_loglik", "auc", "brier", "calibration_slope"]].mean().reset_index(); summ.to_csv(RES / "h2_modern_cv_summary.csv", index=False)143 return wins, summ144145146def h4(d: pd.DataFrame, y: str = "y456_f5"):147 base = f"{y} ~ P_modern_z + {XS} + C(e_regionpol_6C) + C(decade)"148 dd = d.dropna(subset=[y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "e_regionpol_6C"]).copy(); dd["decade"] = (dd.year // 10) * 10149 rng = np.random.default_rng(20261101); rows = []150 countries = dd.country_text_id.unique()151 for rep in range(3):152 perm = rng.permutation(countries); parts = np.array_split(perm, 5) # country-blocked folds153 for i, te_c in enumerate(parts):154 te = dd.country_text_id.isin(te_c); tr, ts = dd[~te], dd[te].copy()155 try:156 lin = smf.glm(base, tr, family=BIN).fit(); ll_lin = oos_loglik_binomial(ts[y], lin.predict(ts))157 h = hinge_fit(tr, y, base, "P_modern_z", BIN, grid=np.quantile(tr.P_modern_z, np.linspace(0.2, 0.8, 7)))158 ts["_hinge"] = np.clip(ts.P_modern_z - h[1], 0, None); ll_h = oos_loglik_binomial(ts[y], h[2].predict(ts))159 sp = smf.glm(base.replace("P_modern_z", "cr(P_modern_z, df=3)"), tr, family=BIN).fit(); ll_s = oos_loglik_binomial(ts[y], sp.predict(ts))160 rows.append({"fold": f"{rep}_{i}", "ll_linear": ll_lin, "ll_hinge": ll_h, "ll_spline": ll_s, "knot": h[1]})161 except Exception as e:162 rows.append({"fold": f"{rep}_{i}", "error": str(e)[:80]})163 r = pd.DataFrame(rows); r.to_csv(RES / "h4_modern_folds.csv", index=False); ok = r.dropna(subset=["ll_linear"])164 # in-sample shape (for the figure): predicted probability vs P at covariate means, linear vs spline165 sp = smf.glm(base.replace("P_modern_z", "cr(P_modern_z, df=3)"), dd, family=BIN).fit(); lin = smf.glm(base, dd, family=BIN).fit()166 grid = np.linspace(dd.P_modern_z.quantile(0.02), dd.P_modern_z.quantile(0.98), 40)167 ref = dd.iloc[[0]].copy()168 for c in ["log_gdppc", "v2x_polyarchy", "t_since", "t_since2", "t_since3"]:169 ref[c] = dd[c].mean()170 ref["e_regionpol_6C"] = dd.e_regionpol_6C.mode()[0]; ref["decade"] = 1970171 g = pd.concat([ref] * len(grid), ignore_index=True); g["P_modern_z"] = grid172 pd.DataFrame({"P_z": grid, "p_linear": lin.predict(g), "p_spline": sp.predict(g)}).to_csv(RES / "h4_modern_curve.csv", index=False)173 return {"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()),174 "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())}175176177def h5(d: pd.DataFrame, y: str = "y456_f5"):178 rows = []179 for col, name in [("e_regionpol_6C", "region"), ("period", "period")]:180 for g in sorted(d[col].dropna().unique()):181 held = d[d[col] != g]; inb = d[d[col] == g]182 fe = "C(year)" if col == "e_regionpol_6C" else "C(e_regionpol_6C) + C(year)"183 try:184 r = fit_formula(f"{y} ~ P_modern_z + {XS} + {fe}", held, BIN, groups="country_text_id"); row = coef_row(r, "P_modern_z", f"H5 leave-{name}-out: {g}")185 rb = fit_formula(f"{y} ~ P_modern_z + {XS} + " + ("C(year)" if col == "e_regionpol_6C" else "C(e_regionpol_6C)"), inb, BIN, groups="country_text_id")186 row.update({"within_block_coef": rb.params["P_modern_z"], "within_block_se": rb.bse["P_modern_z"], "within_block_n": rb._n_used, "within_block_events": int(rb._data_used[y].sum())})187 rows.append(row)188 except Exception as e:189 rows.append({"model": f"H5 {name} {g}", "error": str(e)[:80]})190 out = pd.DataFrame(rows); out.to_csv(RES / "h5_modern_blocks.csv", index=False); return out191192193def robustness(d: pd.DataFrame):194 rows = []195 for h in [1, 3, 10]:196 dh = build_sample(h) if h != 5 else d197 r = fit_formula(f"y456_f{h} ~ P_modern_z + {XS} + {FE}", dh, BIN, groups="country_text_id"); rows.append(coef_row(r, "P_modern_z", f"R: horizon h={h}", {"events": int(r._data_used[f"y456_f{h}"].sum())}))198 for comp in ["y4_f5", "y5_f5", "y6_f5"]:199 r = fit_formula(f"{comp} ~ P_modern_z + {XS} + {FE}", d, BIN, groups="country_text_id"); rows.append(coef_row(r, "P_modern_z", f"R: component {comp}", {"events": int(r._data_used[comp].sum())}))200 for v, lab in [("P_modern_pc1_z", "K = PC1"), ("P_modern_within_year_z", "standardized within year"), ("P_modern_area_z", "D = (log pop + log area)/2")]:201 r = fit_formula(f"y456_f5 ~ {v} + {XS} + {FE}", d, BIN, groups="country_text_id"); rows.append(coef_row(r, v, f"R: {lab}", {"events": int(r._data_used.y456_f5.sum())}))202 r = fit_formula(f"y456_f5 ~ P_modern_z + {XS} + {FE}", d[d.year >= 1900], BIN, groups="country_text_id"); rows.append(coef_row(r, "P_modern_z", "R: years ≥ 1900", {"events": int(r._data_used.y456_f5.sum())}))203 r = fit_formula(f"y456_f5 ~ P_modern_z + log_gdppc + t_since + t_since2 + t_since3 + {FE}", d, BIN, groups="country_text_id"); rows.append(coef_row(r, "P_modern_z", "R: polyarchy not controlled", {"events": int(r._data_used.y456_f5.sum())}))204 r = fit_formula(f"y456_f5 ~ P_modern_z + {XS} + y456_b5 + {FE}", d, BIN, groups="country_text_id"); rows.append(coef_row(r, "P_modern_z", "R: + persistence control", {"events": int(r._data_used.y456_f5.sum())}))205 # Oster / Cinelli-Hazlett style sensitivity on the LPM version of the primary206 try:207 dd = d.dropna(subset=["y456_f5", "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "e_regionpol_6C"]).copy()208 full = smf.ols(f"y456_f5 ~ P_modern_z + {XS} + {FE}", dd).fit(); short = smf.ols(f"y456_f5 ~ P_modern_z + {FE}", dd).fit()209 restricted = smf.ols(f"y456_f5 ~ {XS} + {FE}", dd).fit()210 r2_full, r2_short = full.rsquared, short.rsquared; b_full, b_short = full.params["P_modern_z"], short.params["P_modern_z"]211 rmax = min(1.0, 1.3 * r2_full); delta1 = b_full - (b_short - b_full) * (rmax - r2_full) / max(r2_full - r2_short, 1e-9)212 t = full.tvalues["P_modern_z"]; dfres = full.df_resid; r2_yd = t ** 2 / (t ** 2 + dfres)213 rv = 0.5 * (np.sqrt(r2_yd ** 2 + 4 * r2_yd) - r2_yd)214 rows.append({"model": "Sensitivity (LPM): Oster bias-adjusted β at Rmax=1.3·R² (δ=1); Cinelli–Hazlett RV", "term": "P_modern_z", "coef": delta1, "oster_beta_full": b_full, "oster_beta_short": b_short, "r2_full": r2_full, "r2_short": r2_short, "partial_r2_P": r2_yd, "robustness_value": rv, "n": len(dd)})215 except Exception as e:216 rows.append({"model": "Sensitivity", "error": str(e)[:100]})217 out = pd.DataFrame(rows); out.to_csv(RES / "h2_modern_robustness.csv", index=False); return out218219220def main():221 d = build_sample(5); d.to_parquet(PROC / "modern_h2_sample.parquet", index=False)222 prim = primary(d); prim.to_csv(RES / "h2_modern_primary.csv", index=False)223 wins, summ = oos(d)224 h4s = h4(d); (RES / "h4_modern_summary.json").write_text(json.dumps(h4s, indent=2))225 h5s = h5(d); rob = robustness(d)226 pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70)227 cols = ["model", "term", "coef", "se", "p", "n", "events", "error"]228 print("=== PRIMARY ===\n", prim[[c for c in cols if c in prim.columns]].round(3).to_string(index=False))229 print("\n=== OOS wins ===\n", wins.round(3).to_string(index=False)); print("\n=== OOS summary ===\n", summ.round(3).to_string(index=False))230 print("\n=== H4 ===", json.dumps(h4s, indent=1))231 print("\n=== H5 ===\n", h5s[[c for c in h5s.columns if c in ("model", "coef", "se", "p", "n", "within_block_coef", "within_block_se", "within_block_n", "within_block_events", "error")]].round(3).to_string(index=False))232 print("\n=== ROBUSTNESS ===\n", rob[[c for c in rob.columns if c in cols + ["robustness_value", "partial_r2_P", "oster_beta_full", "oster_beta_short"]]].round(4).to_string(index=False))233234235if __name__ == "__main__":236 main()237