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 + no_prior_event"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 # no prior event observed: spell measured from the country's first observed year (left-censored), with an indicator (A-004)53 first = p.groupby("country_text_id").year.transform("min")54 p["no_prior_event"] = p.t_since.isna().astype(float)55 p["t_since"] = p.t_since.fillna(p.year - first + 1) / 100.056 p["t_since2"] = p.t_since ** 2; p["t_since3"] = p.t_since ** 357 for c in ["P_modern", "D_modern", "K_modern", "P_modern_pc1", "P_modern_within_year", "P_modern_area"]:58 p[c + "_z"] = (p[c] - p[c].mean()) / p[c].std()59 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())60 p["period"] = pd.cut(p.year, [1700, 1899, 1945, 1989, 2100], labels=["pre-1900", "1900-1945", "1946-1989", "1990+"])61 return p626364def primary(d: pd.DataFrame, y: str = "y456_f5"):65 rows = []66 specs = {67 "C1 primary: Y456 ~ P + X + region FE + year FE": (f"{y} ~ P_modern_z + {XS} + {FE}", "P_modern_z"),68 "C2: Y456 ~ P + FE only": (f"{y} ~ P_modern_z + {FE}", "P_modern_z"),69 "C3: Y456 ~ P + K + X + FE (reparametrization of B5: beta_P = beta_D)": (f"{y} ~ P_modern_z + K_modern_z + {XS} + {FE}", "P_modern_z"),70 "B1 baseline population only": (f"{y} ~ D_modern_z + {XS} + {FE}", "D_modern_z"),71 "B2 baseline capacity only": (f"{y} ~ K_modern_z + {XS} + {FE}", "K_modern_z"),72 "B3 baseline development (gdppc) only": (f"{y} ~ log_gdppc + v2x_polyarchy + t_since + t_since2 + t_since3 + {FE}", "log_gdppc"),73 "B4 baseline persistence": (f"{y} ~ y456_b5 + {XS} + {FE}", "y456_b5"),74 "B5 D and K separately": (f"{y} ~ D_modern_z + K_modern_z + {XS} + {FE}", "D_modern_z"),75 "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"),76 "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"),77 "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"),78 "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"),79 }80 for label, (f, term) in specs.items():81 try:82 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())}))83 if label.startswith("B5"):84 rows.append(coef_row(r, "K_modern_z", label))85 from wp35.econometrics.h2_seshat import wald_sum_zero86 wt = wald_sum_zero(r, "D_modern_z", "K_modern_z"); rows.append({"model": "B5 Wald test of beta_D + beta_K = 0 (restriction implied by P)", "term": "D+K", "coef": wt["sum_coef"], "se": wt["sum_se"], "p": wt["wald_p"], "n": r._n_used})87 except Exception as e:88 rows.append({"model": label, "error": str(e)[:100]})89 # country fixed effects: conditional logit via statsmodels (drops never-event countries)90 try:91 from statsmodels.discrete.conditional_models import ConditionalLogit92 dd = d.dropna(subset=[y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since"]).copy()93 dd["decade"] = (dd.year // 10) * 1094 X = pd.concat([dd[["P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "t_since2", "t_since3", "no_prior_event"]], pd.get_dummies(dd.decade, prefix="dec", drop_first=True, dtype=float)], axis=1)95 cl = ConditionalLogit(dd[y].values, X.values, groups=pd.factorize(dd.country_text_id)[0]).fit(disp=0, maxiter=200)96 b, se = cl.params[0], cl.bse[0]97 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())})98 except Exception as e:99 rows.append({"model": "C4 conditional logit", "error": str(e)[:120]})100 # LPM two-way FE101 try:102 dd = d.dropna(subset=[y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since"]).copy()103 for c in [y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "t_since2", "t_since3", "no_prior_event"]:104 dd[c + "_dm"] = dd[c] - dd.groupby("country_text_id")[c].transform("mean")105 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 + no_prior_event_dm + C(year)", dd, None, groups="country_text_id")106 rows.append(coef_row(r, "P_modern_z_dm", "C5: linear probability, country + year FE", {"events": int(dd[y].sum())}))107 except Exception as e:108 rows.append({"model": "C5 LPM", "error": str(e)[:100]})109 return pd.DataFrame(rows)110111112def oos(d: pd.DataFrame, y: str = "y456_f5"):113 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}",114 "B4 persistence": f"{y} ~ y456_b5 + {XS} + {FE}", "B5 D+K": f"{y} ~ D_modern_z + K_modern_z + {XS} + {FE}"}115 rows = []; preds = []116 need = [y, "P_modern_z", "D_modern_z", "K_modern_z", "y456_b5", "log_gdppc", "v2x_polyarchy", "t_since", "e_regionpol_6C", "year"]117 dd = d.dropna(subset=need).copy()118 # forward chaining (year FE replaced by decade FE so that test years have a fitted level: pre-registered deviation logged as A-001)119 fc = {k: v.replace("C(year)", "C(decade)") for k, v in formulas.items()}120 dd["decade"] = (dd.year // 10) * 10121 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)]122 for name, f in fc.items():123 for sname, tr, te in splits:124 try:125 m = smf.glm(f, dd[tr], family=BIN).fit()126 te_d = dd[te].copy(); te_d["decade"] = np.minimum(te_d.decade, dd[tr].decade.max()) # unseen decades mapped to the last training decade127 seen = set(dd[tr].e_regionpol_6C.unique()); te_d.loc[~te_d.e_regionpol_6C.isin(seen), "e_regionpol_6C"] = dd[tr].e_regionpol_6C.mode().iloc[0] # unseen regions (A-004)128 p = m.predict(te_d); r = binary_metrics(te_d[y], p); r.update({"model": name, "fold": sname, "scheme": "forward"}); rows.append(r)129 preds.append(pd.DataFrame({"scheme": "forward", "fold": sname, "model": name, "y": te_d[y].values, "p": np.asarray(p)}))130 except Exception as e:131 rows.append({"model": name, "fold": sname, "scheme": "forward", "error": str(e)[:80]})132 # leave-one-region-out (region FE removed)133 lr = {k: v.replace("C(e_regionpol_6C) + ", "") for k, v in formulas.items()}134 for name, f in lr.items():135 for reg in sorted(dd.e_regionpol_6C.unique()):136 try:137 m = smf.glm(f, dd[dd.e_regionpol_6C != reg], family=BIN).fit(); te_d = dd[dd.e_regionpol_6C == reg]138 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)139 preds.append(pd.DataFrame({"scheme": "loro", "fold": f"loro_{int(reg)}", "model": name, "y": te_d[y].values, "p": np.asarray(p)}))140 except Exception as e:141 rows.append({"model": name, "fold": f"loro_{int(reg)}", "scheme": "loro", "error": str(e)[:80]})142 cv = pd.DataFrame(rows); cv.to_csv(RES / "h2_modern_cv_folds.csv", index=False)143 pd.concat(preds).to_parquet(RES / "h2_modern_cv_predictions.parquet", index=False)144 wins = []145 for scheme in ["forward", "loro"]:146 c = cv[cv.scheme == scheme]147 for b in ["B0 FE+X", "B1 D", "B2 K", "B4 persistence", "B5 D+K"]:148 for metric in ["oos_loglik", "auc"]:149 w = fold_win_share(c, "C1 P", b, metric); w.update({"scheme": scheme, "metric": metric}); wins.append(w)150 wins = pd.DataFrame(wins); wins.to_csv(RES / "h2_modern_cv_wins.csv", index=False)151 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)152 return wins, summ153154155def h4(d: pd.DataFrame, y: str = "y456_f5"):156 base = f"{y} ~ P_modern_z + {XS} + C(e_regionpol_6C) + C(decade)"157 from wp35.econometrics.common import map_unseen_levels158 dd = d.dropna(subset=[y, "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "e_regionpol_6C"]).copy(); dd["decade"] = (dd.year // 10) * 10159 rng = np.random.default_rng(20261101); rows = []160 countries = dd.country_text_id.unique()161 for rep in range(3):162 perm = rng.permutation(countries); parts = np.array_split(perm, 5) # country-blocked folds163 for i, te_c in enumerate(parts):164 te = dd.country_text_id.isin(te_c); tr, ts = dd[~te], map_unseen_levels(dd[~te], dd[te], base)165 try:166 lin = smf.glm(base, tr, family=BIN).fit(); ll_lin = oos_loglik_binomial(ts[y], lin.predict(ts))167 h = hinge_fit(tr, y, base, "P_modern_z", BIN, grid=np.quantile(tr.P_modern_z, np.linspace(0.2, 0.8, 7)))168 ts["_hinge"] = np.clip(ts.P_modern_z - h[1], 0, None); ll_h = oos_loglik_binomial(ts[y], h[2].predict(ts))169 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))170 rows.append({"fold": f"{rep}_{i}", "ll_linear": ll_lin, "ll_hinge": ll_h, "ll_spline": ll_s, "knot": h[1]})171 except Exception as e:172 rows.append({"fold": f"{rep}_{i}", "error": str(e)[:80]})173 r = pd.DataFrame(rows); r.to_csv(RES / "h4_modern_folds.csv", index=False); ok = r.dropna(subset=["ll_linear"])174 # in-sample shape (for the figure): predicted probability vs P at covariate means, linear vs spline175 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()176 grid = np.linspace(dd.P_modern_z.quantile(0.02), dd.P_modern_z.quantile(0.98), 40)177 ref = dd.iloc[[0]].copy()178 for c in ["log_gdppc", "v2x_polyarchy", "t_since", "t_since2", "t_since3"]:179 ref[c] = dd[c].mean()180 ref["e_regionpol_6C"] = dd.e_regionpol_6C.mode()[0]; ref["decade"] = 1970181 g = pd.concat([ref] * len(grid), ignore_index=True); g["P_modern_z"] = grid182 pd.DataFrame({"P_z": grid, "p_linear": lin.predict(g), "p_spline": sp.predict(g)}).to_csv(RES / "h4_modern_curve.csv", index=False)183 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()),184 "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())}185186187def h5(d: pd.DataFrame, y: str = "y456_f5"):188 rows = []189 for col, name in [("e_regionpol_6C", "region"), ("period", "period")]:190 for g in sorted(d[col].dropna().unique()):191 held = d[d[col] != g]; inb = d[d[col] == g]192 fe = "C(year)" if col == "e_regionpol_6C" else "C(e_regionpol_6C) + C(year)"193 try:194 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}")195 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")196 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())})197 rows.append(row)198 except Exception as e:199 rows.append({"model": f"H5 {name} {g}", "error": str(e)[:80]})200 out = pd.DataFrame(rows); out.to_csv(RES / "h5_modern_blocks.csv", index=False); return out201202203def robustness(d: pd.DataFrame):204 rows = []205 for h in [1, 3, 10]:206 dh = build_sample(h) if h != 5 else d207 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())}))208 for comp in ["y4_f5", "y5_f5", "y6_f5"]:209 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())}))210 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")]:211 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())}))212 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())}))213 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())}))214 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())}))215 # Oster / Cinelli-Hazlett style sensitivity on the LPM version of the primary216 try:217 dd = d.dropna(subset=["y456_f5", "P_modern_z", "log_gdppc", "v2x_polyarchy", "t_since", "e_regionpol_6C"]).copy()218 # distinct event-years (not windows) in the estimation sample, for the text219 (RES / "h2_modern_event_counts.json").write_text(json.dumps({"country_years": int(len(dd)), "event_windows": int(dd.y456_f5.sum()), "distinct_event_years_in_sample": int((dd.any_event == 1).sum()), "distinct_event_years_panel": int((d.any_event == 1).sum())}))220 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()221 restricted = smf.ols(f"y456_f5 ~ {XS} + {FE}", dd).fit()222 r2_full, r2_short = full.rsquared, short.rsquared; b_full, b_short = full.params["P_modern_z"], short.params["P_modern_z"]223 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)224 t = full.tvalues["P_modern_z"]; dfres = full.df_resid; r2_yd = t ** 2 / (t ** 2 + dfres)225 rv = 0.5 * (np.sqrt(r2_yd ** 2 + 4 * r2_yd) - r2_yd)226 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)})227 except Exception as e:228 rows.append({"model": "Sensitivity", "error": str(e)[:100]})229 out = pd.DataFrame(rows); out.to_csv(RES / "h2_modern_robustness.csv", index=False); return out230231232def main():233 d = build_sample(5); d.to_parquet(PROC / "modern_h2_sample.parquet", index=False)234 prim = primary(d); prim.to_csv(RES / "h2_modern_primary.csv", index=False)235 wins, summ = oos(d)236 h4s = h4(d); (RES / "h4_modern_summary.json").write_text(json.dumps(h4s, indent=2))237 h5s = h5(d); rob = robustness(d)238 pd.set_option("display.width", 250); pd.set_option("display.max_colwidth", 70)239 cols = ["model", "term", "coef", "se", "p", "n", "events", "error"]240 print("=== PRIMARY ===\n", prim[[c for c in cols if c in prim.columns]].round(3).to_string(index=False))241 print("\n=== OOS wins ===\n", wins.round(3).to_string(index=False)); print("\n=== OOS summary ===\n", summ.round(3).to_string(index=False))242 print("\n=== H4 ===", json.dumps(h4s, indent=1))243 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))244 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))245246247if __name__ == "__main__":248 main()249