# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Tür çeşitliliği ulusal deniz avını istikrarlı kılıyor mu? 168 ülkede portföy etkisi, türler arası eş zamanlılık ve av değişkenliği, 1950–2024
# Does species diversity stabilise national marine catches? Portfolio effect, synchrony among taxa and catch variability in 168 countries, 1950–2024
#
# GOMD Yayınlar · Makale 21 — https://gomd.pirireis.edu.tr/publications.php?id=rp21
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `f/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `f/out/` klasörlerinden okur, sonuçları `f/out/` ve şekilleri `f/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `f/` and runs them in order, each in its own process;
# it reads data from `data/` and `f/out/`, and writes results to `f/out/` and figures to `f/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow
# Çalıştırma / Run:  python rp21_national_catch_portfolio_effect.py
#
# Catch data: FAO Global Capture Production (FishStat, release 2026.1.0), © FAO, https://www.fao.org/fishery/en/collection/capture .
# GOMD, Piri Reis Üniversitesi ile T.C. Çevre, Şehircilik ve İklim Değişikliği Bakanlığı Çevre Yönetimi Genel Müdürlüğü tarafından geliştirilen açık veri platformudur.
# GOMD is an open data platform developed by Piri Reis University together with the Directorate General of Environmental Management of the Ministry of Environment, Urbanisation and Climate Change of Türkiye.
# ------------------------------------------------------------------------------
import os, subprocess, sys

MODULES = {}
ORDER = ['f/f_prep.py', 'f/f_analysis.py', 'f/f_extra.py', 'f/f_fig.py']
SKIP_IF_EXISTS = {'f/f_prep.py': 'f/out/marine.parquet'}   # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package

# ==============================================================================
# f/f_prep.py
# ==============================================================================
MODULES['f/f_prep.py'] = r'''
# FAO Kuresel Av Uretimi (1950-2024): deniz alanlari, balik + kabuklu + yumusakca + diger omurgasiz; ulke x takson x yil
import pandas as pd, numpy as np
D='data/fao/x/'
q=pd.read_csv(D+'Capture_Quantity.csv',dtype={'COUNTRY.UN_CODE':str,'AREA.CODE':str}); q.columns=['un','sp','area','meas','yr','val','st']
sp=pd.read_csv(D+'CL_FI_SPECIES_GROUPS.csv',usecols=['3A_Code','Name_En','Scientific_Name','Major_Group','ISSCAAP_Group_En','Taxonomic_Code']).rename(columns={'3A_Code':'sp','Name_En':'sp_name','ISSCAAP_Group_En':'isscaap','Taxonomic_Code':'tax'})
ar=pd.read_csv(D+'CL_FI_WATERAREA_GROUPS.csv',dtype={'Code':str},usecols=['Code','Name_En','InlandMarine_Group_En']).rename(columns={'Code':'area','Name_En':'area_name','InlandMarine_Group_En':'im'})
co=pd.read_csv(D+'CL_FI_COUNTRY_GROUPS.csv',dtype={'UN_Code':str},usecols=['UN_Code','ISO3_Code','Name_En','Continent_Group_En','EcoClass_Group_En','GeoRegion_Group_En']).rename(columns={'UN_Code':'un','ISO3_Code':'iso','Name_En':'country','Continent_Group_En':'cont','EcoClass_Group_En':'eco','GeoRegion_Group_En':'georeg'})
q=q.merge(sp,on='sp',how='left').merge(ar,on='area',how='left').merge(co,on='un',how='left')
m=q[(q.meas=='Q_tlw')&(q.im=='Marine areas')&q.Major_Group.isin(['PISCES','CRUSTACEA','MOLLUSCA','INVERTEBRATA AQUATICA'])].copy()
m['species_level']=~m.tax.astype(str).str.contains('X')
m['official']=(m.st=='A')
print('rows',len(m),'countries',m.un.nunique(),'taxa',m.sp.nunique(),'years',m.yr.min(),m.yr.max())
print('world Mt',(m.groupby('yr').val.sum()/1e6).iloc[[0,10,20,30,40,50,60,70,74]].round(1).to_dict())
print('species-level share',(m[m.species_level].groupby('yr').val.sum()/m.groupby('yr').val.sum()).iloc[::10].round(3).to_dict())
print('official share',(m[m.official].groupby('yr').val.sum()/m.groupby('yr').val.sum()).iloc[::10].round(3).to_dict())
m[['un','iso','country','cont','eco','georeg','sp','sp_name','isscaap','Major_Group','species_level','official','area','area_name','yr','val']].to_parquet('f/out/marine.parquet')
print(m.groupby('isscaap').val.sum().sort_values(ascending=False).head(15)/1e6)
print(m.eco.value_counts(dropna=False).to_dict()); print(m.cont.value_counts(dropna=False).to_dict())
print(m[m.country.str.contains('Sov|Yugo|Czecho|Serbia and|Sudan|Ethiopia|Netherlands Ant',na=False)].groupby('country').yr.agg(['min','max']))
'''

# ==============================================================================
# f/f_analysis.py
# ==============================================================================
MODULES['f/f_analysis.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.formula.api as smf, statsmodels.api as sm
from scipy import stats
m=pd.read_parquet('f/out/marine.parquet')
WIN={'1950–1974':(1950,1974),'1975–1999':(1975,1999),'2000–2024':(2000,2024)}; WN=list(WIN)
SMALLPEL='Herrings, sardines, anchovies'
R={'n_rows':int(len(m)),'n_countries':int(m.un.nunique()),'n_taxa':int(m.sp.nunique()),'years':[int(m.yr.min()),int(m.yr.max())],
   'world':{int(y):float(v) for y,v in m.groupby('yr').val.sum().items()},
   'sp_share':{int(y):float(v) for y,v in (m[m.species_level].groupby('yr').val.sum()/m.groupby('yr').val.sum()).items()}}
def portfolio(W):
    """W: yil x varlik matrisi (ton). Birinci farklarla portfoy ayristirmasi."""
    W=W.loc[:,(W.sum()>0)]; X=W.sum(axis=1); mT=X.mean(); d=W.diff().dropna(); dT=X.diff().dropna()
    sig=d.std(ddof=1); sT=dT.std(ddof=1); mi=W.mean()
    S=sig.sum(); phi=sT**2/S**2; phi0=(sig**2).sum()/S**2; rho=(phi-phi0)/(1-phi0) if phi0<1 else np.nan
    p=mi/mi.sum(); p=p[p>0]; H=-(p*np.log(p)).sum(); neff=np.exp(H); simp=1/(p**2).sum()
    cvT=sT/mT; cvw=S/mT
    # ortalama-varyans (Taylor) portfoy etkisi
    ok=(mi>0)&(sig>0)
    z=a=np.nan; pe_mv=np.nan
    if ok.sum()>=6:
        r=stats.linregress(np.log(mi[ok]),np.log(sig[ok])); z=r.slope; a=r.intercept; s_pred=np.exp(a+z*np.log(mT)); pe_mv=1-sT/s_pred
    # logaritmik buyume oynakligi (toplam)
    g=np.log(X).diff().dropna(); sdg=g.std(ddof=1)
    return dict(n_assets=int(W.shape[1]),mean=float(mT),cvT=float(cvT),cvw=float(cvw),phi=float(phi),phi0=float(phi0),rho=float(rho),pe=float(1-np.sqrt(phi)),neff=float(neff),simpson=float(simp),top1=float(p.max()),z=float(z),pe_mv=float(pe_mv),sdg=float(sdg))
rows=[]
for (un,g) in m.groupby('un'):
    info=g.iloc[0]
    for wn,(a,b) in WIN.items():
        x=g[(g.yr>=a)&(g.yr<=b)]
        if x.yr.nunique()<(b-a+1): continue
        tot=x.groupby('yr').val.sum()
        if (tot<=0).any() or tot.mean()<1000: continue
        for lev,key in(('taxa','sp'),('group','isscaap')):
            W=x.pivot_table(index='yr',columns=key,values='val',aggfunc='sum').reindex(range(a,b+1)).fillna(0)
            W=W.loc[:,(W.sum()>0)]
            if W.shape[1]<3 or W.sum(axis=1).diff().dropna().std()==0: continue
            if x[x.official].val.sum()/x.val.sum()<0.5: continue
            r=portfolio(W); r.update(un=un,iso=info.iso,country=info.country,cont=info.cont,eco=info.eco,win=wn,level=lev)
            r['sp_share']=float(x[x.species_level].val.sum()/x.val.sum()); r['sp_share_first']=float(x[(x.yr<=a+4)&x.species_level].val.sum()/x[x.yr<=a+4].val.sum()); r['sp_share_last']=float(x[(x.yr>=b-4)&x.species_level].val.sum()/x[x.yr>=b-4].val.sum())
            r['official']=float(x[x.official].val.sum()/x.val.sum()); r['smallpel']=float(x[x.isscaap==SMALLPEL].val.sum()/x.val.sum()); r['n_areas']=int(x.area.nunique())
            rows.append(r)
D=pd.DataFrame(rows); D.to_csv('f/out/country_window.csv',index=False)
T=D[D.level=='taxa'].copy(); Gp=D[D.level=='group'].copy()
R['n_cw']={w:int((T.win==w).sum()) for w in WN}; R['n_cw_total']=int(len(T)); R['n_c']=int(T.un.nunique())
def summ(x):
    return dict(n=int(len(x)),cvT=float(x.cvT.median()),cvw=float(x.cvw.median()),phi=float(x.phi.median()),sqphi=float(np.sqrt(x.phi).median()),pe=float(x.pe.median()),pe_p25=float(x.pe.quantile(.25)),pe_p75=float(x.pe.quantile(.75)),rho=float(x.rho.median()),neff=float(x.neff.median()),n_assets=float(x.n_assets.median()),pe_mv=float(x.pe_mv.median()),z=float(x.z.median()),sdg=float(x.sdg.median()),
                pe_w=float(np.average(x.pe,weights=x['mean'])),cvT_w=float(np.average(x.cvT,weights=x['mean'])),rho_w=float(np.average(x.rho,weights=x['mean'])),share_pos_pe=float((x.pe>0).mean()),share_rho_neg=float((x.rho<0).mean()),pe_mv_pos=float((x.pe_mv>0).mean()))
R['summary']={lev:{w:summ(d[d.win==w]) for w in WN} for lev,d in(('taxa',T),('group',Gp))}
R['summary']['taxa']['all']=summ(T); R['summary']['group']['all']=summ(Gp)
# dengeli panel: uc donemde de bulunan ulkeler
bal=T.groupby('un').win.nunique(); balu=bal[bal==3].index; TB=T[T.un.isin(balu)]; GB=Gp[Gp.un.isin(balu)]
R['balanced']={'n':int(len(balu)),'taxa':{w:summ(TB[TB.win==w]) for w in WN},'group':{w:summ(GB[GB.win==w]) for w in WN}}
# eslestirilmis degisim testleri (ilk -> son donem)
def paired(d,col,f=lambda v:v):
    w=d.pivot(index='un',columns='win',values=col); a=f(w[WN[0]]); b=f(w[WN[2]]); t=stats.wilcoxon(b-a); return dict(med_first=float(a.median()),med_last=float(b.median()),med_diff=float((b-a).median()),p=float(t.pvalue),share_up=float(((b-a)>0).mean()))
R['paired']={lev:{c:paired(d,c) for c in('pe','rho','cvT','cvw','neff','phi','sdg','pe_mv')} for lev,d in(('taxa',TB),('group',GB))}
# ---- regresyonlar: log cvT = 0.5 log phi + log cvw
def regs(d,extra=''):
    d=d.copy(); d['lcv']=np.log(d.cvT); d['lsq']=0.5*np.log(d.phi); d['lcw']=np.log(d.cvw); d['ln']=np.log(d.neff); d['lm']=np.log(d['mean']); d=d.replace([np.inf,-np.inf],np.nan).dropna(subset=['lcv','lsq','lcw','ln','lm','smallpel']).copy(); d['g']=pd.factorize(d.un)[0]
    out={}
    for dep in('lcv','lsq','lcw'):
        mm=smf.ols(dep+' ~ ln + lm + smallpel + C(win)'+extra,d).fit(cov_type='cluster',cov_kwds={'groups':d.g})
        out[dep]={k:dict(b=float(mm.params[k]),se=float(mm.bse[k]),p=float(mm.pvalues[k])) for k in('ln','lm','smallpel')}; out[dep]['r2']=float(mm.rsquared); out[dep]['n']=int(mm.nobs)
    # ulke sabit etkili (ulke ici degisim)
    mm=smf.ols('lcv ~ ln + lm + smallpel + C(win) + C(un)',d).fit(cov_type='cluster',cov_kwds={'groups':d.g}); out['fe']={k:dict(b=float(mm.params[k]),se=float(mm.bse[k]),p=float(mm.pvalues[k])) for k in('ln','lm','smallpel')}
    return out
R['reg']={'taxa':regs(T),'group':regs(Gp)}
# saglamlik: tur duzeyi raporlama payi yuksek ve kararli olanlar
hi=T[(T.sp_share_first>=0.8)&(T.sp_share_last>=0.8)]; R['robust']={'hi_species':dict(n=int(len(hi)),**{k:v for k,v in summ(hi).items() if k in('pe','rho','cvT','neff','pe_mv')}),'reg_hi':regs(hi)['lcv']}
off=T[T.official>=0.9]; R['robust']['official']=dict(n=int(len(off)),**{k:v for k,v in summ(off).items() if k in('pe','rho','cvT','neff')}); R['robust']['reg_off']=regs(off)['lcv']
big=T[T['mean']>=100000]; R['robust']['big']=dict(n=int(len(big)),**{k:v for k,v in summ(big).items() if k in('pe','rho','cvT','neff')}); R['robust']['reg_big']=regs(big)['lcv']
# ---- gelir ve kita
R['by_eco']={str(k):summ(x) for k,x in T[T.win==WN[2]].groupby('eco')}; R['by_cont']={str(k):summ(x) for k,x in T[T.win==WN[2]].groupby('cont')}
# ---- en yuksek/dusuk ulkeler (son donem, ortalama av >= 100 bin ton)
L=T[(T.win==WN[2])&(T['mean']>=100000)].sort_values('cvT'); R['n_big_last']=int(len(L))
f=lambda r:[r.country,float(r['mean']),float(r.cvT),float(r.pe),float(r.neff),float(r.rho),float(r.smallpel),float(r.top1)]
R['most_stable']=[f(r) for _,r in L.head(10).iterrows()]; R['least_stable']=[f(r) for _,r in L.tail(10).iloc[::-1].iterrows()]
R['largest']=[f(r) for _,r in T[T.win==WN[2]].sort_values('mean',ascending=False).head(15).iterrows()]
# ---- kuresel duzey: ulkeler arasi portfoy (dunya toplami) ve FAO alanlari arasi
R['global']={}
for wn,(a,b) in WIN.items():
    x=m[(m.yr>=a)&(m.yr<=b)]
    Wc=x.pivot_table(index='yr',columns='un',values='val',aggfunc='sum').fillna(0); Wa=x.pivot_table(index='yr',columns='area',values='val',aggfunc='sum').fillna(0); Ws=x.pivot_table(index='yr',columns='sp',values='val',aggfunc='sum').fillna(0); Wg=x.pivot_table(index='yr',columns='isscaap',values='val',aggfunc='sum').fillna(0)
    R['global'][wn]={k:{kk:vv for kk,vv in portfolio(W).items() if kk in('cvT','cvw','phi','pe','rho','neff','n_assets','mean','sdg','pe_mv','z')} for k,W in(('country',Wc),('area',Wa),('taxa',Ws),('group',Wg))}
# anchoveta haric dunya
anc=m[m.sp_name.str.contains('Anchoveta',na=False)].sp.unique(); R['anchoveta_codes']=[str(a) for a in anc]
for wn,(a,b) in WIN.items():
    x=m[(m.yr>=a)&(m.yr<=b)&~m.sp.isin(anc)]; Ws=x.pivot_table(index='yr',columns='sp',values='val',aggfunc='sum').fillna(0); R['global'][wn]['taxa_exanch']={kk:vv for kk,vv in portfolio(Ws).items() if kk in('cvT','cvw','phi','pe','rho','neff')}
# ---- Taylor yasasi: tur duzeyi sigma ~ m^z (ulke-donem sabit etkili)
rowsz=[]
for (un,g) in m.groupby('un'):
    for wn,(a,b) in WIN.items():
        if not ((T.un==un)&(T.win==wn)).any(): continue
        x=g[(g.yr>=a)&(g.yr<=b)]; W=x.pivot_table(index='yr',columns='sp',values='val',aggfunc='sum').reindex(range(a,b+1)).fillna(0); mi=W.mean(); sg=W.diff().dropna().std(ddof=1); ok=(mi>0)&(sg>0)
        rowsz.append(pd.DataFrame({'cw':un+wn,'lm':np.log(mi[ok]),'ls':np.log(sg[ok])}))
Z=pd.concat(rowsz); Z['lmd']=Z.lm-Z.groupby('cw').lm.transform('mean'); Z['lsd']=Z.ls-Z.groupby('cw').ls.transform('mean'); r=stats.linregress(Z.lmd,Z.lsd); R['taylor']=dict(z=float(r.slope),se=float(r.stderr),n=int(len(Z)),r2=float(r.rvalue**2))
def clean(o):
    if isinstance(o,dict): return {str(k):clean(v) for k,v in o.items()}
    if isinstance(o,(list,tuple)): return [clean(v) for v in o]
    if isinstance(o,(float,np.floating)): return None if not np.isfinite(o) else float(o)
    if isinstance(o,np.integer): return int(o)
    return o
json.dump(clean(R),open('f/out/results.json','w'),ensure_ascii=False,indent=1)

KS=('n','cvT','cvw','sqphi','pe','rho','neff','n_assets','pe_mv','z','share_rho_neg','pe_mv_pos','cvT_w','pe_w')
pr=lambda d:{k:(round(d[k],3) if d.get(k) is not None else None) for k in KS if k in d}
print('n_cw',R['n_cw'],R['n_c'])
for lev in('taxa','group'):
    for w in WN+['all']: print('SUM',lev,w,pr(clean(R['summary'][lev][w])))
print('BAL n',R['balanced']['n'])
for lev in('taxa','group'):
    for w in WN: print('BAL',lev,w,pr(clean(R['balanced'][lev][w])))
print('PAIRED',json.dumps(clean(R['paired']))[:3000])
print('REG',json.dumps(clean(R['reg']))[:3500])
print('ROB',json.dumps(clean(R['robust']))[:2500])
for k,v in R['by_eco'].items(): print('ECO',k,pr(clean(v)))
for k,v in R['by_cont'].items(): print('CONT',k,pr(clean(v)))
print('STABLE',[[x[0],round(x[1]/1e3),round(x[2],3),round(x[3],2),round(x[4],1),round(x[6],2)] for x in R['most_stable']])
print('UNSTABLE',[[x[0],round(x[1]/1e3),round(x[2],3),round(x[3],2),round(x[4],1),round(x[6],2)] for x in R['least_stable']])
for w in WN: print('GLOBAL',w,{k:{kk:round(vv,3) for kk,vv in v.items() if kk in('cvT','cvw','pe','rho','neff','n_assets','pe_mv')} for k,v in clean(R['global'][w]).items()})
print('TAYLOR',R['taylor'],R['anchoveta_codes'])
'''

# ==============================================================================
# f/f_extra.py
# ==============================================================================
MODULES['f/f_extra.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.formula.api as smf, statsmodels.api as sm
from scipy import stats
m=pd.read_parquet('f/out/marine.parquet'); D=pd.read_csv('f/out/country_window.csv',dtype={'un':str}); R=json.load(open('f/out/results.json'))
WIN={'1950–1974':(1950,1974),'1975–1999':(1975,1999),'2000–2024':(2000,2024)}; WN=list(WIN)
T=D[D.level=='taxa'].copy(); G=D[D.level=='group'].copy(); X={}
# 1. eslestirilmis testler (eksikler atilarak) ve degisimin ayristirilmasi
def paired(d,col):
    w=d.pivot(index='un',columns='win',values=col)[[WN[0],WN[2]]].dropna(); a,b=w[WN[0]],w[WN[2]]; t=stats.wilcoxon(b-a); return dict(n=int(len(w)),med_first=float(a.median()),med_last=float(b.median()),med_diff=float((b-a).median()),p=float(t.pvalue),share_up=float(((b-a)>0).mean()))
for lev,d in(('taxa',T),('group',G)):
    bal=d.groupby('un').win.nunique(); db=d[d.un.isin(bal[bal==3].index)]
    X['paired_'+lev]={c:paired(db,c) for c in('pe','rho','cvT','cvw','neff','phi','sdg','pe_mv','n_assets','sp_share')}
    w=db.assign(lcv=np.log(db.cvT),lsq=0.5*np.log(db.phi),lcw=np.log(db.cvw)).pivot(index='un',columns='win',values=['lcv','lsq','lcw'])
    X['decomp_'+lev]={k:float((w[k][WN[2]]-w[k][WN[0]]).dropna().median()) for k in('lcv','lsq','lcw')}; X['decomp_'+lev]['mean']={k:float((w[k][WN[2]]-w[k][WN[0]]).dropna().mean()) for k in('lcv','lsq','lcw')}
# 2. ek regresyonlar: alan sayisi, log buyume sd, en buyuk takson payi
def prep(d):
    d=d.copy(); d['lcv']=np.log(d.cvT); d['lsd']=np.log(d.sdg); d['ln']=np.log(d.neff); d['lm']=np.log(d['mean']); d['la']=np.log(d.n_areas); d['lsq']=0.5*np.log(d.phi); d['lcw']=np.log(d.cvw)
    d=d.replace([np.inf,-np.inf],np.nan).dropna(subset=['lcv','lsd','ln','lm','la','smallpel','lsq','lcw']); d['g']=pd.factorize(d.un)[0]; return d
def fit(f,d,keys):
    mm=smf.ols(f,d).fit(cov_type='cluster',cov_kwds={'groups':d.g}); o={k:dict(b=float(mm.params[k]),se=float(mm.bse[k]),p=float(mm.pvalues[k])) for k in keys}; o['r2']=float(mm.rsquared); o['n']=int(mm.nobs); return o
dT=prep(T); dG=prep(G)
X['reg_areas']=fit('lcv ~ ln + lm + smallpel + la + C(win)',dT,('ln','lm','smallpel','la'))
X['reg_sdg']=fit('lsd ~ ln + lm + smallpel + C(win)',dT,('ln','lm','smallpel'))
X['reg_top1']=fit('lcv ~ top1 + lm + smallpel + C(win)',dT,('top1','lm','smallpel'))
X['reg_cont']=fit('lcv ~ ln + lm + smallpel + C(win) + C(cont)',dT.dropna(subset=['cont']),('ln','lm','smallpel'))
X['reg_spshare']=fit('lcv ~ ln + lm + smallpel + sp_share + C(win)',dT,('ln','lm','smallpel','sp_share'))
for w in WN: X['reg_'+w]=fit('lcv ~ ln + lm + smallpel',dT[dT.win==w],('ln','lm','smallpel'))
# yalniz buyukluk
X['reg_size_only']=fit('lcv ~ lm + C(win)',dT,('lm',))
# 3. sok sikligi: toplam avda yildan yila %20'den fazla dusus
rows=[]
for (un,g) in m.groupby('un'):
    for wn,(a,b) in WIN.items():
        if not ((T.un==un)&(T.win==wn)).any(): continue
        tot=g[(g.yr>=a)&(g.yr<=b)].groupby('yr').val.sum().reindex(range(a,b+1)); gl=np.log(tot).diff().dropna()
        rows.append(dict(un=un,win=wn,n20=int((gl<np.log(0.8)).sum()),n10=int((gl<np.log(0.9)).sum()),nup20=int((gl>np.log(1.2)).sum()),ny=int(len(gl)),worst=float(gl.min())))
S=pd.DataFrame(rows).merge(T[['un','win','neff','mean','smallpel','cvT','pe']],on=['un','win'])
X['shock']=dict(n=int(len(S)),rate20=float(S.n20.sum()/S.ny.sum()),rate10=float(S.n10.sum()/S.ny.sum()),rateup20=float(S.nup20.sum()/S.ny.sum()),share_any20=float((S.n20>0).mean()),worst_med=float(S.worst.median()))
S['q']=pd.qcut(S.neff,4,labels=False); X['shock']['by_q']=[dict(q=int(k)+1,neff=float(x.neff.median()),rate20=float(x.n20.sum()/x.ny.sum()),rate10=float(x.n10.sum()/x.ny.sum()),cvT=float(x.cvT.median()),n=int(len(x))) for k,x in S.groupby('q')]
S['ln']=np.log(S.neff); S['lm']=np.log(S['mean']); S['g']=pd.factorize(S.un)[0]
mm=sm.GLM(S.n20,sm.add_constant(S[['ln','lm','smallpel']]),family=sm.families.Poisson(),offset=np.log(S.ny)).fit(cov_type='cluster',cov_kwds={'groups':S.g})
X['shock']['pois']={k:dict(irr=float(np.exp(mm.params[k])),lo=float(np.exp(mm.params[k]-1.96*mm.bse[k])),hi=float(np.exp(mm.params[k]+1.96*mm.bse[k])),p=float(mm.pvalues[k])) for k in('ln','lm','smallpel')}
S['wq']=S.win; X['shock']['by_win']={w:dict(rate20=float(x.n20.sum()/x.ny.sum()),rate10=float(x.n10.sum()/x.ny.sum())) for w,x in S.groupby('win')}
# 4. kuresel yillik degisim: dunya toplami log buyume sd, anchoveta paylari
tot=m.groupby('yr').val.sum(); anc=m[m.sp=='VET'].groupby('yr').val.sum().reindex(tot.index).fillna(0); ex=tot-anc
X['world']={w:dict(sd_g=float(np.log(tot.loc[a:b]).diff().std()),sd_g_ex=float(np.log(ex.loc[a:b]).diff().std()),anch_share=float(anc.loc[a:b].sum()/tot.loc[a:b].sum()),mean=float(tot.loc[a:b].mean()),var_share_anch=float(1-(ex.loc[a:b].diff().var()/tot.loc[a:b].diff().var()))) for w,(a,b) in WIN.items()}
X['world']['peak']=[int(tot.idxmax()),float(tot.max())]; X['world']['last']=[int(tot.index[-1]),float(tot.iloc[-1])]
# 5. ISSCAAP gruplarinin oynakligi (dunya, 2000-2024): grup basina cv ve pay
x=m[(m.yr>=2000)&(m.yr<=2024)]; W=x.pivot_table(index='yr',columns='isscaap',values='val',aggfunc='sum').fillna(0); mi=W.mean(); cv=W.diff().dropna().std()/mi
X['groups']=[[k,float(mi[k]),float(mi[k]/mi.sum()),float(cv[k])] for k in mi.sort_values(ascending=False).head(12).index]
# 6. ulke ornekleri: Peru, Norvec, Izlanda, Turkiye son donem + ilk donem
def ex_(c):
    r=T[(T.country==c)].set_index('win'); return {w:{k:float(r.loc[w,k]) for k in('mean','cvT','cvw','pe','neff','rho','smallpel','top1','n_assets')} for w in r.index}
X['cases']={c:ex_(c) for c in('Peru','Norway','Iceland','Türkiye','Japan','China','Indonesia','United States of America','Chile','India')}
# 7. korelasyon: pe ile neff, rho ile neff (Spearman)
X['corr']=dict(pe_neff=float(stats.spearmanr(T.pe,T.neff)[0]),rho_neff=float(stats.spearmanr(T.rho,T.neff)[0]),cv_neff=float(stats.spearmanr(T.cvT,T.neff)[0]),cv_mean=float(stats.spearmanr(T.cvT,T['mean'])[0]),neff_mean=float(stats.spearmanr(T.neff,T['mean'])[0]))
# 8. sifir korelasyon varsayimi altinda beklenen PE ile gozlenen PE
T['pe0']=1-np.sqrt(T.phi0); X['pe0']=dict(med=float(T.pe0.median()),med_obs=float(T.pe.median()),corr=float(np.corrcoef(T.pe0,T.pe)[0,1]),med_gap=float((T.pe-T.pe0).median()),share_within5=float(((T.pe-T.pe0).abs()<0.05).mean()))
R['x']=X; json.dump(R,open('f/out/results.json','w'),ensure_ascii=False,indent=1)
for k,v in X.items(): print(k,json.dumps(v,ensure_ascii=False)[:1300]); print()
'''

# ==============================================================================
# f/f_fig.py
# ==============================================================================
MODULES['f/f_fig.py'] = r'''
# -*- coding: utf-8 -*-
import warnings; warnings.filterwarnings('ignore')
import json, matplotlib; matplotlib.use('Agg'); import matplotlib.pyplot as plt
import numpy as np, pandas as pd
from scipy import stats
plt.rcParams.update({'font.family':'DejaVu Sans','font.size':9,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':300})
R=json.load(open('f/out/results.json')); X=R['x']; D=pd.read_csv('f/out/country_window.csv',dtype={'un':str}); T=D[D.level=='taxa']; G=D[D.level=='group']; m=pd.read_parquet('f/out/marine.parquet')
B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'; Gn='#1e8449'; WN=['1950–1974','1975–1999','2000–2024']; WC={WN[0]:'#5dade2',WN[1]:O,WN[2]:Rd}
# S1 dunya avi
tot=m.groupby('yr').val.sum()/1e6; anc=(m[m.sp=='VET'].groupby('yr').val.sum()/1e6).reindex(tot.index).fillna(0)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
a=ax[0]; a.plot(tot.index,tot,color=B,lw=1.4,label='Toplam'); a.plot(tot.index,tot-anc,color=O,lw=1.2,label='Perulu hamsi dışında'); a.fill_between(tot.index,tot-anc,tot,color='#aed6f1',alpha=.6,lw=0)
a.set_ylabel('Deniz avı (milyon ton)'); a.legend(frameon=False,fontsize=7.5,loc='lower right'); a.set_title('(a) Dünya deniz avı, 1950–2024',fontsize=9,loc='left')
a=ax[1]; d1=tot.diff().dropna(); d2=(tot-anc).diff().dropna(); a.bar(d1.index,d1,color='#aed6f1',width=.9,label='Toplam'); a.plot(d2.index,d2,color=O,lw=1,label='Perulu hamsi dışında'); a.axhline(0,color='#999',lw=.6)
a.set_ylabel('Yıllık değişim (milyon ton)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Yıllık değişim',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('f/fig/Figure_1.png'); plt.close()
# S2 PE ve rho dagilimi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
a=ax[0]; a.hist(T.pe,bins=30,color='#aed6f1',edgecolor='white',lw=.3,label='Takson düzeyi'); a.hist(G.pe,bins=30,histtype='step',color=Rd,lw=1.1,label='ISSCAAP grubu düzeyi'); a.axvline(T.pe.median(),color=B,lw=1,ls='--'); a.axvline(G.pe.median(),color=Rd,lw=1,ls='--')
a.set_xlabel('Portföy etkisi (1 − √φ)'); a.set_ylabel('Ülke-dönem sayısı'); a.legend(frameon=False,fontsize=7.5,loc='upper left'); a.set_title('(a) Portföy etkisi',fontsize=9,loc='left')
a=ax[1]; a.hist(T.rho.clip(-0.3,0.6),bins=40,color='#aed6f1',edgecolor='white',lw=.3); a.axvline(0,color='#999',lw=.7); a.axvline(T.rho.median(),color=B,lw=1,ls='--'); a.set_xlabel('Taksonlar arası ortalama korelasyon'); a.set_ylabel('Ülke-dönem sayısı'); a.set_title('(b) Ortalama korelasyon',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('f/fig/Figure_2.png'); plt.close()
# S3 cesitlilik: es zamanlilik, agirlikli CV, toplam CV
fig,ax=plt.subplots(1,3,figsize=(7.4,2.8))
for a,col,t,yl in((ax[0],np.sqrt(T.phi),'(a) Eş zamanlılık','√φ'),(ax[1],T.cvw,'(b) Taksonların ağırlıklı değişkenliği','CV (ağırlıklı)'),(ax[2],T.cvT,'(c) Toplam avın değişkenliği','CV (toplam)')):
    for w in WN: s=T.win==w; a.loglog(T.neff[s],col[s],'.',ms=3,color=WC[w],alpha=.6,label=w)
    ok=(col>0); r=stats.linregress(np.log(T.neff[ok]),np.log(col[ok])); xx=np.array([1.2,80]); a.loglog(xx,np.exp(r.intercept+r.slope*np.log(xx)),color='#222',lw=1); a.text(0.05,0.06,f"eğim = {r.slope:.2f}".replace('.',',').replace('-','−'),transform=a.transAxes,fontsize=7.5)
    a.set_xlabel('Etkin takson sayısı'); a.set_ylabel(yl); a.set_title(t,fontsize=8.5,loc='left')
ax[0].legend(frameon=False,fontsize=6.5,loc='upper right',markerscale=2)
fig.tight_layout(); fig.savefig('f/fig/Figure_3.png'); plt.close()
# S4 donemlere gore (dengeli panel)
bal=T.groupby('un').win.nunique(); TB=T[T.un.isin(bal[bal==3].index)]; GB=G[G.un.isin(bal[bal==3].index)]
fig,ax=plt.subplots(1,3,figsize=(7.4,2.9))
def bx(a,d,col,pos,c):
    bp=a.boxplot([d[d.win==w][col].dropna() for w in WN],positions=np.arange(3)+pos,widths=.32,patch_artist=True,showfliers=False,medianprops=dict(color='#222'))
    for p_ in bp['boxes']: p_.set(facecolor=c,alpha=.7,lw=.5)
bx(ax[0],TB,'pe',-0.18,'#5dade2'); bx(ax[0],GB,'pe',0.18,'#f5b041'); ax[0].set_title('(a) Portföy etkisi',fontsize=9,loc='left')
bx(ax[1],TB,'cvT',0,'#5dade2'); ax[1].set_title('(b) Toplam avın değişkenliği',fontsize=9,loc='left')
bx(ax[2],TB,'neff',-0.18,'#5dade2'); bx(ax[2],GB,'neff',0.18,'#f5b041'); ax[2].set_yscale('log'); ax[2].set_title('(c) Etkin sayı',fontsize=9,loc='left')
for a in ax: a.set_xticks(range(3)); a.set_xticklabels(WN,fontsize=7)
import matplotlib.patches as mp; ax[0].legend(handles=[mp.Patch(color='#5dade2',alpha=.7,label='Takson'),mp.Patch(color='#f5b041',alpha=.7,label='ISSCAAP grubu')],frameon=False,fontsize=7,loc='lower right')
fig.tight_layout(); fig.savefig('f/fig/Figure_4.png'); plt.close()
# S5 buyukluk olceklenmesi + sok sikligi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
a=ax[0]
for w in WN: s=T.win==w; a.loglog(T['mean'][s],T.cvT[s],'.',ms=3.5,color=WC[w],alpha=.6,label=w)
r=stats.linregress(np.log(T['mean']),np.log(T.cvT)); xx=np.array([1e3,2e7]); a.loglog(xx,np.exp(r.intercept+r.slope*np.log(xx)),color='#222',lw=1,label=f"eğim = {r.slope:.2f}".replace('.',',').replace('-','−'))
a.set_xlabel('Ortalama yıllık av (ton)'); a.set_ylabel('Toplam avın değişkenliği (CV)'); a.legend(frameon=False,fontsize=7,markerscale=2); a.set_title('(a) Değişkenlik ve büyüklük',fontsize=9,loc='left')
a=ax[1]; q=X['shock']['by_q']; xs=np.arange(4); a.bar(xs-0.2,[100*v['rate10'] for v in q],width=.38,color='#f5b041',label="%10'dan fazla düşüş"); a.bar(xs+0.2,[100*v['rate20'] for v in q],width=.38,color=Rd,label="%20'den fazla düşüş")
a.set_xticks(xs); a.set_xticklabels([f"{v['neff']:.1f}".replace('.',',') for v in q]); a.set_xlabel('Etkin takson sayısı (dörtte birlik dilim ortancası)'); a.set_ylabel('Ülke-yılların payı (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Büyük düşüşlerin sıklığı',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('f/fig/Figure_5.png'); plt.close()
'''

if __name__ == "__main__":
    for _d in ('f', 'f/out', 'f/fig'): os.makedirs(_d, exist_ok=True)
    for _m in ORDER:
        with open(_m, 'w', encoding='utf-8') as _f: _f.write(MODULES[_m].lstrip('\n'))
    for _m in ORDER:
        if _m in SKIP_IF_EXISTS and os.path.exists(SKIP_IF_EXISTS[_m]):
            print('== atlandi / skipped (cikti mevcut / output present):', _m); continue
        print('==', _m, flush=True); subprocess.run([sys.executable, _m], check=True)   # her modul ayri surecte calisir / each module runs in its own process
    print('Bitti / done. Sonuclar / results: f/out/results.json')