"""
SHACKLES OF DEMOCRACY — Evidence of Decline: reproducible analysis
===================================================================
Builds a country-year panel (Freedom House + World Bank) and runs:
  PART 1  the decline trend (aggregate + within-country fixed effects, by region)
  PART 2  the economic-pressure -> backsliding panel (two-way FE, placebo)
  PART 3  illustrative forecast + figures

DATA SOURCES (public; fetched once via git clone in the build environment):
  * Freedom House "Freedom in the World" ratings 1972-2025
      repo: https://github.com/xmarquez/democracyData  (data-raw/fh_country_status_1973_2026.xlsx,
            fh_aggregate_scores_2003_2026.xlsx)  -- original: freedomhouse.org
  * World Bank World Development Indicators (real GDP growth, GDP per capita PPP,
      population, trade) repo: https://github.com/datasets/world-development-indicators

A prebuilt copy of the merged data ships as panel.csv, so this script can be run
analysis-only by pointing PANEL at that file. Requires: pandas numpy statsmodels
linearmodels country_converter openpyxl.
"""

# ====================== PART 0: BUILD PANEL ======================
import pandas as pd, numpy as np, country_converter as coco, warnings
warnings.filterwarnings("ignore")
DD="/tmp/dd/data-raw"; WDI="/tmp/wdi/indicators"

# ---------- 1. Freedom House PR/CL 1972-2025 (wide -> long) ----------
raw=pd.read_excel(f"{DD}/fh_country_status_1973_2026.xlsx",sheet_name="Country Ratings, Statuses ",header=None)
years=raw.iloc[1].tolist()            # year under review, every 3 cols
subhdr=raw.iloc[2].tolist()           # PR / CL / Status
rows=[]
data=raw.iloc[3:].reset_index(drop=True)
for c in range(1,len(years)):
    y=years[c]
    if pd.isna(y): continue
    try: y=int(y)
    except: continue
    # this col is PR; c+1 CL; c+2 Status  (only when subhdr[c]=='PR')
    if str(subhdr[c]).strip().upper()!="PR": continue
    block=data.iloc[:,[0,c,c+1,c+2]].copy()
    block.columns=["country","PR","CL","Status"]
    block["year"]=y
    rows.append(block)
fh=pd.concat(rows,ignore_index=True)
fh=fh[fh["country"].notna()]
for col in ["PR","CL"]:
    fh[col]=pd.to_numeric(fh[col],errors="coerce")
fh=fh.dropna(subset=["PR","CL"])
fh["fh_avg"]=(fh["PR"]+fh["CL"])/2                 # 1=most free .. 7=least free
fh["FreeScore"]=(7-fh["fh_avg"])/6*100             # 0..100, higher = freer
fh["iso3"]=coco.convert(fh["country"].tolist(),to="ISO3",not_found=None)
fh=fh[fh["iso3"].notna() & (fh["iso3"]!="not found")]
fh=fh.drop_duplicates(subset=["iso3","year"])
print("FH long:",fh.shape,"years",fh.year.min(),"-",fh.year.max(),"countries",fh.iso3.nunique())

# ---------- 2. World Bank indicators ----------
def wb(ind,name):
    d=pd.read_csv(f"{WDI}/{ind}/data.csv")
    d=d.rename(columns={"Country Code":"iso3","Year":"year","Value":name})[["iso3","year",name]]
    return d
g=wb("ny.gdp.mktp.kd.zg","gdp_growth")
pc=wb("ny.gdp.pcap.pp.cd","gdppc_ppp")
pop=wb("sp.pop.totl","pop")
ex=wb("ne.exp.gnfs.zs","exp_gdp"); im=wb("ne.imp.gnfs.zs","imp_gdp")

m=fh[["iso3","year","country","PR","CL","Status","fh_avg","FreeScore"]]
for d in [g,pc,pop,ex,im]:
    m=m.merge(d,on=["iso3","year"],how="left")
m["trade_open"]=m["exp_gdp"]+m["imp_gdp"]
m["log_gdppc"]=np.log(m["gdppc_ppp"])
m["log_pop"]=np.log(m["pop"])
m=m.sort_values(["iso3","year"])
m["gdp_growth_l1"]=m.groupby("iso3")["gdp_growth"].shift(1)
m["gdp_growth_l2"]=m.groupby("iso3")["gdp_growth"].shift(2)
m["FreeScore_l1"]=m.groupby("iso3")["FreeScore"].shift(1)
m["dFree"]=m["FreeScore"]-m["FreeScore_l1"]
m["gdp_growth_lead1"]=m.groupby("iso3")["gdp_growth"].shift(-1)  # placebo (future)

out="/sessions/serene-beautiful-maxwell/mnt/outputs/sd/panel.csv"
m.to_csv(out,index=False)
print("PANEL:",m.shape,"-> panel.csv")
print("non-null: growth",m.gdp_growth.notna().sum(),"gdppc",m.gdppc_ppp.notna().sum(),"trade",m.trade_open.notna().sum())
print(m[m.year==2024][["country","FreeScore","gdp_growth","log_gdppc"]].head(6).to_string())

# ============== PARTS 1-2: REGRESSIONS (writes results.txt) ==============
import pandas as pd, numpy as np, statsmodels.api as sm, statsmodels.formula.api as smf
from linearmodels.panel import PanelOLS
import country_converter as coco, warnings; warnings.filterwarnings("ignore")
m=pd.read_csv("/sessions/serene-beautiful-maxwell/mnt/outputs/sd/panel.csv")
L=[]; p=lambda s="": L.append(str(s))

# ============ PART 1: DECLINE TREND ============
p("="*66); p("PART 1 — IS DEMOCRACY DECLINING? (Freedom House FreeScore 0-100)"); p("="*66)
yr=m.groupby("year")["FreeScore"].mean()
# T1 full-period aggregate trend, Newey-West (HAC) SEs
d=yr.reset_index(); X=sm.add_constant(d["year"])
r=sm.OLS(d["FreeScore"],X).fit(cov_type="HAC",cov_kwds={"maxlags":4})
p(f"\n[T1] Global mean score ~ year, 1972-2025 (HAC se):")
p(f"  slope = {r.params['year']:+.4f} pts/yr  (se {r.bse['year']:.4f}, t {r.tvalues['year']:.2f}, p {r.pvalues['year']:.4f})")

# T2 RECESSION ERA 2006-2025, within-country (country FE), clustered SEs
rec=m[(m.year>=2006)&m.FreeScore.notna()].copy()
rec["yc"]=rec["year"]-2006
pe=rec.set_index(["iso3","year"])
fe=PanelOLS.from_formula("FreeScore ~ 1 + yc + EntityEffects",pe).fit(cov_type="clustered",cluster_entity=True)
b=fe.params["yc"]; se=fe.std_errors["yc"]
p(f"\n[T2] WITHIN-country trend 2006-2025 (country FE, cluster by country):")
p(f"  slope = {b:+.4f} pts/yr (se {se:.4f}, t {b/se:.2f}, p {fe.pvalues['yc']:.4g})  N={int(fe.nobs)}")
p(f"  => average country lost {abs(b)*19:.1f} points over 2006-2025." if b<0 else "  (rising)")

# T3 by region (use FH region from aggregate file)
reg=pd.read_excel("/tmp/dd/data-raw/fh_aggregate_scores_2003_2026.xlsx",sheet_name="FIW06-26")
reg.columns=[str(c).strip() for c in reg.columns]
reg["iso3"]=coco.convert(reg["Country/Territory"].tolist(),to="ISO3",not_found=None)
reg["iso3"]=reg["iso3"].apply(lambda x: x if isinstance(x,str) and x!="not found" else np.nan)
rmap=reg.dropna(subset=["iso3"]).drop_duplicates("iso3").set_index("iso3")["Region"].to_dict()
rec["region"]=rec["iso3"].map(rmap)
p("\n[T3] Within-country trend by region, 2006-2025 (pts/yr):")
for rg,sub in rec.dropna(subset=["region"]).groupby("region"):
    if sub.iso3.nunique()<5: continue
    pp=sub.set_index(["iso3","year"])
    try:
        f=PanelOLS.from_formula("FreeScore ~ 1 + yc + EntityEffects",pp).fit(cov_type="clustered",cluster_entity=True)
        p(f"  {rg:14s} {f.params['yc']:+.3f}/yr (p {f.pvalues['yc']:.3g}, {sub.iso3.nunique()} countries)")
    except Exception as e: p(f"  {rg}: {e}")

# share of countries declining vs improving 2006->2025
piv=m[m.year.isin([2006,2025])].pivot_table(index="iso3",columns="year",values="FreeScore")
piv=piv.dropna(); ch=piv[2025]-piv[2006]
p(f"\n[T4] Of {len(ch)} countries with data both years: {(ch<0).sum()} declined, {(ch>0).sum()} improved, {(ch==0).sum()} unchanged.")
p(f"     Mean change = {ch.mean():+.2f} pts; median {ch.median():+.2f}.")

# ============ PART 2: PANEL CAUSAL (economic pressure -> backsliding) ============
p("\n"+"="*66); p("PART 2 — DOES ECONOMIC PERFORMANCE PREDICT DEMOCRACY? (2-way FE)"); p("="*66)
panel=m.dropna(subset=["FreeScore","gdp_growth_l1","log_gdppc","log_pop"]).copy()
P=panel.set_index(["iso3","year"])
def show(name,formula,df):
    f=PanelOLS.from_formula(formula,df).fit(cov_type="clustered",cluster_entity=True,cluster_time=True)
    p(f"\n[{name}]  {formula}"); p(f"  N={int(f.nobs)}, entities={df.index.get_level_values(0).nunique()}, within-R2={f.rsquared_within:.3f}")
    for v in f.params.index:
        if v in ("Intercept",): continue
        star="***" if f.pvalues[v]<.01 else "**" if f.pvalues[v]<.05 else "*" if f.pvalues[v]<.1 else ""
        p(f"    {v:18s} {f.params[v]:+.4f} (se {f.std_errors[v]:.4f}, p {f.pvalues[v]:.3g}){star}")
    return f
show("P1 level","FreeScore ~ 1 + gdp_growth_l1 + log_gdppc + log_pop + EntityEffects + TimeEffects",P)
show("P2 +trade","FreeScore ~ 1 + gdp_growth_l1 + log_gdppc + log_pop + trade_open + EntityEffects + TimeEffects",P.reset_index().dropna(subset=["trade_open"]).set_index(["iso3","year"]))
dfc=m.dropna(subset=["dFree","gdp_growth_l1","log_gdppc"]).set_index(["iso3","year"])
show("P3 change","dFree ~ 1 + gdp_growth_l1 + log_gdppc + EntityEffects + TimeEffects",dfc)
# placebo: future growth should NOT predict today's democracy
plc=m.dropna(subset=["FreeScore","gdp_growth_lead1","gdp_growth_l1","log_gdppc","log_pop"]).set_index(["iso3","year"])
show("P4 placebo(+lead)","FreeScore ~ 1 + gdp_growth_l1 + gdp_growth_lead1 + log_gdppc + log_pop + EntityEffects + TimeEffects",plc)

open("/sessions/serene-beautiful-maxwell/mnt/outputs/sd/results.txt","w").write("\n".join(L))
print("\n".join(L))

# ============== PART 3: FORECAST + FIGURES ==============
import pandas as pd, numpy as np, statsmodels.api as sm, matplotlib
matplotlib.use("Agg"); import matplotlib.pyplot as plt
m=pd.read_csv("panel.csv")
INK="#1a1a2e"; RED="#c1121f"; BLU="#3a6ea5"; GLD="#d4a017"
plt.rcParams.update({"font.family":"DejaVu Serif","axes.edgecolor":"#555","figure.dpi":140})

# ---- status share (Free/PF/NF) per year ----
m["status"]=m["Status"].astype(str).str.strip().str[:2].str.upper()
def share(s): 
    t=m[m.year==s.name]; 
    return pd.Series({"Free":(s=="F").mean()})
yr=m.groupby("year")
shares=yr["status"].apply(lambda s:(s=="F").mean()*100)
meanfs=yr["FreeScore"].mean()

# ---- forecast from recession-era aggregate trend (2006-2025) ----
rec=meanfs[(meanfs.index>=2006)].reset_index(); rec.columns=["year","fs"]
X=sm.add_constant(rec["year"]); r=sm.OLS(rec["fs"],X).fit(cov_type="HAC",cov_kwds={"maxlags":4})
slope=r.params["year"]; 
fut=np.arange(2006,2076); Xf=sm.add_constant(fut)
pred=r.get_prediction(Xf); pm=pred.predicted_mean; ci=pred.conf_int(alpha=0.1)
# year crossing Partly-Free floor (33.3) and a "collapse" reference
def cross(level): 
    idx=np.where(pm<=level)[0]; return int(fut[idx[0]]) if len(idx) else None
print(f"Aggregate recession-era slope: {slope:+.3f}/yr (p {r.pvalues['year']:.3g})")
print(f"2025 mean={meanfs.loc[2025]:.1f}; projected 2050={pm[fut==2050][0]:.1f}, 2075={pm[-1]:.1f}")
print(f"Linear extrapolation crosses 'Partly Free' floor (33.3) ~ {cross(33.3)}")

# ===== FIG 1: long trend with two regimes =====
fig,ax=plt.subplots(figsize=(9,5))
ax.plot(meanfs.index,meanfs.values,color=INK,lw=2.4)
ax.axvspan(2006,2025,color=RED,alpha=0.08)
pk=meanfs.loc[2006]
ax.annotate(f"2006 peak ({pk:.1f})",(2006,pk),xytext=(1992,pk-4),color=RED,fontsize=9,
            arrowprops=dict(arrowstyle="->",color=RED))
ax.annotate(f"2025 ({meanfs.loc[2025]:.1f})",(2025,meanfs.loc[2025]),xytext=(2026,meanfs.loc[2025]-6),
            color=RED,fontsize=9)
ax.set_title("Global average freedom: the third wave, then the recession\nFreedom House score (0–100, higher = freer), country mean by year",fontsize=11,color=INK)
ax.set_ylabel("Mean Freedom Score"); ax.set_xlim(1972,2027); ax.set_ylim(41,66); ax.grid(alpha=.25)
fig.tight_layout(); fig.savefig("fig1_trend.png")

# ===== FIG 2: within-country regional slopes =====
regs={"Middle East":-0.581,"Eurasia":-0.502,"Africa":-0.389,"Americas":-0.367,"Europe":-0.304,"Asia":0.101}
fig,ax=plt.subplots(figsize=(8,4.4))
k=list(regs);v=[regs[x] for x in k]
ax.barh(k,v,color=[RED if x<0 else BLU for x in v])
ax.axvline(0,color=INK,lw=.8)
ax.set_title("Within-country change in freedom, 2006–2025 (points/year, country fixed effects)",fontsize=10.5,color=INK)
ax.set_xlabel("Annual change in Freedom Score"); ax.grid(axis="x",alpha=.25)
fig.tight_layout(); fig.savefig("fig2_regions.png")

# ===== FIG 3: forecast fan =====
fig,ax=plt.subplots(figsize=(9,5))
ax.plot(meanfs.index,meanfs.values,color=INK,lw=2,label="Observed")
fmask=fut>=2025
ax.plot(fut[fmask],pm[fmask],"--",color=RED,lw=2,label="Linear extrapolation of 2006–25 rate")
ax.fill_between(fut[fmask],ci[fmask,0],ci[fmask,1],color=RED,alpha=.12,label="90% band")
ax.axhline(75,color=GLD,ls=":",lw=1.3); ax.text(1973,76,"“Free” threshold",color=GLD,fontsize=8)
ax.axhline(33.3,color="#7a1f1f",ls=":",lw=1.3); ax.text(1973,34.3,"“Not Free” threshold",color="#7a1f1f",fontsize=8)
ax.set_title("If the recession-era rate simply continued…\n(illustrative extrapolation — NOT a prediction; see caveats)",fontsize=11,color=INK)
ax.set_ylabel("Mean Freedom Score"); ax.set_xlim(1972,2075); ax.set_ylim(20,75)
ax.legend(fontsize=8,loc="lower left"); ax.grid(alpha=.25)
fig.tight_layout(); fig.savefig("fig3_forecast.png")
print("charts saved")

# ====== PART 4: V-DEM + EXTRA IVs (build_v2) ======
import pandas as pd, numpy as np, pyreadr, warnings; warnings.filterwarnings("ignore")
WDI="/tmp/wdi/indicators"
m=pd.read_csv("panel.csv")
# --- V-Dem (slim) ---
v=pyreadr.read_r("/tmp/dd/data/vdem_simple.rda"); v=v[list(v.keys())[0]]
v=v[(v.year>=1972)][["country_text_id","year","v2x_libdem","v2x_polyarchy"]].rename(columns={"country_text_id":"iso3"})
v["libdem"]=v["v2x_libdem"]*100; v["polyarchy"]=v["v2x_polyarchy"]*100
v=v.dropna(subset=["iso3"]).drop_duplicates(["iso3","year"])
m=m.merge(v[["iso3","year","libdem","polyarchy"]],on=["iso3","year"],how="left")
# --- extra WDI IVs ---
def wb(ind,name):
    d=pd.read_csv(f"{WDI}/{ind}/data.csv").rename(columns={"Country Code":"iso3","Year":"year","Value":name})
    return d[["iso3","year",name]]
for ind,name in [("fp.cpi.totl.zg","inflation"),("sl.uem.totl.zs","unemp"),("ny.gdp.pcap.kd.zg","pcgrowth")]:
    m=m.merge(wb(ind,name),on=["iso3","year"],how="left")
m["inflation"]=m["inflation"].clip(-50,200)   # tame hyperinflation outliers
m=m.sort_values(["iso3","year"])
for c in ["inflation","unemp","pcgrowth"]:
    m[c+"_l1"]=m.groupby("iso3")[c].shift(1)
m["libdem_l1"]=m.groupby("iso3")["libdem"].shift(1); m["dlibdem"]=m["libdem"]-m["libdem_l1"]
m.to_csv("panel_v2.csv",index=False)
print("panel_v2:",m.shape)
print("non-null: libdem",m.libdem.notna().sum(),"inflation",m.inflation.notna().sum(),"unemp",m.unemp.notna().sum())
# FH vs V-Dem correlation
cc=m.dropna(subset=["FreeScore","libdem"])
print(f"corr(FH FreeScore, V-Dem libdem) = {cc.FreeScore.corr(cc.libdem):.3f}  (N={len(cc)})")

# ====== PART 5: V-DEM REPLICATION, SPEED-OF-DECLINE, EXPANDED PANEL (writes results_v2.txt) ======
import pandas as pd, numpy as np
from linearmodels.panel import PanelOLS
import warnings; warnings.filterwarnings("ignore")
m=pd.read_csv("panel_v2.csv"); L=[]; p=lambda s="":L.append(str(s))
def fe(formula,df,cl_time=False):
    return PanelOLS.from_formula(formula,df).fit(cov_type="clustered",cluster_entity=True,cluster_time=cl_time)
def slope(dv,df,var="yc"):
    f=fe(f"{dv} ~ 1 + {var} + EntityEffects",df.set_index(["iso3","year"]))
    return f.params[var],f.std_errors[var],f.pvalues[var],int(f.nobs)

# ===== A. V-DEM REPLICATION OF DECLINE =====
p("="*64); p("A. DECLINE IS NOT A FREEDOM HOUSE ARTIFACT (V-Dem replication)"); p("="*64)
for dv,lab,yend in [("FreeScore","Freedom House",2025),("libdem","V-Dem Liberal Dem",2024),("polyarchy","V-Dem Electoral Dem",2024)]:
    d=m[(m.year>=2006)&(m.year<=yend)&m[dv].notna()].copy(); d["yc"]=d.year-2006
    b,se,pv,n=slope(dv,d); 
    gm=m.groupby("year")[dv].mean(); peak=int(gm[(gm.index>=1990)].idxmax())
    p(f"  {lab:20s} within-country {b:+.3f}/yr (p {pv:.2g}); global mean peak={peak}, then {gm.loc[2006]:.1f}->{gm.loc[yend]:.1f}")

# ===== B. SPEED OF DECLINE: constant or accelerating? =====
p("\n"+"="*64); p("B. IS THE SPEED OF DECLINE CONSTANT SINCE 2006?"); p("="*64)
for dv,lab,yend in [("FreeScore","Freedom House",2025),("libdem","V-Dem LibDem",2024)]:
    d=m[(m.year>=2006)&(m.year<=yend)&m[dv].notna()].copy(); d["yc"]=d.year-2006; d["yc2"]=d["yc"]**2
    q=fe(f"{dv} ~ 1 + yc + yc2 + EntityEffects",d.set_index(["iso3","year"]))
    accel="ACCELERATING (concave down)" if q.params['yc2']<0 and q.pvalues['yc2']<.1 else ("DECELERATING" if q.params['yc2']>0 and q.pvalues['yc2']<.1 else "~constant (quadratic n.s.)")
    p(f"\n  [{lab}] quadratic: yc={q.params['yc']:+.3f} (p{q.pvalues['yc']:.2g}), yc^2={q.params['yc2']:+.4f} (p{q.pvalues['yc2']:.2g}) -> {accel}")
    # segmented slopes
    segs=[(2006,2012),(2013,2018),(2019,yend)]
    out=[]
    for a,b in segs:
        s=m[(m.year>=a)&(m.year<=b)&m[dv].notna()].copy(); s["yc"]=s.year-a
        bb,se,pv,n=slope(dv,s); out.append(f"{a}-{b}: {bb:+.3f}/yr (p{pv:.2g})")
    p("    segments  "+" | ".join(out))

# ===== C. WHERE IS THE ECONOMIC SIGNAL? (expanded multivariate) =====
p("\n"+"="*64); p("C. LOCATING THE CHANNEL: expanded two-way FE panel"); p("="*64)
def run(name,dv,rhs,need):
    df=m.dropna(subset=[dv]+need).set_index(["iso3","year"])
    f=fe(f"{dv} ~ 1 + {rhs} + EntityEffects + TimeEffects",df,cl_time=True)
    p(f"\n[{name}] DV={dv}  N={int(f.nobs)}, entities={df.index.get_level_values(0).nunique()}, within-R2={f.rsquared_within:.3f}")
    for v in f.params.index:
        if v=="Intercept": continue
        st="***" if f.pvalues[v]<.01 else "**" if f.pvalues[v]<.05 else "*" if f.pvalues[v]<.1 else ""
        p(f"    {v:16s} {f.params[v]:+.4f} (se {f.std_errors[v]:.4f}, p {f.pvalues[v]:.3g}){st}")
core="gdp_growth_l1 + inflation_l1 + trade_open + log_gdppc + log_pop"
run("C1 FH core",       "FreeScore",core,["gdp_growth_l1","inflation_l1","trade_open","log_gdppc","log_pop"])
run("C2 FH +unemp",     "FreeScore",core+" + unemp_l1",["gdp_growth_l1","inflation_l1","trade_open","log_gdppc","log_pop","unemp_l1"])
run("C3 V-Dem core",    "libdem",   core,["gdp_growth_l1","inflation_l1","trade_open","log_gdppc","log_pop"])
open("results_v2.txt","w").write("\n".join(L)); print("\n".join(L))

# ====== PART 6: ROBUSTNESS FIGURES ======
import pandas as pd, numpy as np, matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
m=pd.read_csv("panel_v2.csv")
INK="#1a1a2e"; RED="#c1121f"; BLU="#3a6ea5"; GLD="#d4a017"
plt.rcParams.update({"font.family":"DejaVu Serif","axes.edgecolor":"#555","figure.dpi":140})

# FIG 4: FH vs V-Dem (two independent measures, both decline)
fh=m.groupby("year")["FreeScore"].mean(); ld=m.groupby("year")["libdem"].mean()
fig,ax=plt.subplots(figsize=(9,5))
ax.plot(fh.index,fh.values,color=INK,lw=2.3,label="Freedom House (peak 2006)")
ax.plot(ld[ld.index<=2024].index,ld[ld.index<=2024].values,color=BLU,lw=2.3,label="V-Dem Liberal Democracy (peak 2012)")
for s,yr,c in [(fh,2006,INK),(ld,2012,BLU)]:
    ax.scatter([yr],[s.loc[yr]],color=RED,zorder=5,s=30)
ax.set_xlim(1972,2025); ax.grid(alpha=.25); ax.legend(fontsize=9,loc="lower right")
ax.set_ylabel("Global mean (0–100, higher = freer)")
ax.set_title("Two independent measures, one story (r = 0.92)\nGlobal average democracy rose, peaked, and has fallen since",fontsize=11,color=INK)
fig.tight_layout(); fig.savefig("fig4_fh_vs_vdem.png")

# FIG 5: SPEED — segmented within-country slopes (acceleration)
periods=["2006–2012","2013–2018","2019–2025"]
fhsl=[-0.176,-0.366,-0.380]; vdsl=[0.197,-0.081,-0.546]
x=np.arange(3); w=0.38
fig,ax=plt.subplots(figsize=(9,5))
ax.bar(x-w/2,fhsl,w,color=INK,label="Freedom House")
ax.bar(x+w/2,vdsl,w,color=BLU,label="V-Dem Liberal Democracy")
ax.axhline(0,color="#555",lw=.8)
for i,(a,b) in enumerate(zip(fhsl,vdsl)):
    ax.text(i-w/2,a-0.03 if a<0 else a+0.01,f"{a:+.2f}",ha="center",va="top" if a<0 else "bottom",fontsize=8.5,color=INK)
    ax.text(i+w/2,b-0.03 if b<0 else b+0.01,f"{b:+.2f}",ha="center",va="top" if b<0 else "bottom",fontsize=8.5,color=BLU)
ax.set_xticks(x); ax.set_xticklabels(periods); ax.grid(axis="y",alpha=.25)
ax.set_ylabel("Within-country change (points / year)")
ax.set_title("The decline is accelerating, not constant\nWithin-country rate of change by period (negative = losing freedom)",fontsize=11,color=INK)
ax.legend(fontsize=9,loc="lower left")
fig.tight_layout(); fig.savefig("fig5_speed.png")
print("fig4, fig5 saved")

# quick verification of segment slopes via independent mean-of-per-country method (FH)
from numpy.polynomial import polynomial as P
def seg(dv,a,b):
    d=m[(m.year>=a)&(m.year<=b)&m[dv].notna()]; sl=[]
    for iso,g in d.groupby("iso3"):
        if g.year.nunique()>=4: sl.append(np.polyfit(g.year,g[dv],1)[0])
    return np.mean(sl)
print("verify FH seg (per-country mean):",[round(seg("FreeScore",a,b),3) for a,b in [(2006,2012),(2013,2018),(2019,2025)]])

# ====== PART 7: SECURITY/CONFLICT — military burden panel + global battle deaths ======
import pandas as pd, numpy as np
from linearmodels.panel import PanelOLS
import warnings; warnings.filterwarnings("ignore")
m=pd.read_csv("panel_v2.csv")
mil=pd.read_csv("/tmp/wdi/indicators/ms.mil.xpnd.gd.zs/data.csv").rename(
    columns={"Country Code":"iso3","Year":"year","Value":"milex"})[["iso3","year","milex"]]
gdp=pd.read_csv("/tmp/wb_gdp/data/gdp.csv").rename(columns={"Country Code":"iso3","Year":"year","Value":"gdp_usd"})[["iso3","year","gdp_usd"]]
m=m.merge(mil,on=["iso3","year"],how="left").merge(gdp,on=["iso3","year"],how="left")
m=m.sort_values(["iso3","year"])
m["milex_l1"]=m.groupby("iso3")["milex"].shift(1)
m["d_milex"]=m.groupby("iso3")["milex"].diff()
L=[]; p=lambda s="":L.append(str(s))

# ===== GLOBAL: security stress vs democracy & its speed =====
p("="*64); p("D. GLOBAL SECURITY STRESS vs DEMOCRATIC DECLINE"); p("="*64)
g=m.dropna(subset=["milex","gdp_usd"]).copy()
# GDP-weighted world military burden (% of world GDP)
world=g.groupby("year").apply(lambda d:np.average(d.milex,weights=d.gdp_usd))
fs=m.groupby("year")["FreeScore"].mean(); ld=m.groupby("year")["libdem"].mean()
both=pd.concat({"milburden":world,"FreeScore":fs,"libdem":ld},axis=1).dropna()
both["dFS"]=both["FreeScore"].diff()  # annual SPEED of change
p(f"\nWorld military burden (% world GDP): 2000={both.milburden.get(2000,np.nan):.2f}, 2014={both.milburden.get(2014,np.nan):.2f}, 2024={both.milburden.get(2024,np.nan):.2f}")
# correlate security stress with the SPEED of decline, recession era
rec=both[both.index>=2006]
p(f"corr(world mil burden, annual change in FreeScore), 2006-2024: {rec['milburden'].corr(rec['dFS']):+.2f}  (negative => more security stress, faster decline)")
for a,b in [(2006,2012),(2013,2018),(2019,2024)]:
    s=both[(both.index>=a)&(both.index<=b)]
    p(f"  {a}-{b}: mean mil burden {s.milburden.mean():.2f}% | mean annual FreeScore change {s.dFS.mean():+.2f}/yr")

# ===== PANEL: does country security stress predict erosion? (2-way FE) =====
p("\n"+"="*64); p("E. WITHIN-COUNTRY: MILITARY/SECURITY BURDEN -> DEMOCRACY (2-way FE)"); p("="*64)
def run(name,dv,rhs,need):
    df=m.dropna(subset=[dv]+need).set_index(["iso3","year"])
    f=PanelOLS.from_formula(f"{dv} ~ 1 + {rhs} + EntityEffects + TimeEffects",df).fit(cov_type="clustered",cluster_entity=True,cluster_time=True)
    p(f"\n[{name}] DV={dv} N={int(f.nobs)} entities={df.index.get_level_values(0).nunique()} withinR2={f.rsquared_within:.3f}")
    for v in f.params.index:
        if v=="Intercept":continue
        st="***" if f.pvalues[v]<.01 else "**" if f.pvalues[v]<.05 else "*" if f.pvalues[v]<.1 else ""
        p(f"    {v:14s} {f.params[v]:+.4f} (se {f.std_errors[v]:.4f}, p {f.pvalues[v]:.3g}){st}")
run("E1 level","FreeScore","milex_l1 + gdp_growth_l1 + log_gdppc + log_pop",["milex_l1","gdp_growth_l1","log_gdppc","log_pop"])
run("E2 V-Dem","libdem","milex_l1 + gdp_growth_l1 + log_gdppc + log_pop",["milex_l1","gdp_growth_l1","log_gdppc","log_pop"])
run("E3 change","dFree","d_milex + gdp_growth_l1 + log_gdppc",["dFree","d_milex","gdp_growth_l1","log_gdppc"])
# recent era only (2013+), where decline is fast
mr=m[m.year>=2013]
df=mr.dropna(subset=["FreeScore","milex_l1","gdp_growth_l1","log_gdppc","log_pop"]).set_index(["iso3","year"])
f=PanelOLS.from_formula("FreeScore ~ 1 + milex_l1 + gdp_growth_l1 + log_gdppc + log_pop + EntityEffects + TimeEffects",df).fit(cov_type="clustered",cluster_entity=True,cluster_time=True)
p(f"\n[E4 2013+ era] DV=FreeScore N={int(f.nobs)}: milex_l1 {f.params['milex_l1']:+.4f} (p {f.pvalues['milex_l1']:.3g})")
m.to_csv("panel_v3.csv",index=False); open("results_v3.txt","w").write("\n".join(L)); print("\n".join(L))

# --- global battle deaths vs decline speed + fig6 ---
import pandas as pd, numpy as np, matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
m=pd.read_csv("panel_v3.csv")
# World battle deaths from state-based conflicts (UCDP via OWID), OWID_WRL
wd={1972:299083,1973:222940,1974:261646,1975:135684,1976:106722,1977:45550,1978:102573,1979:157289,1980:125492,
1981:197263,1982:234509,1983:228286,1984:249266,1985:228466,1986:235400,1987:194655,1988:152503,1989:55184,
1990:80297,1991:70353,1992:53406,1993:44949,1994:38517,1995:36742,1996:28885,1997:40340,1998:40218,1999:81047,
2000:78671,2001:23592,2002:21106,2003:23194,2004:19646,2005:12322,2006:20163,2007:19214,2008:28741,2009:34606,
2010:21088,2011:25080,2012:74071,2013:93488,2014:115972,2015:104435,2016:90624,2017:72023,2018:55031,2019:52610,
2020:73228,2021:199789,2022:276893,2023:131061,2024:128439}
bd=pd.Series(wd,name="deaths"); 
fs=m.groupby("year")["FreeScore"].mean(); ld=m.groupby("year")["libdem"].mean()
df=pd.concat({"deaths":bd,"FreeScore":fs,"libdem":ld},axis=1).dropna()
df["dFS"]=df["FreeScore"].diff(); df["logdeaths"]=np.log10(df["deaths"])
L=[];p=lambda s="":L.append(str(s))
p("="*64);p("F. GLOBAL SECURITY ISSUES (battle deaths) vs SPEED OF DECLINE");p("="*64)
rec=df[df.index>=2006]
p(f"\nWorld battle deaths: 2005={int(df.deaths[2005])}, 2014={int(df.deaths[2014])}, 2022={int(df.deaths[2022])}")
p(f"corr(log battle deaths, annual FreeScore change), 2006-2024 = {rec.logdeaths.corr(rec.dFS):+.2f}")
p(f"  (negative => years with more violence are years of faster democratic decline)")
p("\nPeriod | mean battle deaths/yr | mean decline speed (FreeScore/yr)")
for a,b in [(2000,2005),(2006,2012),(2013,2018),(2019,2024)]:
    s=df[(df.index>=a)&(df.index<=b)]
    p(f"  {a}-{b}:  {int(s.deaths.mean()):>9,}   {s.dFS.mean():+.2f}")
# lag: does a violent year precede faster decline next year?
rec2=df[df.index>=2005].copy(); rec2["dFS_next"]=rec2["FreeScore"].diff().shift(-1)
p(f"\ncorr(log deaths year t, decline speed year t+1), 2005-2023 = {rec2['logdeaths'].corr(rec2['dFS_next']):+.2f}")
open("results_v4.txt","w").write("\n".join(L)); print("\n".join(L))

# ===== FIG 6: global violence vs democracy =====
INK="#1a1a2e"; RED="#c1121f"; BLU="#3a6ea5"
plt.rcParams.update({"font.family":"DejaVu Serif","figure.dpi":140})
fig,ax=plt.subplots(figsize=(9.2,5))
ax.fill_between(df.index,df.deaths/1000,color=RED,alpha=.22,step="mid")
ax.plot(df.index,df.deaths/1000,color=RED,lw=1.3,label="Battle deaths (state-based conflicts, thousands)")
ax.set_ylabel("Global battle deaths (thousands)",color=RED); ax.tick_params(axis="y",labelcolor=RED)
ax.set_ylim(0,300)
ax2=ax.twinx()
ax2.plot(fs.index,fs.values,color=INK,lw=2.4,label="Freedom House score (right)")
ax2.scatter([2006],[fs.loc[2006]],color="k",zorder=5,s=25)
ax2.set_ylabel("Global mean Freedom Score",color=INK); ax2.set_ylim(40,66)
ax.axvspan(2011,2024,color="grey",alpha=.06)
ax.set_xlim(1990,2024)
ax.set_title("Global violence and democratic decline rose together after ~2011\nUCDP battle deaths vs Freedom House score (annual, world)",fontsize=11,color=INK)
l1,la1=ax.get_legend_handles_labels(); l2,la2=ax2.get_legend_handles_labels()
ax.legend(l1+l2,la1+la2,fontsize=8,loc="upper left")
fig.tight_layout(); fig.savefig("fig6_security.png"); print("fig6 saved")

# ====== PART 8: POWER TRANSITION THEORY (relative power, parity, contestation) ======
import pyreadr, pandas as pd, numpy as np, matplotlib; matplotlib.use("Agg"); import matplotlib.pyplot as plt
WDI="/tmp/wdi/indicators"; L=[]; p=lambda s="":L.append(str(s))
def wb(ind,name): return pd.read_csv(f"{WDI}/{ind}/data.csv").rename(columns={"Country Code":"iso3","Year":"year","Value":name})[["iso3","year",name]]
def ratio(df,val):
    w=df[df.iso3.isin(["USA","CHN"])].pivot(index="year",columns="iso3",values=val).dropna(); return (w["CHN"]/w["USA"])
g=wb("ny.gdp.mktp.cd","gdp"); pc=wb("ny.gdp.pcap.pp.cd","pc"); pop=wb("sp.pop.totl","popn")
ppp=pc.merge(pop,on=["iso3","year"]); ppp["gdp"]=ppp["pc"]*ppp["popn"]
mx=wb("ms.mil.xpnd.gd.zs","mx").merge(g,on=["iso3","year"]); mx["spend"]=mx["mx"]/100*mx["gdp"]
nmc=pyreadr.read_r("/tmp/ps/data/cow_nmc.rda")["cow_nmc"]; nmc["iso3"]=nmc.ccode.map({2:"USA",710:"CHN"})
r_cinc=ratio(nmc.rename(columns={}).assign(iso3=nmc.iso3)[nmc.iso3.notna()][["iso3","year","cinc"]],"cinc")
r_gdp=ratio(g,"gdp"); r_ppp=ratio(ppp,"gdp"); r_mil=ratio(mx,"spend")

p("="*66); p("G. POWER TRANSITION: CHINA/US RELATIVE POWER & PARITY (ratio, 1=parity)"); p("="*66)
def first_par(s): x=s[s>=1]; return int(x.index.min()) if len(x) else None
for nm,s in [("CINC composite",r_cinc),("GDP (PPP)",r_ppp),("GDP (market)",r_gdp),("Military spending",r_mil)]:
    p(f"  {nm:18s} 2000={s.get(2000,np.nan):.2f} 2010={s.get(2010,np.nan):.2f} 2016={s.get(2016,np.nan):.2f} 2023={s.get(2023,np.nan):.2f} | parity from {first_par(s)}")
gp=["USA","CHN","JPN","DEU","GBR","FRA","RUS","IND"]; gg=g[g.iso3.isin(gp)].pivot(index="year",columns="iso3",values="gdp")
usshare=gg["USA"]/gg.sum(axis=1); chshare=gg["CHN"]/gg.sum(axis=1)
p(f"\n  US share of 8-great-power GDP stays ~40% (1990:{usshare.get(1990)*100:.0f}% -> 2023:{usshare.get(2023)*100:.0f}%); China 2%->{chshare.get(2023)*100:.0f}%.")
p("  => 'US decline' is RELATIVE TO CHINA and measure-dependent: parity reached on PPP/CINC (~2014/1990s), NOT on market GDP/military.")

# Dissatisfaction proxy: active strategic rivalries + armed conflicts per year
tss=pyreadr.read_r("/tmp/ps/data/tss_rivalries.rda")["tss_rivalries"].copy()
tss["end"]=tss["end"].replace({9999:2024}); 
yrs=range(1946,2025)
nriv=pd.Series({y:int(((tss.start<=y)&(tss.end>=y)).sum()) for y in yrs})
usc=tss[((tss.ccode1==2)&(tss.ccode2==710))|((tss.ccode1==710)&(tss.ccode2==2))]
acd=pyreadr.read_r("/tmp/ps/data/ucdp_acd.rda")["ucdp_acd"]
nconf=acd.groupby("year")["conflict_id"].nunique()
p(f"\nDissatisfaction/contestation proxies: active strategic rivalries 2000={nriv.get(2000)} 2010={nriv.get(2010)} 2020={nriv.get(2020)}")
p(f"  US-China strategic rivalry on record: {len(usc)>0} (type: positional/ideological)")
p(f"  Active armed conflicts (UCDP): 2005={int(nconf.get(2005))} 2014={int(nconf.get(2014))} 2023={int(nconf.get(2023))}")

# Alignment: parity-approach vs democracy & violence (descriptive)
m=pd.read_csv("panel_v3.csv"); fs=m.groupby("year")["FreeScore"].mean()
wd={2005:12322,2006:20163,2007:19214,2008:28741,2009:34606,2010:21088,2011:25080,2012:74071,2013:93488,2014:115972,2015:104435,2016:90624,2017:72023,2018:55031,2019:52610,2020:73228,2021:199789,2022:276893,2023:131061,2024:128439}
bd=pd.Series(wd)
al=pd.concat({"ppp_ratio":r_ppp,"FreeScore":fs,"deaths":bd,"nconf":nconf},axis=1).loc[1995:2024]
p(f"\nDescriptive alignment 2008-2024 (entering the parity zone):")
p(f"  corr(China/US PPP power ratio, global FreeScore) = {al.loc[2000:].ppp_ratio.corr(al.loc[2000:].FreeScore):+.2f} (negative: closer to parity, less democracy)")
p(f"  corr(China/US PPP power ratio, active conflicts) = {al.loc[2000:].ppp_ratio.corr(al.loc[2000:].nconf):+.2f}")
open("results_v5.txt","w").write("\n".join(L)); print("\n".join(L))

# ===== FIG 7: relative power & parity =====
INK="#1a1a2e";RED="#c1121f";BLU="#3a6ea5";GLD="#d4a017";GRN="#2a7f4f"
plt.rcParams.update({"font.family":"DejaVu Serif","figure.dpi":140})
fig,ax=plt.subplots(figsize=(9.2,5.2))
for s,lab,c in [(r_cinc.loc[1950:],"CINC composite",GLD),(r_ppp.loc[1950:],"GDP (PPP)",RED),(r_gdp.loc[1960:],"GDP (market)",INK),(r_mil.loc[1989:],"Military spending",BLU)]:
    ax.plot(s.index,s.values,lw=2,label=lab,color=c)
ax.axhline(1,color="k",ls="--",lw=1); ax.text(1952,1.03,"PARITY",fontsize=8)
ax.axhspan(0.8,1.2,color="red",alpha=.06)
ax.set_ylim(0,2); ax.set_xlim(1950,2024); ax.grid(alpha=.25)
ax.set_ylabel("China ÷ United States"); ax.legend(fontsize=8.5,loc="upper left")
ax.set_title("How close is the challenger? China/US relative power by measure\nParity (=1) reached on PPP & CINC; not yet on market GDP or military spending",fontsize=11,color=INK)
fig.tight_layout(); fig.savefig("fig7_power_parity.png")

# ===== FIG 8: entering the parity zone vs democracy & violence =====
fig,ax=plt.subplots(figsize=(9.2,5.2))
ax.plot(r_ppp.loc[1990:].index,r_ppp.loc[1990:].values,color=RED,lw=2.4,label="China/US power ratio, PPP (left)")
ax.axhline(1,color=RED,ls=":",lw=1); ax.set_ylabel("China/US power ratio (PPP)",color=RED); ax.tick_params(axis="y",labelcolor=RED); ax.set_ylim(0,1.4)
ax.axvspan(2008,2024,color="grey",alpha=.07)
ax2=ax.twinx(); ax2.plot(fs.loc[1990:].index,fs.loc[1990:].values,color=INK,lw=2.4,label="Global Freedom score (right)")
ax2.set_ylabel("Global mean Freedom Score",color=INK); ax2.set_ylim(50,66)
ax.set_xlim(1990,2024); ax.set_title("Entering the parity zone as democracy declines\nChallenger approaches the hegemon (PPP parity ~2014) as global freedom falls",fontsize=11,color=INK)
h1,l1=ax.get_legend_handles_labels(); h2,l2=ax2.get_legend_handles_labels(); ax.legend(h1+h2,l1+l2,fontsize=8.5,loc="lower left")
fig.tight_layout(); fig.savefig("fig8_parity_democracy.png"); print("figs saved")
