# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Limanlar firmalar gibi mi büyür? 2.065 dünya limanında Gibrat yasası, büyüme oynaklığının ölçeklenmesi ve sıra hareketliliği, 2019–2025
# Do ports grow like firms? Gibrat's law, the scaling of growth volatility and rank mobility at 2,065 world ports, 2019–2025
#
# GOMD Yayınlar · Makale 18 — https://gomd.pirireis.edu.tr/publications.php?id=rp18
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `g/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `g/out/` klasörlerinden okur, sonuçları `g/out/` ve şekilleri `g/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `g/` and runs them in order, each in its own process;
# it reads data from `data/` and `g/out/`, and writes results to `g/out/` and figures to `g/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow
# Çalıştırma / Run:  python rp18_ports_gibrat_law_growth_volatility.py
#
# Port data: IMF PortWatch / UN Global Platform / University of Oxford (https://portwatch.imf.org), daily port activity and trade estimates derived from AIS; retrieved from the official open ArcGIS API on 13 November 2026.
# 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 = ['g/g_prep.py', 'g/g_analysis.py', 'g/g_extra.py', 'g/g_fig.py']
SKIP_IF_EXISTS = {}   # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package

# ==============================================================================
# g/g_prep.py
# ==============================================================================
MODULES['g/g_prep.py'] = r'''
import pandas as pd, numpy as np
d=pd.read_parquet('data/pwd/daily.parquet'); d['thr']=d['import']+d['export']; d['cthr']=d.import_container+d.export_container
d['yr']=d.date.dt.year; d['doy']=d.date.dt.dayofyear; d['par']=(d.date.dt.dayofyear%2)   # tek/cift gun
d['q']=d.date.dt.to_period('Q').astype(str); d['m']=d.date.dt.to_period('M').astype(str)
A=d.groupby(['portid','yr']).agg(calls=('portcalls','sum'),thr=('thr','sum'),imp=('import','sum'),exp=('export','sum'),cthr=('cthr','sum'),ccont=('portcalls_container','sum'),cbulk=('portcalls_dry_bulk','sum'),ctank=('portcalls_tanker','sum'),cgen=('portcalls_general_cargo','sum'),croro=('portcalls_roro','sum'),nd=('thr','size')).reset_index()
A.to_parquet('g/out/annual.parquet')
S=d.groupby(['portid','yr','par']).agg(calls=('portcalls','sum'),thr=('thr','sum')).reset_index(); S.to_parquet('g/out/annual_split.parquet')
# Oca-Eki (10 ay) karsilastirmasi 2026 icin
J=d[d.date.dt.month<=10].groupby(['portid','yr']).agg(calls=('portcalls','sum'),thr=('thr','sum')).reset_index(); J.to_parquet('g/out/janoct.parquet')
Mo=d.groupby(['portid','m']).agg(calls=('portcalls','sum'),thr=('thr','sum')).reset_index(); Mo.to_parquet('g/out/monthly.parquet')
print(A.groupby('yr')[['calls','thr']].agg(['sum',lambda x:(x>0).sum()]))
'''

# ==============================================================================
# g/g_analysis.py
# ==============================================================================
MODULES['g/g_analysis.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm, statsmodels.formula.api as smf
from scipy import stats, optimize
A=pd.read_parquet('g/out/annual.parquet'); S=pd.read_parquet('g/out/annual_split.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid')
WB=pd.read_csv('data/wb_country.csv') if True else None
R={}
FULL=list(range(2019,2026))
W={v:A[A.yr.isin(FULL)].pivot(index='portid',columns='yr',values=v) for v in('thr','calls','cthr','ccont','cbulk','ctank','cgen','croro')}
R['n_ports']=int(len(W['thr'])); R['n_days']=int(A.nd.sum())
R['world']={int(y):dict(thr=float(W['thr'][y].sum()),calls=float(W['calls'][y].sum()),n_pos=int((W['thr'][y]>0).sum())) for y in FULL}
# ornek: 2019-2025 her yil en az 12 ugrak ve pozitif yuk
def sample(v,minc=12):
    ok=(W['calls'][FULL].min(axis=1)>=minc)&(W[v][FULL].min(axis=1)>0); return W[v][ok]
T=sample('thr'); Cc=sample('calls'); R['n_sample']=int(len(T)); R['sample_share_thr']=float(T[2019].sum()/W['thr'][2019].sum())
def panel(X):
    L=np.log(X); rows=[]
    for y in FULL[1:]:
        rows.append(pd.DataFrame({'portid':X.index,'yr':y,'g':(L[y]-L[y-1]).values,'s':L[y-1].values,'s2':(L[y-2].values if y-2 in L else np.nan)}))
    return pd.concat(rows)
PT=panel(T); PC=panel(Cc)
PT['iso']=PT.portid.map(P.ISO3); PC['iso']=PC.portid.map(P.ISO3)
R['desc']={}
for nm,X,PP in(('thr',T,PT),('calls',Cc,PC)):
    g=PP.g; R['desc'][nm]=dict(n=int(len(g)),mean=float(g.mean()),sd=float(g.std()),med=float(g.median()),iqr=float(g.quantile(.75)-g.quantile(.25)),kurt=float(stats.kurtosis(g,fisher=False)),skew=float(stats.skew(g)),
        size_med=float(X[2019].median()),size_p10=float(X[2019].quantile(.1)),size_p90=float(X[2019].quantile(.9)),size_max=float(X[2019].max()))
# ---- 1. Gibrat regresyonlari
def ols(PP,f,cl='portid'):
    m=smf.ols(f,PP).fit(cov_type='cluster',cov_kwds={'groups':PP[cl]}); return m
R['gibrat']={}
for nm,PP in(('thr',PT),('calls',PC)):
    m=ols(PP,'g ~ s + C(yr)'); R['gibrat'][nm+'_ols']=dict(b=float(m.params['s']),se=float(m.bse['s']),p=float(m.pvalues['s']),n=int(m.nobs))
    # IV: s_{t-1} icin arac s_{t-2}
    d=PP.dropna(subset=['s2']).copy()
    fs=smf.ols('s ~ s2 + C(yr)',d).fit(); d['shat']=fs.fittedvalues
    from linearmodels.iv import IV2SLS
    iv=IV2SLS.from_formula('g ~ 1 + C(yr) + [s ~ s2]',d).fit(cov_type='clustered',clusters=d.portid)
    R['gibrat'][nm+'_iv']=dict(b=float(iv.params['s']),se=float(iv.std_errors['s']),p=float(iv.pvalues['s']),n=int(iv.nobs),fs_b=float(fs.params['s2']))
    m=ols(d,'g ~ s + C(yr)'); R['gibrat'][nm+'_ols_same']=dict(b=float(m.params['s']),se=float(m.bse['s']))
# bolunmus ornek: buyukluk tek gunlerden, buyume cift gunlerden
def split(v,ids):
    a=S[S.par==1].pivot(index='portid',columns='yr',values=v).loc[ids,FULL]; b=S[S.par==0].pivot(index='portid',columns='yr',values=v).loc[ids,FULL]
    ok=(a.min(axis=1)>0)&(b.min(axis=1)>0); a,b=np.log(a[ok]),np.log(b[ok]); rows=[]
    for y in FULL[1:]:
        rows.append(pd.DataFrame({'portid':a.index,'yr':y,'g_b':(b[y]-b[y-1]).values,'s_a':a[y-1].values,'g_a':(a[y]-a[y-1]).values,'s_b':b[y-1].values}))
    return pd.concat(rows)
for nm,X in(('thr',T),('calls',Cc)):
    d=split(nm,X.index)
    m=ols(d,'g_b ~ s_a + C(yr)'); m2=ols(d,'g_a ~ s_b + C(yr)'); m3=ols(d,'g_a ~ s_a + C(yr)')
    R['gibrat'][nm+'_split']=dict(b=float(m.params['s_a']),se=float(m.bse['s_a']),p=float(m.pvalues['s_a']),b2=float(m2.params['s_b']),se2=float(m2.bse['s_b']),b_same=float(m3.params['s_a']),se_same=float(m3.bse['s_a']),n=int(m.nobs),nports=int(d.portid.nunique()))
# uzun fark 2019->2025, baslangic buyuklugu 2019; IV ve split
for nm,X in(('thr',T),('calls',Cc)):
    d=pd.DataFrame({'g':np.log(X[2025]/X[2019]),'s':np.log(X[2019])}); m=smf.ols('g~s',d).fit(cov_type='HC1')
    R['gibrat'][nm+'_long']=dict(b=float(m.params['s']),se=float(m.bse['s']),p=float(m.pvalues['s']),n=int(m.nobs))
    # ortalama uclar: (2024+2025)/(2019+2020)
    d=pd.DataFrame({'g':np.log((X[2024]+X[2025])/(X[2019]+X[2020])),'s':np.log((X[2019]+X[2020])/2)})
# yillara gore beta (split)
d=split('thr',T.index); R['gibrat']['thr_split_year']={int(y):dict(b=float(m.params['s_a']),se=float(m.bse['s_a'])) for y in FULL[1:] for m in[smf.ols('g_b ~ s_a',d[d.yr==y]).fit(cov_type='HC1')]}
# buyukluk dilimleri: ortalama buyume (split ile)
d['bin']=d.groupby('yr').s_a.transform(lambda x: pd.qcut(x,10,labels=False))
R['gibrat']['thr_decile']=[dict(dec=int(k)+1,size=float(np.exp(x.s_a.median())),g_mean=float(x.g_b.mean()),g_med=float(x.g_b.median()),sd=float(x.g_b.std()),iqr=float(x.g_b.quantile(.75)-x.g_b.quantile(.25)),n=int(len(x))) for k,x in d.groupby('bin')]
dc=split('calls',Cc.index); dc['bin']=dc.groupby('yr').s_a.transform(lambda x: pd.qcut(x,10,labels=False))
R['gibrat']['calls_decile']=[dict(dec=int(k)+1,size=float(np.exp(x.s_a.median())*2),g_mean=float(x.g_b.mean()),g_med=float(x.g_b.median()),sd=float(x.g_b.std()),iqr=float(x.g_b.quantile(.75)-x.g_b.quantile(.25)),n=int(len(x))) for k,x in dc.groupby('bin')]
# ---- 2. Oynaklik olceklenmesi: sigma(g|S) ~ S^-beta
R['scaling']={}
def scal(PP,nb=20):
    PP=PP.copy(); m=smf.ols('g ~ s + C(yr)',PP).fit(); PP['e']=m.resid
    PP['bin']=pd.qcut(PP.s,nb,labels=False); B=PP.groupby('bin').agg(s=('s','median'),sd=('e','std'),mad=('e',lambda x:1.4826*np.median(np.abs(x-np.median(x)))),n=('e','size'))
    r1=stats.linregress(B.s,np.log(B.sd)); r2=stats.linregress(B.s,np.log(B.mad))
    # istasyon duzeyi: limanin kendi buyume sd'si ~ ortalama log buyukluk
    q=PP.groupby('portid').agg(s=('s','mean'),sd=('g','std')); q=q[q.sd>0]; r3=sm.OLS(np.log(q.sd),sm.add_constant(q.s)).fit(cov_type='HC1')
    # ML: e ~ N(0, exp(a - b s)^2)
    def nll(p): sg=np.exp(p[0]-p[1]*PP.s.values); return np.sum(np.log(sg)+0.5*(PP.e.values/sg)**2)
    o=optimize.minimize(nll,[0,0.2],method='Nelder-Mead')
    # Laplace ML: e ~ Laplace(0, exp(a - b s))
    def nll2(p): sg=np.exp(p[0]-p[1]*PP.s.values); return np.sum(np.log(2*sg)+np.abs(PP.e.values)/sg)
    o2=optimize.minimize(nll2,[0,0.2],method='Nelder-Mead')
    return dict(beta_sd=float(-r1.slope),se_sd=float(r1.stderr),r2_sd=float(r1.rvalue**2),beta_mad=float(-r2.slope),se_mad=float(r2.stderr),beta_port=float(-r3.params['s']),se_port=float(r3.bse['s']),beta_ml=float(o.x[1]),beta_lap=float(o2.x[1]),bins=[[float(a),float(b),float(c)] for a,b,c in zip(B.s,B.sd,B.mad)]),PP
R['scaling']['thr'],PTe=scal(PT); R['scaling']['calls'],PCe=scal(PC)
# yuk oynakligi ugrak sayisina gore (bilesik surec): thr buyumesi sd ~ calls^-beta
PT2=PT.merge(PC[['portid','yr','s']].rename(columns={'s':'sc'}),on=['portid','yr'],how='inner'); m=smf.ols('g ~ s + C(yr)',PT2).fit(); PT2['e']=m.resid
PT2['bin']=pd.qcut(PT2.sc,20,labels=False); B=PT2.groupby('bin').agg(s=('sc','median'),sd=('e','std')); r=stats.linregress(B.s,np.log(B.sd)); R['scaling']['thr_by_calls']=dict(beta=float(-r.slope),se=float(r.stderr),r2=float(r.rvalue**2),bins=[[float(a),float(b)] for a,b in zip(B.s,B.sd)])
# bootstrap (liman bloklari) beta_sd icin
rng=np.random.default_rng(3)
def bsb(PP,nb=20,B=400):
    ids=PP.portid.unique(); grp={k:v for k,v in PP.groupby('portid')}; out=[]
    for _ in range(B):
        s=pd.concat([grp[i] for i in rng.choice(ids,len(ids))]); s['bin']=pd.qcut(s.s,nb,labels=False,duplicates='drop'); q=s.groupby('bin').agg(s=('s','median'),sd=('e','std')); out.append(-stats.linregress(q.s,np.log(q.sd)).slope)
    return [float(np.percentile(out,2.5)),float(np.percentile(out,97.5))]
R['scaling']['thr']['ci']=bsb(PTe); R['scaling']['calls']['ci']=bsb(PCe)
# ---- 3. Buyume orani dagilimi: Laplace vs Gauss (olceklenmis artiklar)
R['dist']={}
for nm,PP,b in(('thr',PTe,R['scaling']['thr']['beta_sd']),('calls',PCe,R['scaling']['calls']['beta_sd'])):
    PP['bin']=pd.qcut(PP.s,20,labels=False); z=PP.e/PP.groupby('bin').e.transform('std'); z=z-z.mean()
    lg=stats.norm.logpdf(z,*stats.norm.fit(z)).sum(); ll=stats.laplace.logpdf(z,*stats.laplace.fit(z)).sum()
    # Subbotin (gennorm) sekil parametresi
    bta,loc,sc=stats.gennorm.fit(z,floc=0); lsub=stats.gennorm.logpdf(z,bta,loc,sc).sum()
    R['dist'][nm]=dict(n=int(len(z)),aic_norm=float(-2*lg+4),aic_lap=float(-2*ll+4),aic_sub=float(-2*lsub+4),subbotin_b=float(bta),kurt=float(stats.kurtosis(z,fisher=False)),p_gt3=float((np.abs(z)>3).mean()),hist=[[float(a),float(b)] for a,b in zip(*[(lambda h,e:(0.5*(e[1:]+e[:-1]),h))(*np.histogram(z,bins=np.linspace(-6,6,61),density=True))][0])])
# ---- 4. Buyukluk dagilimi: lognormal, Pareto kuyrugu
R['size']={}
for y in(2019,2022,2025):
    x=W['thr'][y]; x=x[x>0].sort_values(ascending=False).values
    lx=np.log(x); out=dict(n=int(len(x)),mu=float(lx.mean()),sig=float(lx.std()),skew=float(stats.skew(lx)),kurt=float(stats.kurtosis(lx,fisher=False)),gini=float((2*np.arange(1,len(x)+1)[::-1]*np.sort(x)[::-1]).sum()/(len(x)*x.sum())-(len(x)+1)/len(x)) )
    xs=np.sort(x); n_=len(xs); out['gini']=float((2*(np.arange(1,n_+1)*xs).sum())/(n_*xs.sum())-(n_+1)/n_)
    out['top10']=float(x[:10].sum()/x.sum()); out['top50']=float(x[:50].sum()/x.sum()); out['top100']=float(x[:100].sum()/x.sum()); out['hhi']=float(((x/x.sum())**2).sum()*1e4)
    for k in(50,100,200,400):
        t=x[:k]; out['hill_%d'%k]=float(1/np.mean(np.log(t/t[-1])[:-1])); 
        r=stats.linregress(np.log(t),np.log(np.arange(1,k+1)-0.5)); out['rank_%d'%k]=float(-r.slope)
    R['size'][str(y)]=out
x=W['thr'][2025]; x=x[x>0].values
# Clauset tarzi: xmin taramasi ile KS en kucuk
xs=np.sort(x)[::-1]; best=None
for k in range(30,800,10):
    t=xs[:k]; xm=t[-1]; a=1+ (k)/np.sum(np.log(t/xm)); cdf=1-(np.sort(t)/xm)**(-(a-1)); ks=np.max(np.abs(cdf-np.arange(1,k+1)/k))
    if best is None or ks<best[0]: best=(ks,k,a,xm)
# lognormal vs pareto kuyruk karsilastirmasi (Vuong benzeri LR) xmin ustunde
ks,k,a,xm=best; t=xs[:k]; llp=np.sum(np.log((a-1)/xm)-a*np.log(t/xm))
mu,sg=np.log(x).mean(),np.log(x).std()
# kesilmis lognormal
def nll(p):
    m,s=p; 
    if s<=0: return 1e18
    return -np.sum(stats.lognorm.logpdf(t,s,scale=np.exp(m))-np.log(stats.lognorm.sf(xm,s,scale=np.exp(m))))
o=optimize.minimize(nll,[mu,sg],method='Nelder-Mead'); lll=-o.fun
li=(np.log((a-1)/xm)-a*np.log(t/xm))-(stats.lognorm.logpdf(t,o.x[1],scale=np.exp(o.x[0]))-np.log(stats.lognorm.sf(xm,o.x[1],scale=np.exp(o.x[0]))))
v=li.sum()/(np.sqrt(k)*li.std()); R['size']['tail2025']=dict(ks=float(ks),k=int(k),alpha_ccdf=float(a-1),xmin=float(xm),ll_pareto=float(llp),ll_lognorm=float(lll),vuong=float(v),p=float(2*(1-stats.norm.cdf(abs(v)))))
sw=stats.shapiro(np.log(x)); jb=stats.jarque_bera(np.log(x)); R['size']['lognorm_test2025']=dict(shapiro_w=float(sw[0]),shapiro_p=float(sw[1]),jb=float(jb[0]),jb_p=float(jb[1]))
# ---- 5. Sira hareketliligi
X=W['thr'].loc[T.index]; rk=X.rank(ascending=False)
R['rank']=dict(spearman_19_25=float(stats.spearmanr(X[2019],X[2025])[0]),spearman_yoy={int(y):float(stats.spearmanr(X[y-1],X[y])[0]) for y in FULL[1:]},pearson_log_19_25=float(np.corrcoef(np.log(X[2019]),np.log(X[2025]))[0,1]))
for k in(10,20,50,100):
    a=set(rk.index[rk[2019]<=k]); b=set(rk.index[rk[2025]<=k]); R['rank']['stay_top%d'%k]=float(len(a&b)/k)
q19=pd.qcut(X[2019],5,labels=False); q25=pd.qcut(X[2025],5,labels=False); tm=pd.crosstab(q19,q25,normalize='index'); R['rank']['trans']=tm.values.tolist(); R['rank']['stay_diag']=float(np.mean(np.diag(tm.values)))
R['rank']['abs_rank_change_med']=float((rk[2025]-rk[2019]).abs().median()); 
# sira degisimi buyukluge gore: |dlog rank| 
dr=(np.log(rk[2025])-np.log(rk[2019])).abs(); R['rank']['dlogrank_by_q']=[float(dr[q19==i].median()) for i in range(5)]
top=X.sort_values(2025,ascending=False).head(15); R['rank']['top15']=[[P.portname[i],P.country[i],float(top.loc[i,2019]),float(top.loc[i,2025]),int(rk.loc[i,2019]),int(rk.loc[i,2025])] for i in top.index]
# ---- 6. Buyumenin kaliciligi (seri korelasyon) ve yas/tur farklari
d=PT.sort_values(['portid','yr']); d['gl']=d.groupby('portid').g.shift(1); dd=d.dropna(subset=['gl'])
m=ols(dd,'g ~ gl + s + C(yr)'); R['persist']=dict(rho=float(m.params['gl']),se=float(m.bse['gl']),p=float(m.pvalues['gl']),n=int(m.nobs))
dS=split('thr',T.index).sort_values(['portid','yr']); dS['gl_a']=dS.groupby('portid').g_a.shift(1); q=dS.dropna(subset=['gl_a']); m=ols(q,'g_b ~ gl_a + s_a + C(yr)'); R['persist']['rho_split']=float(m.params['gl_a']); R['persist']['se_split']=float(m.bse['gl_a'])
# liman turu: baskin gemi turu (2019 ugraklari)
ty=pd.DataFrame({k:W[v][2019] for k,v in(('Konteyner','ccont'),('Kuru dökme','cbulk'),('Tanker','ctank'),('Genel kargo','cgen'),('Ro-Ro','croro'))}); dom=ty.idxmax(axis=1); shr=ty.max(axis=1)/ty.sum(axis=1)
dS['type']=dS.portid.map(dom); dS['cont']=dS.portid.map(P.continent); R['hetero']={}
for k,x in dS.groupby('type'):
    m=ols(x,'g_b ~ s_a + C(yr)'); x=x.copy(); x['e']=m.resid; x['bin']=pd.qcut(x.s_a,10,labels=False); B=x.groupby('bin').agg(s=('s_a','median'),sd=('e','std')); r=stats.linregress(B.s,np.log(B.sd))
    R['hetero'][k]=dict(n=int(x.portid.nunique()),b=float(m.params['s_a']),se=float(m.bse['s_a']),beta=float(-r.slope),beta_se=float(r.stderr),g_sd=float(x.g_b.std()))
R['hetero_region']={}
for k,x in dS.groupby('cont'):
    m=ols(x,'g_b ~ s_a + C(yr)'); x=x.copy(); x['e']=m.resid; x['bin']=pd.qcut(x.s_a,10,labels=False); B=x.groupby('bin').agg(s=('s_a','median'),sd=('e','std')); r=stats.linregress(B.s,np.log(B.sd))
    R['hetero_region'][k]=dict(n=int(x.portid.nunique()),b=float(m.params['s_a']),se=float(m.bse['s_a']),beta=float(-r.slope),beta_se=float(r.stderr))
# ---- 7. Ulke toplami: taneciklilik (Gabaix): ulke buyume oynakligi ~ liman Herfindahl'i
Xc=X.copy(); Xc['iso']=Xc.index.map(P.ISO3); cs=Xc.groupby('iso')[FULL].sum(); h=Xc.groupby('iso').apply(lambda d:((d[2019]/d[2019].sum())**2).sum()); npc=Xc.groupby('iso').size()
gc=np.log(cs).diff(axis=1).iloc[:,1:]; q=pd.DataFrame({'sd':gc.std(axis=1),'h':h,'n':npc,'s':np.log(cs[2019])}); q=q[(q.n>=3)&(q.sd>0)]
m=smf.ols('np.log(sd) ~ np.log(h) + s',q).fit(cov_type='HC1'); R['granular']=dict(n=int(len(q)),b_h=float(m.params['np.log(h)']),se_h=float(m.bse['np.log(h)']),p_h=float(m.pvalues['np.log(h)']),b_s=float(m.params['s']),p_s=float(m.pvalues['s']),r2=float(m.rsquared))
# kuresel toplam buyumenin taneciklilik payi: en buyuk 100 limanin ozgun soklari
gw=np.log(X.sum()).diff().dropna(); gi=np.log(X).diff(axis=1).iloc[:,1:]; w=(X.shift(1,axis=1).iloc[:,1:]/X.shift(1,axis=1).iloc[:,1:].sum())
idio=gi-gi.median(); top100=rk[2019]<=100; G=(w[top100]*idio[top100]).sum(); R['granular']['world_g']={int(k):float(v) for k,v in gw.items()}; R['granular']['granular_resid']={int(k):float(v) for k,v in G.items()}; R['granular']['corr']=float(np.corrcoef(gw.values,G.values)[0,1])
# ---- 8. Saglamlik: esik 52 ugrak; 2026 Oca-Eki
T52=sample('thr',52); d=split('thr',T52.index); m=ols(d,'g_b ~ s_a + C(yr)'); x=panel(T52); mm=smf.ols('g ~ s + C(yr)',x).fit(); x['e']=mm.resid; x['bin']=pd.qcut(x.s,20,labels=False); B=x.groupby('bin').agg(s=('s','median'),sd=('e','std')); r=stats.linregress(B.s,np.log(B.sd))
R['robust']=dict(n52=int(len(T52)),b_split52=float(m.params['s_a']),se52=float(m.bse['s_a']),beta52=float(-r.slope),beta52_se=float(r.stderr))
J=pd.read_parquet('g/out/janoct.parquet').pivot(index='portid',columns='yr',values='thr').loc[T.index]; ok=(J[2025]>0)&(J[2026]>0); d=pd.DataFrame({'g':np.log(J[2026]/J[2025])[ok],'s':np.log(T[2024])[ok]}); m=smf.ols('g~s',d).fit(cov_type='HC1'); R['robust']['b_2026']=float(m.params['s']); R['robust']['se_2026']=float(m.bse['s']); R['robust']['g2026_world']=float(np.log(J[2026].sum()/J[2025].sum()))
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('g/out/results.json','w'),ensure_ascii=False,indent=1)
PTe.to_parquet('g/out/panel_thr.parquet'); PCe.to_parquet('g/out/panel_calls.parquet')
for k in R:
    s=json.dumps(clean(R[k]),ensure_ascii=False)
    print(k, s[:1800] if k not in('dist',) else {a:{x:y for x,y in b.items() if x!='hist'} for a,b in R[k].items()}); print()
'''

# ==============================================================================
# g/g_extra.py
# ==============================================================================
MODULES['g/g_extra.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.formula.api as smf
from scipy import stats
R=json.load(open('g/out/results.json')); A=pd.read_parquet('g/out/annual.parquet'); Mo=pd.read_parquet('g/out/monthly.parquet'); PT=pd.read_parquet('g/out/panel_thr.parquet')
ids=PT.portid.unique(); X={}
# (a) ufuk: aylik (yildan yila ayni ay), ceyreklik (yildan yila), yillik
Mo=Mo[Mo.portid.isin(ids)].copy(); Mo['per']=pd.PeriodIndex(Mo.m,freq='M'); Mo=Mo[Mo.per<='2025-12']
def horizon(df,key,lag,nm):
    w=df.pivot(index='portid',columns=key,values='thr'); w=w.where(w>0); L=np.log(w); g=(L-L.shift(lag,axis=1)).stack().rename('g'); s=L.shift(lag,axis=1).stack().rename('s'); d=pd.concat([g,s],axis=1).dropna().reset_index()
    d['t']=d[key].astype(str); d['e']=d.g-d.groupby('t').g.transform('mean'); m=smf.ols('e ~ s',d).fit(); d['e']=m.resid
    d['bin']=pd.qcut(d.s,20,labels=False); B=d.groupby('bin').agg(s=('s','median'),sd=('e','std')); r=stats.linregress(B.s,np.log(B.sd))
    X[nm]=dict(n=int(len(d)),sd=float(d.g.std()),beta=float(-r.slope),se=float(r.stderr),r2=float(r.rvalue**2))
horizon(Mo,'per',12,'monthly_yoy')
Mo['q']=Mo.per.dt.asfreq('Q'); Q=Mo.groupby(['portid','q']).thr.sum().reset_index(); horizon(Q,'q',4,'quarterly_yoy')
Mo['h']=Mo.per.dt.year.astype(str)+'H'+((Mo.per.dt.month-1)//6+1).astype(str); H=Mo.groupby(['portid','h']).thr.sum().reset_index(); horizon(H,'h',2,'half_yoy')
X['annual']=dict(beta=R['scaling']['thr']['beta_sd'],sd=R['desc']['thr']['sd'])
# (b) ithalat/ihracat ayri
for v in('imp','exp'):
    w=A[A.yr<=2025].pivot(index='portid',columns='yr',values=v).loc[ids]; w=w[(w>0).all(axis=1)]; L=np.log(w); rows=[]
    for y in range(2020,2026): rows.append(pd.DataFrame({'portid':w.index,'yr':y,'g':(L[y]-L[y-1]).values,'s':L[y-1].values,'s2':(L[y-2].values if y-2 in L else np.nan)}))
    d=pd.concat(rows); m=smf.ols('g ~ s + C(yr)',d).fit(); d['e']=m.resid; d['bin']=pd.qcut(d.s,20,labels=False); B=d.groupby('bin').agg(s=('s','median'),sd=('e','std')); r=stats.linregress(B.s,np.log(B.sd))
    from linearmodels.iv import IV2SLS
    dd=d.dropna(); iv=IV2SLS.from_formula('g ~ 1 + C(yr) + [s ~ s2]',dd).fit(cov_type='clustered',clusters=dd.portid)
    X[v]=dict(n=int(w.shape[0]),b_ols=float(m.params['s']),b_iv=float(iv.params['s']),se_iv=float(iv.std_errors['s']),beta=float(-r.slope),se=float(r.stderr),sd=float(d.g.std()))
# (c) cesitlenme: gemi turu entropisi (2019 ugraklari) ve yon dengesi; liman oynakligi ~ buyukluk + cesitlenme
W={k:A[A.yr==2019].set_index('portid')[k] for k in('ccont','cbulk','ctank','cgen','croro','imp','exp','calls','thr')}
ty=pd.DataFrame({k:W[k] for k in('ccont','cbulk','ctank','cgen','croro')}).loc[ids]; sh=ty.div(ty.sum(axis=1),axis=0)
ent=-(sh*np.log(sh.where(sh>0))).sum(axis=1); neff=np.exp(ent)
q=PT.groupby('portid').agg(sd=('g','std'),s=('s','mean')); q['neff']=neff; q['bal']=1-(W['imp'].loc[q.index]-W['exp'].loc[q.index]).abs()/(W['imp'].loc[q.index]+W['exp'].loc[q.index]); q['lc']=np.log(W['calls'].loc[q.index]); q=q[(q.sd>0)&q.neff.notna()]
q['avg']=q.s-q.lc   # log ugrak basina yuk
m=smf.ols('np.log(sd) ~ s + np.log(neff)',q).fit(cov_type='HC1'); m2=smf.ols('np.log(sd) ~ s + np.log(neff) + bal',q).fit(cov_type='HC1'); m3=smf.ols('np.log(sd) ~ lc + avg + np.log(neff) + bal',q).fit(cov_type='HC1')
X['divers']=dict(n=int(len(q)),neff_mean=float(q.neff.mean()),neff_p10=float(q.neff.quantile(.1)),neff_p90=float(q.neff.quantile(.9)),b_s=float(m.params['s']),b_neff=float(m.params['np.log(neff)']),se_neff=float(m.bse['np.log(neff)']),p_neff=float(m.pvalues['np.log(neff)']),r2=float(m.rsquared),
   b_s2=float(m2.params['s']),b_neff2=float(m2.params['np.log(neff)']),se_neff2=float(m2.bse['np.log(neff)']),b_bal=float(m2.params['bal']),se_bal=float(m2.bse['bal']),p_bal=float(m2.pvalues['bal']),r2_2=float(m2.rsquared),
   b_lc=float(m3.params['lc']),se_lc=float(m3.bse['lc']),b_avg=float(m3.params['avg']),se_avg=float(m3.bse['avg']),b_neff3=float(m3.params['np.log(neff)']),se_neff3=float(m3.bse['np.log(neff)']),r2_3=float(m3.rsquared),
   corr_s_neff=float(np.corrcoef(q.s,np.log(q.neff))[0,1]))
# etkisi: tek turden (neff=1) bes ture (neff~3) gecis
X['divers']['effect_1to3']=float(np.exp(m2.params['np.log(neff)']*np.log(3))-1)
# (d) en buyuk sicramalar
PT['name']=PT.portid.map(pd.read_csv('data/pw/ports_global.csv').set_index('portid').fullname)
PT=PT.reset_index(drop=True); big=PT[PT.s>np.log(2e7)].copy(); big['ag']=big.g.abs(); X['big_jumps']=dict(n=int(len(big)),p_abs20=float((big.ag>0.2).mean()),p_abs50=float((big.ag>0.5).mean()),top=[[str(r['name']),int(r['yr']),float(r['g'])] for _,r in big.sort_values('ag',ascending=False).head(8).iterrows()])
allp=PT; X['all_jumps']=dict(p_abs20=float((allp.g.abs()>0.2).mean()),p_abs50=float((allp.g.abs()>0.5).mean()),p_half=float((allp.g<np.log(0.5)).mean()),p_double=float((allp.g>np.log(2)).mean()))
R['x']=X; json.dump(R,open('g/out/results.json','w'),ensure_ascii=False,indent=1); print(json.dumps(X,ensure_ascii=False,indent=0)[:3500])
'''

# ==============================================================================
# g/g_fig.py
# ==============================================================================
MODULES['g/g_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('g/out/results.json')); A=pd.read_parquet('g/out/annual.parquet')
B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'
# S1 buyukluk dagilimi: log histogram + sira-buyukluk
x=A[(A.yr==2025)&(A.thr>0)].thr.values; lx=np.log10(x)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; a.hist(lx,bins=45,density=True,color='#aed6f1',edgecolor='white',lw=.3); g=np.linspace(lx.min(),lx.max(),200); a.plot(g,stats.norm.pdf(g,lx.mean(),lx.std()),color=B,lw=1.3,label='Lognormal uyum')
a.set_xlabel('Yıllık yük, 2025 (log₁₀ ton)'); a.set_ylabel('Yoğunluk'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Büyüklük dağılımı',fontsize=9,loc='left')
a=ax[1]; xs=np.sort(x)[::-1]; r=np.arange(1,len(xs)+1); a.loglog(xs,r,'.',ms=2,color=B); t=R['size']['tail2025']; k=t['k']; al=t['alpha_ccdf']; xx=np.array([t['xmin'],xs[0]]); a.loglog(xx,k*(xx/t['xmin'])**(-al),color=Rd,lw=1.2,label=f"Pareto kuyruğu, üs = {al:.2f}".replace('.',','))
m,s=np.log(x).mean(),np.log(x).std(); gg=np.logspace(np.log10(xs[-1]),np.log10(xs[0]),200); a.loglog(gg,len(xs)*stats.norm.sf((np.log(gg)-m)/s),color=O,lw=1,ls='--',label='Lognormal')
a.set_xlabel('Yıllık yük (ton)'); a.set_ylabel('Sıra'); a.legend(frameon=False,fontsize=7.5,loc='lower left'); a.set_title('(b) Sıra–büyüklük ilişkisi',fontsize=9,loc='left'); a.set_ylim(.8,3000)
fig.tight_layout(); fig.savefig('g/fig/Figure_1.png'); plt.close()
# S2 Gibrat: buyukluk dilimlerine gore ortalama buyume (split) + yillara gore beta
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; d=R['gibrat']['thr_decile']; a.errorbar([q['size'] for q in d],[100*q['g_mean'] for q in d],yerr=[196*q['sd']/np.sqrt(q['n']) for q in d],fmt='o-',color=B,ms=3.5,lw=1.1,capsize=2,label='Ortalama'); a.plot([q['size'] for q in d],[100*q['g_med'] for q in d],'s--',color=O,ms=3,lw=.9,label='Ortanca')
a.axhline(0,color='#999',lw=.6); a.set_xscale('log'); a.set_xlabel('Önceki yıl yükü (ton; yarı örneklem)'); a.set_ylabel('Yıllık log büyüme × 100'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Büyüklük onda birliklerine göre büyüme',fontsize=9,loc='left')
a=ax[1]; y=R['gibrat']['thr_split_year']; ys=sorted(y); a.errorbar([int(k) for k in ys],[y[k]['b'] for k in ys],yerr=[1.96*y[k]['se'] for k in ys],fmt='o',color=B,ms=4,capsize=2); a.axhline(0,color='#999',lw=.6); a.axhline(R['gibrat']['thr_split']['b'],color=Rd,lw=.9,ls='--',label='Havuzlanmış'); a.set_ylabel('Büyüklük katsayısı (β)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Yıllara göre büyüklük katsayısı',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('g/fig/Figure_2.png'); plt.close()
# S3 oynaklik olceklenmesi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
for a,k,t,xl in((ax[0],'thr','(a) Yük','Önceki yıl yükü (ton)'),(ax[1],'calls','(b) Gemi uğrağı','Önceki yıl uğrak sayısı')):
    b=np.array(R['scaling'][k]['bins']); a.loglog(np.exp(b[:,0]),b[:,1],'o',color=B,ms=4,label='Standart sapma'); a.loglog(np.exp(b[:,0]),b[:,2],'s',color=O,ms=3.5,label='Sağlam ölçü (MAD)')
    be=R['scaling'][k]['beta_sd']; r=stats.linregress(b[:,0],np.log(b[:,1])); a.loglog(np.exp(b[:,0]),np.exp(r.intercept+r.slope*b[:,0]),color=Rd,lw=1,label=f"β = {be:.2f}".replace('.',','))
    a.loglog(np.exp(b[:,0]),np.exp(np.log(b[0,1])-0.5*(b[:,0]-b[0,0])),color='#999',lw=.8,ls=':',label='Bağımsız birimler (β = 0,5)')
    a.set_xlabel(xl); a.set_ylabel('Büyüme oranının yayılımı'); a.set_title(t,fontsize=9,loc='left'); a.legend(frameon=False,fontsize=7)
fig.tight_layout(); fig.savefig('g/fig/Figure_3.png'); plt.close()
# S4 buyume dagilimi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
for a,k,t in((ax[0],'thr','(a) Yük'),(ax[1],'calls','(b) Gemi uğrağı')):
    h=np.array(R['dist'][k]['hist']); m=h[:,1]>0; a.semilogy(h[m,0],h[m,1],'o',color=B,ms=3); g=np.linspace(-6,6,300); a.semilogy(g,stats.laplace.pdf(g,0,1/np.sqrt(2)),color=Rd,lw=1.1,label='Laplace'); a.semilogy(g,stats.norm.pdf(g),color='#999',lw=1,ls='--',label='Normal')
    a.set_ylim(1e-4,2); a.set_xlabel('Ölçeklenmiş büyüme artığı'); a.set_ylabel('Yoğunluk'); a.set_title(t,fontsize=9,loc='left'); a.legend(frameon=False,fontsize=7.5)
fig.tight_layout(); fig.savefig('g/fig/Figure_4.png'); plt.close()
# S5 sira hareketliligi: 2019 vs 2025 sira; gecis matrisi
W=A[A.yr.isin([2019,2025])].pivot(index='portid',columns='yr',values='thr'); W=W[(W>0).all(axis=1)]; rk=W.rank(ascending=False)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.2))
a=ax[0]; a.loglog(rk[2019],rk[2025],'.',ms=2,color=B,alpha=.6); a.plot([1,2000],[1,2000],color='#999',lw=.7); a.set_xlabel('Sıra, 2019'); a.set_ylabel('Sıra, 2025'); a.set_title('(a) Yük sıralaması',fontsize=9,loc='left')
a=ax[1]; T=np.array(R['rank']['trans']); im=a.imshow(100*T,cmap='Blues',vmin=0,vmax=100)
for i in range(5):
    for j in range(5): a.text(j,i,f"{100*T[i,j]:.0f}",ha='center',va='center',fontsize=8,color='white' if T[i,j]>.5 else '#222')
a.set_xticks(range(5)); a.set_yticks(range(5)); a.set_xticklabels(['1','2','3','4','5']); a.set_yticklabels(['1','2','3','4','5']); a.set_xlabel('Beşte birlik dilim, 2025'); a.set_ylabel('Beşte birlik dilim, 2019'); a.set_title('(b) Geçiş matrisi (%)',fontsize=9,loc='left')
for s in a.spines.values(): s.set_visible(False)
fig.tight_layout(); fig.savefig('g/fig/Figure_5.png'); plt.close()
'''

if __name__ == "__main__":
    for _d in ('g', 'g/out', 'g/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: g/out/results.json')