import pandas as pd, numpy as np
import matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from linearmodels.panel import PanelOLS
plt.rcParams['font.family']='DejaVu Serif'
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','Country Name':'cname','Year':'year','Value':name})
    d['code']=d['code'].str.upper(); d['year']=d['year'].astype(int); return d[['code','cname','year',name]]
mil=wdi('ms.mil.xpnd.gd.zs','mil'); edu=wdi('se.xpd.totl.gd.zs','edu')
print("military years:",mil.year.min(),"-",mil.year.max(),"; education years:",edu.year.min(),"-",edu.year.max())
m=mil.merge(edu[['code','year','edu']],on=['code','year'],how='inner').dropna(subset=['mil','edu'])
# drop tiny outliers
m=m[(m.mil<25)&(m.edu<15)]
out=[]
def w(s): out.append(s); print(s)
w("GUNS vs BUTTER: does a fall in military %GDP go with a RISE in education %GDP? (within-country)")
w("N=%d, %d countries, %d-%d"%(len(m),m.code.nunique(),m.year.min(),m.year.max()))
pf=m.set_index(['code','year'])
fe=PanelOLS.from_formula('edu ~ mil + EntityEffects + TimeEffects',data=pf).fit(cov_type='clustered',cluster_entity=True)
w("\n[Two-way FE] education ~ military: b=%.4f (se %.4f, p=%.3g)"%(fe.params['mil'],fe.std_errors['mil'],fe.pvalues['mil']))
w("   (negative b => when military share falls, education share rises = reallocation/peace dividend)")
# first-difference robustness
m=m.sort_values(['code','year']); m['d_mil']=m.groupby('code')['mil'].diff(); m['d_edu']=m.groupby('code')['edu'].diff()
fd=smf.ols('d_edu ~ d_mil',data=m.dropna(subset=['d_mil','d_edu'])).fit(cov_type='cluster',cov_kwds={'groups':m.dropna(subset=['d_mil','d_edu'])['code']})
w("[First differences] d_education ~ d_military: b=%.4f (p=%.3g)"%(fd.params['d_mil'],fd.pvalues['d_mil']))
open('/tmp/results_gunsbutter.txt','w').write("\n".join(out))

# Descriptive chart: military %GDP fall-and-rise for selected countries
sel={'DEU':'Germany','SWE':'Sweden','POL':'Poland','FRA':'France','GBR':'United Kingdom'}
fig,ax=plt.subplots(figsize=(7.2,4.6))
colors={'DEU':'#c0504d','SWE':'#4f81bd','POL':'#9bbb59','FRA':'#8064a2','GBR':'#7f7f7f'}
for c,nm in sel.items():
    d=mil[(mil.code==c)&(mil.year>=1990)].sort_values('year')
    if len(d): ax.plot(d.year,d.mil,label=nm,color=colors[c],lw=1.6)
ax.axhline(2.0,color='#bbb',ls='--',lw=0.9); ax.text(1991,2.06,'NATO 2% guideline',fontsize=8,color='#888')
ax.set_xlabel('Year',fontsize=10); ax.set_ylabel('Military spending (% of GDP)',fontsize=10)
ax.set_title('The peace dividend, and its reversal',fontsize=12.5,pad=10)
ax.legend(frameon=False,fontsize=9,loc='upper right')
for s in ['top','right']: ax.spines[s].set_visible(False)
plt.tight_layout(); plt.savefig('/tmp/fig_peace_dividend.png',dpi=170)
print("\nGermany mil%GDP last yrs:"); print(mil[(mil.code=='DEU')&(mil.year>=2018)][['year','mil']].to_string(index=False))
