import pandas as pd, numpy as np
import statsmodels.formula.api as smf
from linearmodels.panel import PanelOLS
p=pd.read_csv('/sessions/serene-beautiful-maxwell/mnt/Shackles of Democracy/Regression - Evidence of Decline/glob_panel.csv')
p=p.dropna(subset=['ge','trade','log_gdppc']).copy()
def line(s): print(s)
line("="*70)
line("MODEL 21: does the connected world push governments toward CAPABLE,")
line("ECONOMICALLY EFFECTIVE government?  DV = WGI Government Effectiveness")
line("N=%d, %d countries, %d-%d"%(len(p),p.code.nunique(),p.year.min(),p.year.max()))
line("corr(gov_effectiveness, trade openness)=%.3f ; corr(ge, log GDPpc)=%.3f"%(p['ge'].corr(p.trade),p['ge'].corr(p.log_gdppc)))

# pooled OLS with year FE, cluster by country
m1=smf.ols('ge ~ trade + log_gdppc + C(year)',data=p).fit(cov_type='cluster',cov_kwds={'groups':p['code']})
line("\n[A] ge ~ trade + log GDPpc, year FE, cluster by country (across countries)")
for v in ['trade','log_gdppc']: line("   %-10s b=%.4f  p=%.3g"%(v,m1.params[v],m1.pvalues[v]))

# two-way FE (within country)
pf=p.set_index(['code','year'])
m2=PanelOLS.from_formula('ge ~ trade + log_gdppc + EntityEffects + TimeEffects',data=pf).fit(cov_type='clustered',cluster_entity=True)
line("\n[B] TWO-WAY FE (country+year): within a country, does rising openness track rising effectiveness?")
for v in ['trade','log_gdppc']: line("   %-10s b=%.4f  p=%.3g"%(v,m2.params[v],m2.pvalues[v]))

# broadband
pb=p.dropna(subset=['bband'])
m3=smf.ols('ge ~ bband + log_gdppc + C(year)',data=pb).fit(cov_type='cluster',cov_kwds={'groups':pb['code']})
line("\n[C] ge ~ broadband + log GDPpc, year FE (information connection), N=%d"%len(pb))
for v in ['bband','log_gdppc']: line("   %-10s b=%.4f  p=%.3g"%(v,m3.params[v],m3.pvalues[v]))

# NON-DEMOCRACY subsample (poly<0.5)
nd=p[p.poly<0.5].copy()
m4=smf.ols('ge ~ trade + log_gdppc + C(year)',data=nd).fit(cov_type='cluster',cov_kwds={'groups':nd['code']})
line("\n[D] NON-DEMOCRACIES ONLY (polyarchy<0.5): N=%d, %d countries"%(len(nd),nd.code.nunique()))
for v in ['trade','log_gdppc']: line("   %-10s b=%.4f  p=%.3g"%(v,m4.params[v],m4.pvalues[v]))
ndb=nd.dropna(subset=['bband'])
m5=smf.ols('ge ~ bband + log_gdppc + C(year)',data=ndb).fit(cov_type='cluster',cov_kwds={'groups':ndb['code']})
line("   [non-dem broadband] bband b=%.4f p=%.3g (N=%d)"%(m5.params['bband'],m5.pvalues['bband'],len(ndb)))
