import pandas as pd, numpy as np
import statsmodels.formula.api as smf
from linearmodels.panel import PanelOLS
import pyreadr

def parse_wgi(path,name):
    raw=pd.read_csv(path,header=None,dtype=str,encoding='latin-1')
    years=raw.iloc[0].tolist(); sub=raw.iloc[1].tolist(); recs=[]
    for r in range(2,len(raw)):
        row=raw.iloc[r]; code=row[1]
        if not isinstance(code,str) or len(code)!=3: continue
        for c in range(2,len(row)):
            if sub[c]=='Estimate':
                try: y=int(float(years[c])); v=float(row[c])
                except: continue
                recs.append((code.strip().upper(),y,v))
    return pd.DataFrame(recs,columns=['code','year',name])

corr=parse_wgi('/tmp/wgirepo/Corruption.csv','cc')      # control of corruption
gov =parse_wgi('/tmp/wgirepo/Gov Effectiveness.csv','ge')

def wdi(ind,name):
    d=pd.read_csv(f'/tmp/wdi/indicators/{ind}/data.csv'); d.columns=[c.strip() for c in d.columns]
    d=d.rename(columns={'Country Code':'code','Year':'year','Value':name})
    d['code']=d['code'].str.upper(); d['year']=d['year'].astype(int); return d[['code','year',name]]
exp=wdi('ne.exp.gnfs.zs','exp'); imp=wdi('ne.imp.gnfs.zs','imp')
bb=wdi('it.net.bbnd.p2','bband'); gdp=wdi('ny.gdp.pcap.pp.cd','gdppc')

vd=pyreadr.read_r('/tmp/demdata/data/vdem_simple.rda')['vdem_simple']
vd=vd[['country_text_id','year','v2x_polyarchy']].copy()
vd['code']=vd['country_text_id'].str.upper(); vd['year']=vd['year'].astype(int)
vd=vd.rename(columns={'v2x_polyarchy':'poly'})[['code','year','poly']]

m=corr.merge(gov,on=['code','year'],how='outer').merge(exp,on=['code','year'],how='left').merge(imp,on=['code','year'],how='left')
m=m.merge(bb,on=['code','year'],how='left').merge(gdp,on=['code','year'],how='left').merge(vd,on=['code','year'],how='left')
m['trade']=m['exp']+m['imp']
m['log_gdppc']=np.log(m['gdppc'])
p=m.dropna(subset=['cc','trade','log_gdppc']).copy()
# standardized connectedness composite (trade + broadband; broadband often missing early -> z on available)
for c in ['cc','ge','trade','bband','log_gdppc','poly']:
    p[c+'_z']=(p[c]-p[c].mean())/p[c].std()
out=[]
def w(s): out.append(s); print(s)
w("="*68); w("GLOBALIZATION / CONNECTEDNESS -> GOVERNANCE (replicating Sandholtz-Gray)")
w("DV=WGI Control of Corruption. N=%d, %d countries, %d-%d"%(len(p),p.code.nunique(),p.year.min(),p.year.max()))
w("corr(corruption_control, trade_openness)=%.3f ; corr(cc, log_gdppc)=%.3f"%(p.cc.corr(p.trade),p.cc.corr(p.log_gdppc)))

m1=smf.ols('cc ~ trade + log_gdppc + C(year)',data=p).fit(cov_type='cluster',cov_kwds={'groups':p['code']})
w("\n[M1] CC ~ trade openness + log GDPpc, year FE, cluster by country")
for v in ['trade','log_gdppc']: w("   %-10s b=%.4f (p=%.3g)"%(v,m1.params[v],m1.pvalues[v]))
w("   R2=%.3f"%m1.rsquared)

pf=p.set_index(['code','year'])
m2=PanelOLS.from_formula('cc ~ trade + log_gdppc + EntityEffects + TimeEffects',data=pf).fit(cov_type='clustered',cluster_entity=True)
w("\n[M2] TWO-WAY FE (country+year) -- WITHIN country: does rising openness track rising clean govt?")
for v in ['trade','log_gdppc']: w("   %-10s b=%.4f (p=%.3g)"%(v,m2.params[v],m2.pvalues[v]))

# broadband robustness (subset with bband)
pb=p.dropna(subset=['bband'])
m3=smf.ols('cc ~ bband + log_gdppc + C(year)',data=pb).fit(cov_type='cluster',cov_kwds={'groups':pb['code']})
w("\n[M3] DV=CC ~ broadband penetration + log GDPpc, year FE (N=%d): bband b=%.4f (p=%.3g)"%(len(pb),m3.params['bband'],m3.pvalues['bband']))

# KEY: non-democracy subsample (polyarchy < 0.5)
nd=p.dropna(subset=['poly'])
ndsub=nd[nd.poly<0.5]
m4=smf.ols('cc ~ trade + log_gdppc + C(year)',data=ndsub).fit(cov_type='cluster',cov_kwds={'groups':ndsub['code']})
w("\n[M4] NON-DEMOCRACY subsample (polyarchy<0.5), N=%d, %d countries: does connection discipline autocracies?"%(len(ndsub),ndsub.code.nunique()))
for v in ['trade','log_gdppc']: w("   %-10s b=%.4f (p=%.3g)"%(v,m4.params[v],m4.pvalues[v]))

# standardized betas (full sample)
mz=smf.ols('cc_z ~ trade_z + log_gdppc_z + C(year)',data=p).fit()
w("\n[std betas] trade=%.3f  income=%.3f"%(mz.params['trade_z'],mz.params['log_gdppc_z']))
# robustness DV gov effectiveness
mge=smf.ols('ge ~ trade + log_gdppc + C(year)',data=p.dropna(subset=['ge'])).fit(cov_type='cluster',cov_kwds={'groups':p.dropna(subset=['ge'])['code']})
w("[robust DV=Gov Effectiveness] trade b=%.4f (p=%.3g)"%(mge.params['trade'],mge.pvalues['trade']))
p.to_csv('/tmp/glob_panel.csv',index=False)
open('/tmp/results_glob.txt','w').write("\n".join(out))
print("\nSAVED")
