import pandas as pd, numpy as np
import statsmodels.formula.api as smf
from linearmodels.panel import PanelOLS
p=pd.read_csv('/tmp/rol_panel.csv')
out=[]
def w(s): out.append(s); print(s)

w("="*70)
w("FUKUYAMA THREE-PILLARS TEST: does the rule of law need the ballot?")
w("DV = WGI Rule of Law (perception index, ~ -2.5..2.5). N=%d, %d countries, %d-%d"%(
    len(p),p.code.nunique(),p.year.min(),p.year.max()))
w("="*70)

# Model 1: polyarchy alone (+year FE)
m1=smf.ols('rle ~ polyarchy + C(year)',data=p).fit(cov_type='HC1')
w("\n[M1] RoL ~ electoral democracy (polyarchy) alone, year FE")
w("   polyarchy b=%.3f (se %.3f, p=%.1e)  R2=%.3f"%(m1.params['polyarchy'],m1.bse['polyarchy'],m1.pvalues['polyarchy'],m1.rsquared))

# Model 2: add income
m2=smf.ols('rle ~ polyarchy + log_gdppc + C(year)',data=p).fit(cov_type='HC1')
w("\n[M2] RoL ~ polyarchy + log GDP pc, year FE  (HC1)")
for v in ['polyarchy','log_gdppc']:
    w("   %-10s b=%.3f (se %.3f, p=%.1e)"%(v,m2.params[v],m2.bse[v],m2.pvalues[v]))
w("   R2=%.3f"%m2.rsquared)

# Standardized betas (compare magnitudes)
m2z=smf.ols('rle_z ~ polyarchy_z + log_gdppc_z + C(year)',data=p).fit(cov_type='HC1')
w("\n[M2-standardized] beta polyarchy=%.3f  beta log_gdppc=%.3f  -> income carries %.1fx the ballot"%(
    m2z.params['polyarchy_z'],m2z.params['log_gdppc_z'],m2z.params['log_gdppc_z']/m2z.params['polyarchy_z']))

# Model 3: add fiscal capacity (tax/GDP) where available
p3=p.dropna(subset=['tax'])
m3=smf.ols('rle ~ polyarchy + log_gdppc + tax + C(year)',data=p3).fit(cov_type='HC1')
w("\n[M3] + fiscal capacity (tax %%GDP), N=%d"%len(p3))
for v in ['polyarchy','log_gdppc','tax']:
    w("   %-10s b=%.3f (se %.3f, p=%.1e)"%(v,m3.params[v],m3.bse[v],m3.pvalues[v]))

# Model 4: TWO-WAY FE (country + year): within-country
pf=p.set_index(['code','year'])
m4=PanelOLS.from_formula('rle ~ polyarchy + log_gdppc + EntityEffects + TimeEffects',data=pf).fit(cov_type='clustered',cluster_entity=True)
w("\n[M4] TWO-WAY FE (country+year), cluster by country  -- WITHIN-country")
for v in ['polyarchy','log_gdppc']:
    w("   %-10s b=%.3f (se %.3f, p=%.3f)"%(v,m4.params[v],m4.std_errors[v],m4.pvalues[v]))
w("   within R2=%.3f"%m4.rsquared_within)

# Robustness DV: Gov Effectiveness (note shared-source caveat)
m5=smf.ols('gee ~ polyarchy + log_gdppc + C(year)',data=p).fit(cov_type='HC1')
w("\n[M5 robustness] DV=Gov Effectiveness: polyarchy b=%.3f (p=%.1e), log_gdppc b=%.3f (p=%.1e)"%(
    m5.params['polyarchy'],m5.pvalues['polyarchy'],m5.params['log_gdppc'],m5.pvalues['log_gdppc']))

# ---- The affirmative cases: high rule of law, low electoral democracy ----
latest=p[p.year>=2010].groupby('code').agg(rle=('rle','mean'),polyarchy=('polyarchy','mean'),gdppc=('gdppc','mean')).dropna()
rol_hi=latest.rle.quantile(0.75); poly_lo=latest.polyarchy.median()
cases=latest[(latest.rle>=rol_hi)&(latest.polyarchy<=poly_lo)].sort_values('rle',ascending=False)
w("\n"+"="*70)
w("HIGH rule of law (top quartile, RoL>=%.2f) BUT LOW electoral democracy (<= median polyarchy %.2f), 2010-12 avg:"%(rol_hi,poly_lo))
import pycountry
def nm(c):
    try: return pycountry.countries.get(alpha_3=c).name
    except: return c
try:
    for c,row in cases.iterrows():
        w("   %-5s RoL=%+.2f  polyarchy=%.2f"%(c,row.rle,row.polyarchy))
except Exception as e:
    for c,row in cases.iterrows():
        w("   %-5s RoL=%+.2f  polyarchy=%.2f"%(c,row.rle,row.polyarchy))
w("   COUNT of such 'liberal-institutions-without-the-ballot' cases: %d"%len(cases))

# share of variance: partial R2 of polyarchy beyond income
base=smf.ols('rle ~ log_gdppc + C(year)',data=p).fit()
full=smf.ols('rle ~ log_gdppc + polyarchy + C(year)',data=p).fit()
w("\nIncome+year alone R2=%.3f ; adding the ballot raises it to R2=%.3f (gain %.3f)"%(
    base.rsquared,full.rsquared,full.rsquared-base.rsquared))
base2=smf.ols('rle ~ polyarchy + C(year)',data=p).fit()
w("Ballot+year alone R2=%.3f ; adding income raises it to R2=%.3f (gain %.3f)"%(
    base2.rsquared,full.rsquared,full.rsquared-base2.rsquared))

open('/tmp/results_rol.txt','w').write("\n".join(out))
print("\nSAVED /tmp/results_rol.txt")
