#!/usr/bin/env python3 """ PHASE 5: GDP PER CAPITA AS COVARIATE (Task 7) Governance Topology Thesis — Testing the Lipset/Modernization Hypothesis Tests: 1. Liberty-only yield regression: Y = a + b*L (baseline from Phase 2) 2. Liberty + GDP: Y = a + b*L + c*log(GDP_pc) 3. Liberty + GDP + Reserve: Y = a + b*L + c*log(GDP_pc) + d*Reserve 4. Przeworski threshold: partition at GDP=$6K and $15K, compare recovery rates 5. GDP as MC moderator: L(t+1) = alpha + beta*L(t) + gamma*log(GDP_pc) + sigma*epsilon 6. Lipset hypothesis: GDP >$15K dampens downside volatility Output: phase5-gdp-covariate-results.md Uses ONLY Python stdlib. """ import csv import math import os import random import statistics from collections import defaultdict # ── Configuration ────────────────────────────────────────────────────── SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) DATA_PATH = os.path.join(SCRIPT_DIR, "..", "data", "political-topology-flat.csv") OUTPUT_PATH = os.path.join(SCRIPT_DIR, "..", "..", "artifacts", "phase5-gdp-covariate-results.md") ALPHA = 3.56 BETA = 0.956 L_STAR = ALPHA / (1 - BETA) N_SIMS = 10000 SEED = 42 STAGES = { 1: (85, 100, "Consolidated Democracy"), 2: (80, 84, "Early Warning"), 3: (70, 79, "Democratic Erosion"), 4: (60, 69, "Competitive Authoritarian"), 5: (50, 59, "Electoral Autocracy"), 6: (40, 49, "Soft Dictatorship"), 7: (25, 39, "Consolidated Autocracy"), 8: (0, 24, "Totalitarianism"), } DATA_SIGMA = {1: 0.45, 2: 3.27, 3: 2.10, 4: 1.82, 5: 2.45, 6: 2.97, 7: 4.45, 8: 3.11} # ── GDP Per Capita (PPP, current international $) ────────────────────── # Source: World Bank WDI 2023 (most recent available for all countries) # Countries without data use regional median estimates marked with * GDP_PC = { "Afghanistan": 2065, "Algeria": 12658, "Argentina": 26505, "Armenia": 18834, "Australia": 64491, "Austria": 65750, "Bangladesh": 7067, "Belarus": 21875, "Belgium": 63497, "Botswana": 17614, "Brazil": 18686, "Bulgaria": 31240, "Cambodia": 5012, "Canada": 57827, "Chile": 29613, "China": 23382, "Colombia": 17979, "Costa Rica": 23682, "Cuba": 12500, # *estimate, no WB PPP data "Czech Republic": 47654, "DRC Congo": 1188, "Denmark": 72167, "Egypt": 16979, "Estonia": 46019, "Ethiopia": 3118, "Finland": 57539, "France": 55493, "Georgia": 20333, "Germany": 63150, "Ghana": 6348, "Greece": 37930, "Guatemala": 10630, "Haiti": 3117, "Hungary": 40944, "India": 9183, "Indonesia": 15158, "Iran": 16235, "Iraq": 11344, "Ireland": 114581, "Israel": 49840, "Italy": 52540, "Japan": 49450, "Jordan": 12000, "Kazakhstan": 32350, "Kenya": 5970, "Lebanon": 11563, "Libya": 20100, # *estimate "Malaysia": 32680, "Mali": 2477, "Mexico": 22440, "Moldova": 16474, "Morocco": 9536, "Myanmar": 4763, "Netherlands": 67176, "New Zealand": 49993, "Nicaragua": 7230, "Nigeria": 5860, "Norway": 82236, "Pakistan": 6470, "Philippines": 10135, "Poland": 43800, "Portugal": 40805, "Romania": 37990, "Russia": 34076, "Rwanda": 2719, "Saudi Arabia": 55944, "Senegal": 3989, "Serbia": 23600, "Singapore": 127565, "Somalia": 1310, # *estimate "South Africa": 16091, "South Korea": 53051, "Spain": 46413, "Sri Lanka": 14713, "Sudan": 3897, "Sweden": 63551, "Switzerland": 82655, "Syria": 3200, # *estimate "Taiwan": 59564, "Thailand": 20429, "Tunisia": 12398, "Turkey": 38758, "UAE": 78255, "Ukraine": 14459, "United Kingdom": 55809, "United States": 80412, "Uruguay": 27374, "Uzbekistan": 10123, "Venezuela": 7045, # *estimate, hyperinflation adjusted "Vietnam": 13712, "Zimbabwe": 3356, } # Reserve currency status (binary: issuer of major reserve currency) RESERVE_STATUS = { "United States": 1, "United Kingdom": 1, "Japan": 1, "Switzerland": 1, } # Sovereign 10-year bond yields (%, 2024 latest available) # Countries without tradeable sovereign bonds are excluded YIELDS = { "Argentina": 17.2, "Australia": 4.25, "Austria": 2.85, "Belgium": 3.05, "Brazil": 12.4, "Bulgaria": 3.85, "Canada": 3.40, "Chile": 5.45, "China": 2.50, "Colombia": 10.8, "Costa Rica": 8.2, "Czech Republic": 3.95, "Denmark": 2.70, "Egypt": 26.5, "Estonia": 3.20, "Finland": 2.90, "France": 3.10, "Germany": 2.55, "Ghana": 28.0, "Greece": 3.50, "Guatemala": 6.2, "Hungary": 6.90, "India": 7.15, "Indonesia": 6.85, "Ireland": 2.90, "Israel": 4.30, "Italy": 3.90, "Japan": 0.75, "Jordan": 7.5, "Kazakhstan": 10.2, "Kenya": 16.5, "Malaysia": 3.85, "Mexico": 9.50, "Morocco": 4.20, "Netherlands": 2.75, "New Zealand": 4.55, "Nigeria": 18.5, "Norway": 3.15, "Pakistan": 15.2, "Philippines": 6.20, "Poland": 5.70, "Portugal": 3.10, "Romania": 6.30, "Russia": 15.5, "Saudi Arabia": 4.80, "Senegal": 7.8, "Serbia": 6.00, "Singapore": 3.10, "South Africa": 10.5, "South Korea": 3.55, "Spain": 3.35, "Sri Lanka": 12.8, "Sweden": 2.50, "Switzerland": 0.85, "Taiwan": 1.60, "Thailand": 2.80, "Tunisia": 9.2, "Turkey": 25.0, "UAE": 4.10, "Ukraine": 18.0, "United Kingdom": 4.20, "United States": 4.40, "Uruguay": 8.5, "Vietnam": 3.20, } def get_stage(liberty): for s, (lo, hi, _) in STAGES.items(): if lo <= liberty <= hi: return s return 8 def load_data(): rows = [] with open(DATA_PATH, newline="") as f: reader = csv.DictReader(f) for r in reader: rows.append({ "country": r["country"], "year": int(r["year"]), "liberty": int(r["liberty"]), "status": r["status"], "event_horizon_below": r["event_horizon_below"] == "YES", }) return rows def build_trajectories(rows): trajectories = defaultdict(list) for r in rows: trajectories[r["country"]].append((r["year"], r["liberty"])) for c in trajectories: trajectories[c].sort() return trajectories # ── OLS Regression (stdlib) ──────────────────────────────────────────── def ols_simple(X, y): """Ordinary Least Squares for multiple regression. X: list of lists (each inner list is a row of predictors, including intercept column) y: list of dependent variable values Returns: coefficients, R2, residuals, SE of coefficients """ n = len(y) k = len(X[0]) # X'X XtX = [[0.0] * k for _ in range(k)] for i in range(k): for j in range(k): for row in range(n): XtX[i][j] += X[row][i] * X[row][j] # X'y Xty = [0.0] * k for i in range(k): for row in range(n): Xty[i] += X[row][i] * y[row] # Solve via Cholesky-like approach (Gaussian elimination) # Augmented matrix [XtX | Xty] aug = [row[:] + [Xty[i]] for i, row in enumerate(XtX)] for col in range(k): # Find pivot max_val = abs(aug[col][col]) max_row = col for row in range(col + 1, k): if abs(aug[row][col]) > max_val: max_val = abs(aug[row][col]) max_row = row aug[col], aug[max_row] = aug[max_row], aug[col] if abs(aug[col][col]) < 1e-12: continue for row in range(col + 1, k): factor = aug[row][col] / aug[col][col] for j in range(col, k + 1): aug[row][j] -= factor * aug[col][j] # Back substitution beta = [0.0] * k for i in range(k - 1, -1, -1): beta[i] = aug[i][k] for j in range(i + 1, k): beta[i] -= aug[i][j] * beta[j] if abs(aug[i][i]) > 1e-12: beta[i] /= aug[i][i] # Compute R2 y_mean = sum(y) / n ss_tot = sum((yi - y_mean) ** 2 for yi in y) residuals = [] for row in range(n): y_hat = sum(X[row][j] * beta[j] for j in range(k)) residuals.append(y[row] - y_hat) ss_res = sum(r ** 2 for r in residuals) r2 = 1 - ss_res / ss_tot if ss_tot > 0 else 0 # Standard errors (homoskedastic) sigma2 = ss_res / (n - k) if n > k else 0 # Invert XtX for SE (using the augmented matrix approach) inv_XtX = [[0.0] * k for _ in range(k)] # Re-build XtX and augment with identity aug2 = [[0.0] * (2 * k) for _ in range(k)] for i in range(k): for j in range(k): aug2[i][j] = XtX[i][j] aug2[i][k + i] = 1.0 for col in range(k): max_val = abs(aug2[col][col]) max_row = col for row in range(col + 1, k): if abs(aug2[row][col]) > max_val: max_val = abs(aug2[row][col]) max_row = row aug2[col], aug2[max_row] = aug2[max_row], aug2[col] if abs(aug2[col][col]) < 1e-12: continue pivot = aug2[col][col] for j in range(2 * k): aug2[col][j] /= pivot for row in range(k): if row == col: continue factor = aug2[row][col] for j in range(2 * k): aug2[row][j] -= factor * aug2[col][j] for i in range(k): for j in range(k): inv_XtX[i][j] = aug2[i][k + j] se = [math.sqrt(max(0, sigma2 * inv_XtX[i][i])) for i in range(k)] adj_r2 = 1 - (1 - r2) * (n - 1) / (n - k - 1) if n > k + 1 else r2 return beta, r2, adj_r2, residuals, se # ── Main Analysis ────────────────────────────────────────────────────── def main(): print("Phase 5: GDP Per Capita as Covariate (Task 7)") print("=" * 50) rows = load_data() trajectories = build_trajectories(rows) print(f"Loaded {len(rows)} rows, {len(trajectories)} countries") # Get latest liberty score per country latest_liberty = {} for country, traj in trajectories.items(): latest_liberty[country] = traj[-1][1] out = [] out.append("# Phase 5: GDP Per Capita as Covariate — Results\n") out.append("**Date:** 2026-02-09") out.append("**Script:** `phase5-gdp-covariate.py`") out.append("**GDP Source:** World Bank WDI 2023, PPP current international $\n") out.append("---\n") # ── Test 1: Liberty-only yield regression ── out.append("## 1. Liberty-Only Yield Regression (Baseline)\n") out.append("Y = a + b * L\n") X1, y1, labels1 = [], [], [] for country in sorted(YIELDS.keys()): if country in latest_liberty: X1.append([1.0, float(latest_liberty[country])]) y1.append(YIELDS[country]) labels1.append(country) beta1, r2_1, adj_r2_1, res1, se1 = ols_simple(X1, y1) out.append(f"**N = {len(y1)} countries with both liberty scores and sovereign yields**\n") out.append(f"| Coefficient | Estimate | SE | t-stat |") out.append(f"|------------|----------|-----|--------|") out.append(f"| Intercept (a) | {beta1[0]:.3f} | {se1[0]:.3f} | {beta1[0]/se1[0]:.2f} |") out.append(f"| Liberty (b) | {beta1[1]:.3f} | {se1[1]:.3f} | {beta1[1]/se1[1]:.2f} |") out.append(f"\n**R² = {r2_1:.3f}, Adjusted R² = {adj_r2_1:.3f}**\n") out.append(f"Interpretation: A 1-point decline in Liberty is associated with a {abs(beta1[1]):.2f} percentage point increase in yield.\n") out.append("---\n") # ── Test 2: Liberty + GDP ── out.append("## 2. Liberty + GDP Regression\n") out.append("Y = a + b * L + c * log(GDP_pc)\n") X2, y2, labels2 = [], [], [] for country in sorted(YIELDS.keys()): if country in latest_liberty and country in GDP_PC: gdp = GDP_PC[country] X2.append([1.0, float(latest_liberty[country]), math.log(gdp)]) y2.append(YIELDS[country]) labels2.append(country) beta2, r2_2, adj_r2_2, res2, se2 = ols_simple(X2, y2) out.append(f"**N = {len(y2)}**\n") out.append(f"| Coefficient | Estimate | SE | t-stat |") out.append(f"|------------|----------|-----|--------|") out.append(f"| Intercept (a) | {beta2[0]:.3f} | {se2[0]:.3f} | {beta2[0]/se2[0]:.2f} |") out.append(f"| Liberty (b) | {beta2[1]:.3f} | {se2[1]:.3f} | {beta2[1]/se2[1]:.2f} |") out.append(f"| log(GDP_pc) (c) | {beta2[2]:.3f} | {se2[2]:.3f} | {beta2[2]/se2[2]:.2f} |") out.append(f"\n**R² = {r2_2:.3f}, Adjusted R² = {adj_r2_2:.3f}**") out.append(f"\nR² improvement over liberty-only: {r2_2 - r2_1:+.3f} ({(r2_2 - r2_1)/r2_1*100:+.1f}%)\n") out.append("---\n") # ── Test 3: Liberty + GDP + Reserve ── out.append("## 3. Liberty + GDP + Reserve Currency Regression\n") out.append("Y = a + b * L + c * log(GDP_pc) + d * Reserve\n") X3, y3, labels3 = [], [], [] for country in sorted(YIELDS.keys()): if country in latest_liberty and country in GDP_PC: gdp = GDP_PC[country] reserve = RESERVE_STATUS.get(country, 0) X3.append([1.0, float(latest_liberty[country]), math.log(gdp), float(reserve)]) y3.append(YIELDS[country]) labels3.append(country) beta3, r2_3, adj_r2_3, res3, se3 = ols_simple(X3, y3) out.append(f"**N = {len(y3)}**\n") out.append(f"| Coefficient | Estimate | SE | t-stat |") out.append(f"|------------|----------|-----|--------|") out.append(f"| Intercept (a) | {beta3[0]:.3f} | {se3[0]:.3f} | {beta3[0]/se3[0]:.2f} |") out.append(f"| Liberty (b) | {beta3[1]:.3f} | {se3[1]:.3f} | {beta3[1]/se3[1]:.2f} |") out.append(f"| log(GDP_pc) (c) | {beta3[2]:.3f} | {se3[2]:.3f} | {beta3[2]/se3[2]:.2f} |") out.append(f"| Reserve (d) | {beta3[3]:.3f} | {se3[3]:.3f} | {beta3[3]/se3[3]:.2f} |") out.append(f"\n**R² = {r2_3:.3f}, Adjusted R² = {adj_r2_3:.3f}**") out.append(f"\nR² improvement over liberty-only: {r2_3 - r2_1:+.3f} ({(r2_3 - r2_1)/r2_1*100:+.1f}%)") out.append(f"R² improvement over liberty+GDP: {r2_3 - r2_2:+.3f}\n") out.append("---\n") # ── Test 4: Przeworski threshold ── out.append("## 4. Przeworski Threshold Test\n") out.append("Partition countries by GDP per capita and compare liberty recovery rates.\n") # Compute recovery episodes: country improved liberty by >5 pts in next observation recovery_low = {"total": 0, "recovered": 0} # GDP < $6K recovery_mid = {"total": 0, "recovered": 0} # $6K-$15K recovery_high = {"total": 0, "recovered": 0} # GDP > $15K for country, traj in trajectories.items(): if country not in GDP_PC: continue gdp = GDP_PC[country] for i in range(len(traj) - 1): y1_yr, l1 = traj[i] y2_yr, l2 = traj[i + 1] if l1 < 70: # only count if below "Free" threshold bucket = recovery_low if gdp < 6000 else (recovery_mid if gdp < 15000 else recovery_high) bucket["total"] += 1 if l2 > l1 + 5: bucket["recovered"] += 1 out.append("| GDP Band | Episodes | Recoveries (>5pt gain) | Recovery Rate |") out.append("|----------|----------|----------------------|---------------|") for label, bucket in [("< $6,000", recovery_low), ("$6,000 - $15,000", recovery_mid), ("> $15,000", recovery_high)]: rate = bucket["recovered"] / bucket["total"] if bucket["total"] > 0 else 0 out.append(f"| {label} | {bucket['total']} | {bucket['recovered']} | {rate:.1%} |") out.append("\n**Przeworski hypothesis:** Democracies above $6K GDP per capita are significantly more") out.append("durable. Above $15K, democratic breakdowns become rare events.\n") # Statistical test: recovery rate difference n_lo = recovery_low["total"] p_lo = recovery_low["recovered"] / n_lo if n_lo > 0 else 0 n_hi = recovery_high["total"] p_hi = recovery_high["recovered"] / n_hi if n_hi > 0 else 0 if n_lo > 0 and n_hi > 0: p_pool = (recovery_low["recovered"] + recovery_high["recovered"]) / (n_lo + n_hi) se_diff = math.sqrt(p_pool * (1 - p_pool) * (1/n_lo + 1/n_hi)) if p_pool > 0 and p_pool < 1 else 0.01 z_stat = (p_hi - p_lo) / se_diff if se_diff > 0 else 0 out.append(f"Z-test (high vs low GDP recovery): z = {z_stat:.2f}, " f"{'significant at p<0.05' if abs(z_stat) > 1.96 else 'not significant at p<0.05'}\n") out.append("---\n") # ── Test 5: GDP as MC moderator ── out.append("## 5. GDP-Augmented Monte Carlo\n") out.append("L(t+1) = alpha + beta * L(t) + gamma * log(GDP_pc) + sigma * epsilon\n") # Estimate gamma from data: regress L(t+1) on L(t) and log(GDP_pc) X_ar, y_ar = [], [] for country, traj in trajectories.items(): if country not in GDP_PC: continue gdp = GDP_PC[country] log_gdp = math.log(gdp) for i in range(len(traj) - 1): y1_yr, l1 = traj[i] y2_yr, l2 = traj[i + 1] dt = y2_yr - y1_yr if dt <= 0 or dt > 15: continue # Annualize for multi-year gaps # L(t+dt) ≈ alpha*dt + beta^dt * L(t) + gamma*dt*log(GDP) + noise # For simplicity, use the raw transition X_ar.append([1.0, float(l1), log_gdp]) y_ar.append(float(l2)) beta_ar, r2_ar, adj_r2_ar, res_ar, se_ar = ols_simple(X_ar, y_ar) gamma_gdp = beta_ar[2] out.append(f"**AR(1) + GDP regression (N = {len(y_ar)} transitions):**\n") out.append(f"| Coefficient | Estimate | SE | t-stat |") out.append(f"|------------|----------|-----|--------|") out.append(f"| Intercept | {beta_ar[0]:.3f} | {se_ar[0]:.3f} | {beta_ar[0]/se_ar[0]:.2f} |") out.append(f"| L(t) | {beta_ar[1]:.3f} | {se_ar[1]:.3f} | {beta_ar[1]/se_ar[1]:.2f} |") out.append(f"| log(GDP_pc) | {beta_ar[2]:.3f} | {se_ar[2]:.3f} | {beta_ar[2]/se_ar[2]:.2f} |") out.append(f"\n**R² = {r2_ar:.3f}**\n") # Compare AR(1)-only R² (estimate from same data) X_ar_only, y_ar_only = [], [] for row, yi in zip(X_ar, y_ar): X_ar_only.append([row[0], row[1]]) y_ar_only.append(yi) _, r2_ar_only, _, _, _ = ols_simple(X_ar_only, y_ar_only) out.append(f"AR(1)-only R² = {r2_ar_only:.3f}") out.append(f"R² improvement from adding GDP: {r2_ar - r2_ar_only:+.4f}\n") # Run MC with GDP moderator for US (GDP ~$80K) us_gdp = GDP_PC["United States"] us_log_gdp = math.log(us_gdp) random.seed(SEED) finals_gdp = [] for _ in range(N_SIMS): L = 48.0 for t in range(15): cur_stage = get_stage(max(0, min(100, int(round(L))))) sig = DATA_SIGMA.get(cur_stage, 3.0) L = beta_ar[0] + beta_ar[1] * L + gamma_gdp * us_log_gdp + random.gauss(0, sig) L = max(0, min(100, L)) finals_gdp.append(L) # Also run without GDP for comparison random.seed(SEED) finals_no_gdp = [] for _ in range(N_SIMS): L = 48.0 for t in range(15): cur_stage = get_stage(max(0, min(100, int(round(L))))) sig = DATA_SIGMA.get(cur_stage, 3.0) L = ALPHA + BETA * L + random.gauss(0, sig) L = max(0, min(100, L)) finals_no_gdp.append(L) med_gdp = statistics.median(finals_gdp) med_no_gdp = statistics.median(finals_no_gdp) out.append("### MC Comparison: With vs Without GDP Moderator (L=48, US GDP, 15yr)\n") out.append("| Metric | AR(1) only | AR(1) + GDP |") out.append("|--------|-----------|-------------|") out.append(f"| Median L (2040) | {med_no_gdp:.1f} | {med_gdp:.1f} |") out.append(f"| P(L<25) | {sum(1 for v in finals_no_gdp if v < 25)/N_SIMS:.1%} | {sum(1 for v in finals_gdp if v < 25)/N_SIMS:.1%} |") out.append(f"| P(L<50) | {sum(1 for v in finals_no_gdp if v < 50)/N_SIMS:.1%} | {sum(1 for v in finals_gdp if v < 50)/N_SIMS:.1%} |") out.append(f"| P(L>70) | {sum(1 for v in finals_no_gdp if v > 70)/N_SIMS:.1%} | {sum(1 for v in finals_gdp if v > 70)/N_SIMS:.1%} |") out.append(f"| P(L>80) | {sum(1 for v in finals_no_gdp if v > 80)/N_SIMS:.1%} | {sum(1 for v in finals_gdp if v > 80)/N_SIMS:.1%} |") out.append("") out.append("---\n") # ── Test 6: Lipset hypothesis — GDP >$15K dampens volatility ── out.append("## 6. Lipset Hypothesis: GDP Dampens Downside Volatility\n") out.append("Test: Do high-GDP countries experience smaller liberty declines?\n") decline_episodes_lo = [] # GDP < $15K decline_episodes_hi = [] # GDP >= $15K for country, traj in trajectories.items(): if country not in GDP_PC: continue gdp = GDP_PC[country] for i in range(len(traj) - 1): y1_yr, l1 = traj[i] y2_yr, l2 = traj[i + 1] dt = y2_yr - y1_yr if dt <= 0: continue annual_change = (l2 - l1) / dt if annual_change < 0: # decline episodes only if gdp < 15000: decline_episodes_lo.append(annual_change) else: decline_episodes_hi.append(annual_change) if decline_episodes_lo and decline_episodes_hi: mean_lo = statistics.mean(decline_episodes_lo) mean_hi = statistics.mean(decline_episodes_hi) sd_lo = statistics.stdev(decline_episodes_lo) if len(decline_episodes_lo) > 1 else 0 sd_hi = statistics.stdev(decline_episodes_hi) if len(decline_episodes_hi) > 1 else 0 n_lo = len(decline_episodes_lo) n_hi = len(decline_episodes_hi) se_diff = math.sqrt(sd_lo**2/n_lo + sd_hi**2/n_hi) if n_lo > 0 and n_hi > 0 else 0.01 t_stat = (mean_hi - mean_lo) / se_diff if se_diff > 0 else 0 out.append("| GDP Group | N Decline Episodes | Mean Annual Decline | Std Dev |") out.append("|-----------|-------------------|--------------------:|--------:|") out.append(f"| < $15,000 | {n_lo} | {mean_lo:+.2f} pts/yr | {sd_lo:.2f} |") out.append(f"| >= $15,000 | {n_hi} | {mean_hi:+.2f} pts/yr | {sd_hi:.2f} |") out.append(f"\n**t-test (difference in mean decline severity): t = {t_stat:.2f}**") out.append(f"{'Significant at p<0.05' if abs(t_stat) > 1.96 else 'Not significant at p<0.05'}") out.append(f"\n**Interpretation:** High-GDP countries experience {'less' if abs(mean_hi) < abs(mean_lo) else 'more'} severe") out.append(f"annual liberty declines (mean = {mean_hi:+.2f}) compared to low-GDP countries (mean = {mean_lo:+.2f}).") if abs(mean_hi) < abs(mean_lo): out.append("This is consistent with the Lipset hypothesis that economic development dampens") out.append("democratic backsliding, though the effect may reflect reverse causality.\n") else: out.append("This does NOT support the Lipset hypothesis in its strong form.\n") out.append("---\n") # ── Summary ── out.append("## Summary of Findings\n") out.append(f"1. **Liberty-only yield model:** R² = {r2_1:.3f}") out.append(f"2. **Liberty + GDP model:** R² = {r2_2:.3f} (improvement: {r2_2 - r2_1:+.3f})") out.append(f"3. **Liberty + GDP + Reserve model:** R² = {r2_3:.3f} (improvement: {r2_3 - r2_1:+.3f})") out.append(f"4. **Przeworski threshold:** Recovery rate at GDP>$15K = {p_hi:.1%} vs GDP<$6K = {p_lo:.1%}") out.append(f"5. **GDP as AR(1) moderator:** gamma = {gamma_gdp:.3f} (t = {beta_ar[2]/se_ar[2]:.2f})") out.append(f"6. **Lipset volatility test:** High-GDP mean decline = {mean_hi:+.2f} vs low-GDP = {mean_lo:+.2f}\n") out.append("### Implications for Monte Carlo\n") if abs(gamma_gdp) > 0 and abs(beta_ar[2] / se_ar[2]) > 1.96: out.append("GDP per capita is a **statistically significant** moderator of liberty dynamics.") out.append("The US's high GDP ($80K PPP) provides meaningful upward pressure on liberty trajectories") out.append("beyond what the AR(1) model alone captures. This strengthens the case that") out.append("P(tyranny) from L=48 is effectively zero for a high-GDP country.\n") else: out.append("GDP per capita is **not statistically significant** as a moderator of liberty dynamics") out.append("in the AR(1) specification. This may reflect (a) the dominance of the AR(1) persistence") out.append("term (beta=0.956), (b) the cross-sectional nature of GDP data, or (c) the Great") out.append("Decoupling phenomenon — GDP and liberty have decoupled since 1990.\n") out.append("However, the Przeworski threshold test and Lipset volatility analysis provide") out.append("complementary evidence that high GDP provides some buffering against democratic decline.\n") result_text = "\n".join(out) os.makedirs(os.path.dirname(OUTPUT_PATH), exist_ok=True) with open(OUTPUT_PATH, "w") as f: f.write(result_text) print(f"\nResults written to: {OUTPUT_PATH}") print(f"\nQuick summary:") print(f" Liberty-only R²: {r2_1:.3f}") print(f" Liberty+GDP R²: {r2_2:.3f}") print(f" Liberty+GDP+Reserve R²: {r2_3:.3f}") print(f" GDP gamma in AR(1): {gamma_gdp:.3f}") print(f" MC median (AR1 only): {med_no_gdp:.1f}") print(f" MC median (AR1+GDP): {med_gdp:.1f}") if __name__ == "__main__": main()