# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Uğrak sayısı neyi ölçer? Konteyner limanlarında uğrak ve yük döngülerinin ayrışması ve dünya limanlarında ortak hareket, 2019–2026
# What do port calls measure? The divergence of call and cargo cycles at container ports and common movement across world ports, 2019–2026
#
# GOMD Yayınlar · Makale 22 — https://gomd.pirireis.edu.tr/publications.php?id=rp22
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `w/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `w/out/` klasörlerinden okur, sonuçları `w/out/` ve şekilleri `w/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `w/` and runs them in order, each in its own process;
# it reads data from `data/` and `w/out/`, and writes results to `w/out/` and figures to `w/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow xlrd
# Çalıştırma / Run:  python rp22_port_calls_vs_cargo_container_divergence.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. Global Supply Chain Pressure Index: Federal Reserve Bank of New York (https://www.newyorkfed.org/research/policy/gscpi).
# 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 = ['w/w_analysis.py', 'w/w_f2.py', 'w/w_extra.py', 'w/w_extra2.py', 'w/w_extra3.py', 'w/w_cont.py', 'w/w_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

# ==============================================================================
# w/w_analysis.py
# ==============================================================================
MODULES['w/w_analysis.py'] = r'''
# Dunya liman dongusu: kuresel / bolgesel / ulkesel / ozgun varyans ayristirmasi (PortWatch aylik)
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm
d=pd.read_parquet('data/pwd/daily.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid'); wb=pd.read_csv('data/wb_country.csv').set_index('iso3')
d=d[d.date<='2026-10-31'].copy(); d['thr']=d['import']+d['export']; d['m']=d.date.dt.to_period('M')
def panel(col):
    M=d.groupby(['portid','m'])[col].sum().unstack('portid'); return M
R={}
def dh(M):  # Davis-Haltiwanger yildan yila buyume
    a=M; b=M.shift(12); g=2*(a-b)/(a+b); return g.iloc[12:]
def run(col,minmean,label):
    M=panel(col); keep=(M.iloc[:12].sum()>0)&(M.mean()>=minmean)&((M==0).mean()<0.25); M=M.loc[:,keep]
    G=dh(M).fillna(0); T,N=G.shape; ids=G.columns
    w=M.mean()/M.mean().sum()
    Z=(G-G.mean())/G.std(ddof=0); Z=Z.fillna(0)
    # 1. kuresel faktor: standartlastirilmis serilerin 1. temel bileseni
    def pc(Zm,k=1):
        U,S,Vt=np.linalg.svd(Zm.values,full_matrices=False); F=U[:,:k]*np.sqrt(len(Zm)); ev=S**2/ (Zm.shape[0]*Zm.shape[1]); return F,ev
    F,ev=pc(Z,5); Fg=F[:,0]; 
    agg=(G*w).sum(axis=1); 
    if np.corrcoef(Fg,agg)[0,1]<0: Fg=-Fg
    # Marchenko-Pastur ust siniri (standart seri, sigma2=1): (1+sqrt(T/N))^2 / T  (ozdegerler toplam varyans payi cinsinden)
    q=T/N; mp=(1+np.sqrt(q))**2/T
    reg=P.continent.reindex(ids); iso=P.ISO3.reindex(ids)
    def resid(Y,Xm):
        Xc=sm.add_constant(Xm); B=np.linalg.lstsq(Xc,Y,rcond=None)[0]; return Y-Xc@B
    Zv=Z.values; E1=resid(Zv,Fg[:,None])
    # 2. bolgesel faktorler: kuresel faktorden arindirilmis artiklarin bolge ici 1. bileseni
    Freg={}; 
    for r in reg.dropna().unique():
        idx=np.where(reg.values==r)[0]
        if len(idx)<10: continue
        U,S,Vt=np.linalg.svd(E1[:,idx],full_matrices=False); f=U[:,0]*np.sqrt(T); 
        if np.corrcoef(f,E1[:,idx].mean(1))[0,1]<0: f=-f
        Freg[r]=f
    E2=E1.copy()
    for r,f in Freg.items():
        idx=np.where(reg.values==r)[0]; E2[:,idx]=resid(E1[:,idx],f[:,None])
    # 3. ulke faktorleri (>=8 liman)
    Fc={}; E3=E2.copy()
    for c in iso.dropna().unique():
        idx=np.where(iso.values==c)[0]
        if len(idx)<8: continue
        U,S,Vt=np.linalg.svd(E2[:,idx],full_matrices=False); f=U[:,0]*np.sqrt(T); Fc[c]=f; E3[:,idx]=resid(E2[:,idx],f[:,None])
    v0=Zv.var(0); sg=1-E1.var(0)/v0; sr=(E1.var(0)-E2.var(0))/v0; sc=(E2.var(0)-E3.var(0))/v0; si=E3.var(0)/v0
    S=pd.DataFrame({'glob':sg,'reg':sr,'cty':sc,'idio':si,'w':w.values,'mean':M.mean().values,'cont':reg.values,'iso':iso.values,'has_cty':[c in Fc for c in iso.values]},index=ids)
    S['load']=[np.corrcoef(Fg,Zv[:,j])[0,1] for j in range(N)]
    # 4. sifir dagilimi: her limanin serisini dairesel kaydir (oz iliskiyi korur, es zamanliligi bozar)
    rng=np.random.default_rng(11); null_g=[]; null_ev=[]
    for k in range(200):
        Zs=np.column_stack([np.roll(Zv[:,j],rng.integers(1,T)) for j in range(N)]); U,Sv,Vt=np.linalg.svd(Zs,full_matrices=False); f=U[:,0]*np.sqrt(T); null_ev.append(Sv[0]**2/(T*N)); null_g.append(np.median(1-resid(Zs,f[:,None]).var(0)/Zs.var(0)))
    out=dict(label=label,N=int(N),T=int(T),share_world=float(M.sum().sum()/panel(col).sum().sum()),ev=[float(x) for x in ev[:5]],mp=float(mp),null_ev1=float(np.mean(null_ev)),null_ev1_p99=float(np.percentile(null_ev,99)),null_glob_med=float(np.mean(null_g)),
        med={k:float(S[k].median()) for k in('glob','reg','cty','idio')},mean={k:float(S[k].mean()) for k in('glob','reg','cty','idio')},wmean={k:float(np.average(S[k],weights=S.w)) for k in('glob','reg','cty','idio')},
        p90_glob=float(S.glob.quantile(.9)),share_glob_gt10=float((S.glob>0.10).mean()),share_glob_gt25=float((S.glob>0.25).mean()),neg_load=float((S.load<0).mean()),n_cty=int(len(Fc)),n_reg=int(len(Freg)),
        n_with_cty=int(S.has_cty.sum()),cty_share_among=float(S[S.has_cty].cty.median()))
    out['by_cont']={str(r):dict(n=int(len(g)),glob=float(g.glob.median()),reg=float(g.reg.median()),cty=float(g.cty.median()),idio=float(g.idio.median()),wglob=float(np.average(g.glob,weights=g.w))) for r,g in S.groupby('cont')}
    S['q']=pd.qcut(S['mean'],5,labels=False); out['by_size']=[dict(q=int(k)+1,size=float(g['mean'].median()),glob=float(g.glob.median()),reg=float(g.reg.median()),cty=float(g.cty.median()),idio=float(g.idio.median())) for k,g in S.groupby('q')]
    # 5. toplam buyumenin ayristirilmasi: ortak bilesen ve tanecikli (ozgun) bilesen
    Gv=G.values; mu=Gv.mean(0); sd=Gv.std(0); Xf=sm.add_constant(Fg); B=np.linalg.lstsq(Xf,Gv,rcond=None)[0]; common=(Xf@B); idio=Gv-common
    aggc=(common*w.values).sum(1); aggi=(idio*w.values).sum(1); va=agg.var(ddof=0)
    out['agg']=dict(sd=float(agg.std()),share_common=float(np.var(aggc)/va),share_idio=float(np.var(aggi)/va),cov=float(2*np.cov(aggc,aggi,ddof=0)[0,1]/va),corr_F=float(np.corrcoef(Fg,agg)[0,1]),herf=float((w**2).sum()),top10w=float(w.sort_values(ascending=False).head(10).sum()),top100w=float(w.sort_values(ascending=False).head(100).sum()),
        expected_idio_sd=float(np.sqrt(((w.values*idio.std(0))**2).sum())),idio_sd=float(aggi.std()))
    # en buyuk 10 limanin ozgun katkisi
    top=w.sort_values(ascending=False).head(10).index; j=[list(ids).index(t) for t in top]; out['agg']['share_idio_top10']=float(np.var((idio[:,j]*w.values[j]).sum(1))/va)
    # 6. zaman icinde es zamanlilik: 36 aylik kayan pencere, ortalama ikili korelasyon ve PC1 payi
    roll=[]
    for s in range(0,T-36+1,3):
        Zw=Zv[s:s+36]; Zw=(Zw-Zw.mean(0))/np.where(Zw.std(0)>0,Zw.std(0),1); U,Sv,Vt=np.linalg.svd(Zw,full_matrices=False); e1=Sv[0]**2/(36*N); 
        # ortalama ikili korelasyon = (var(ortalama z)*N - 1)/(N-1)
        mz=Zw.mean(1); rbar=(mz.var()*N-1)/(N-1); roll.append([str(G.index[s+35]),float(e1),float(rbar)])
    out['roll']=roll
    mz=Z.values.mean(1); out['rbar_all']=float((mz.var()*N-1)/(N-1))
    # covid haric
    ok=np.array([(p_.year>2021) for p_ in G.index]); Zx=Zv[ok]; Zx=(Zx-Zx.mean(0))/np.where(Zx.std(0)>0,Zx.std(0),1); U,Sv,Vt=np.linalg.svd(Zx,full_matrices=False); fx=U[:,0]*np.sqrt(ok.sum()); out['post2021']=dict(T=int(ok.sum()),ev1=float(Sv[0]**2/(ok.sum()*N)),glob_med=float(np.median(1-resid(Zx,fx[:,None]).var(0))),rbar=float((Zx.mean(1).var()*N-1)/(N-1)),mp=float((1+np.sqrt(ok.sum()/N))**2/ok.sum()))
    S.to_csv(f'w/out/shares_{label}.csv'); pd.DataFrame({'m':G.index.astype(str),'Fg':Fg,'agg':agg.values,**{'reg_'+k:v for k,v in Freg.items()}}).to_csv(f'w/out/factors_{label}.csv',index=False)
    return out,S,G,Fg,Freg,agg
R['calls'],Sc,Gc,Fc_,Fregc,aggc=run('portcalls',10,'calls')
R['thr'],St,Gt,Ft,Fregt,aggt=run('thr',1000,'thr')
R['cont']=run('portcalls_container',5,'container')[0]; R['bulk']=run('portcalls_dry_bulk',5,'bulk')[0]; R['tank']=run('portcalls_tanker',5,'tanker')[0]
R['corr_F_calls_thr']=float(np.corrcoef(Fc_,Ft)[0,1])
# 7. belirleyiciler: kuresel pay ~ buyukluk, konteyner payi, gelir, kita
ty=d.groupby('portid')[['portcalls_container','portcalls_dry_bulk','portcalls_tanker','portcalls_general_cargo','portcalls_roro']].sum(); sh=ty.div(ty.sum(axis=1),axis=0)
import statsmodels.formula.api as smf
D=Sc.join(sh).copy(); D['inc']=D.iso.map(wb.income); D['ls']=np.log(D['mean']); D=D.dropna(subset=['inc','cont','portcalls_container']); D['g_']=pd.factorize(D.iso)[0]
m_=smf.ols('glob ~ ls + portcalls_container + portcalls_tanker + portcalls_dry_bulk + portcalls_roro + C(inc) + C(cont)',D).fit(cov_type='cluster',cov_kwds={'groups':D.g_})
R['det']={k:dict(b=float(m_.params[k]),se=float(m_.bse[k]),p=float(m_.pvalues[k])) for k in m_.params.index if k!='Intercept'}; R['det']['_r2']=float(m_.rsquared); R['det']['_n']=int(m_.nobs)
# 8. bolgesel donguler arasi oncul-ardil iliski: bolge toplam buyumeleri (ugrak) capraz korelasyon
REG=['Asia & Pacific','Europe','North America','South America','Africa']
ra={}
M=panel('portcalls'); Mt=panel('thr'); Mc=panel('portcalls_container'); Mi=d.groupby(['portid','m']).import_container.sum().unstack('portid'); Me=d.groupby(['portid','m']).export_container.sum().unstack('portid')
def regsum(Mx,r): 
    ids=[i for i in Mx.columns if P.continent.get(i)==r]; s=Mx[ids].sum(axis=1); return (np.log(s)-np.log(s.shift(12))).iloc[12:]
RA=pd.DataFrame({r:regsum(M,r) for r in REG}); RT=pd.DataFrame({r:regsum(Mt,r) for r in REG})
R['reg_corr']={a:{b:float(RA[a].corr(RA[b])) for b in REG} for a in REG}
def xcorr(a,b,L=6): return {str(k):float(a.corr(b.shift(k))) for k in range(-L,L+1)}   # k>0: b'nin k ay oncesi, yani b oncul
ASIAEXP=(np.log(Me[[i for i in Me.columns if P.continent.get(i)=='Asia & Pacific']].sum(axis=1))); ASIAEXP=(ASIAEXP-ASIAEXP.shift(12)).iloc[12:]
def impg(r): s=np.log(Mi[[i for i in Mi.columns if P.continent.get(i)==r]].sum(axis=1)); return (s-s.shift(12)).iloc[12:]
R['leadlag']={'eu_import_vs_asia_export':xcorr(impg('Europe'),ASIAEXP),'na_import_vs_asia_export':xcorr(impg('North America'),ASIAEXP),'eu_calls_vs_asia_calls':xcorr(RA['Europe'],RA['Asia & Pacific']),'na_calls_vs_asia_calls':xcorr(RA['North America'],RA['Asia & Pacific'])}
# 9. kuresel faktorun zaman serisi ozeti
F=pd.Series(Fc_,index=Gc.index.astype(str)); A=pd.Series(aggc.values,index=Gc.index.astype(str))
R['factor']=dict(min=[F.idxmin(),float(F.min())],max=[F.idxmax(),float(F.max())],agg_min=[A.idxmin(),float(A.min())],agg_max=[A.idxmax(),float(A.max())],ac1=float(F.autocorr(1)),ac12=float(F.autocorr(12)),agg_by_year={str(y):float(A[[k for k in A.index if k.startswith(str(y))]].mean()) for y in range(2020,2027)},F_by_year={str(y):float(F[[k for k in F.index if k.startswith(str(y))]].mean()) for y in range(2020,2027)})
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('w/out/results.json','w'),ensure_ascii=False,indent=1)
rd=lambda d_:{k:(round(v,3) if isinstance(v,float) else v) for k,v in d_.items()}
for k in('calls','thr','cont','bulk','tank'):
    o=R[k]; print(k,'N',o['N'],'T',o['T'],'ev',[round(x,3) for x in o['ev']],'mp',round(o['mp'],3),'null ev1',round(o['null_ev1'],3),round(o['null_ev1_p99'],3),'null glob',round(o['null_glob_med'],3)); print('  med',rd(o['med']),'w',rd(o['wmean']),'gt10',round(o['share_glob_gt10'],2),'gt25',round(o['share_glob_gt25'],2),'ncty',o['n_cty'],o['n_with_cty'],round(o['cty_share_among'],3))
    print('  agg',rd(o['agg'])); print('  post2021',rd(o['post2021']),'rbar',round(o['rbar_all'],4)); print('  cont',{a:rd(b) for a,b in o['by_cont'].items()}); print('  size',[rd(x) for x in o['by_size']]); print('  roll',[(r[0],round(r[1],3),round(r[2],3)) for r in o['roll'][::3]])
print('corrF',R['corr_F_calls_thr']); print('det',{k:(rd(v) if isinstance(v,dict) else v) for k,v in R['det'].items()}); print('regcorr',R['reg_corr']); print('leadlag',{k:{a:round(b,2) for a,b in v.items()} for k,v in R['leadlag'].items()}); print(R['factor'])
'''

# ==============================================================================
# w/w_f2.py
# ==============================================================================
MODULES['w/w_f2.py'] = r'''
# Konteyner ugrak etkeni ile genel liman etkeninin karsilastirilmasi (etkenler arasi korelasyon, tepe ve dip aylari, one gecme-gecikme)
import pandas as pd, numpy as np, json
c=pd.read_csv('w/out/factors_calls.csv'); k=pd.read_csv('w/out/factors_container.csv'); t=pd.read_csv('w/out/factors_thr.csv'); b=pd.read_csv('w/out/factors_bulk.csv'); ta=pd.read_csv('w/out/factors_tanker.csv')
R=json.load(open('w/out/results.json'))
o={}
for n,f in(('cont',k),('thr',t),('bulk',b),('tank',ta)): o['F_calls_'+n]=float(np.corrcoef(c.Fg,f.Fg)[0,1]); o['agg_calls_'+n]=float(np.corrcoef(c['agg'],f['agg'])[0,1])
o['cont_min']=[k.m[k.Fg.idxmin()],float(k.Fg.min())]; o['cont_max']=[k.m[k.Fg.idxmax()],float(k.Fg.max())]
o['cont_agg_min']=[k.m[k['agg'].idxmin()],float(k['agg'].min())]; o['cont_agg_max']=[k.m[k['agg'].idxmax()],float(k['agg'].max())]
# l > 0: genel etken konteyner etkeninin l ay onunde
o['ll_cont_calls']={str(l):float(np.corrcoef(c.Fg.values[max(0,-l):len(c)-max(0,l)],k.Fg.values[max(0,l):len(k)-max(0,-l)])[0,1]) for l in range(-12,13,3)}
k['y']=k.m.str[:4]; o['cont_agg_by_year']={y:float(v) for y,v in k.groupby('y')['agg'].mean().items()}; o['cont_F_by_year']={y:float(v) for y,v in k.groupby('y').Fg.mean().items()}
R['f2']=o; json.dump(R,open('w/out/results.json','w'),ensure_ascii=False,indent=1); print(o)
'''

# ==============================================================================
# w/w_extra.py
# ==============================================================================
MODULES['w/w_extra.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
import statsmodels.api as sm
d=pd.read_parquet('data/pwd/daily.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid')
d=d[d.date<='2026-10-31'].copy(); d['thr']=d['import']+d['export']; d['m']=d.date.dt.to_period('M'); R=json.load(open('w/out/results.json')); X={}
def build(col,minmean,end=None,start=None):
    M=d.groupby(['portid','m'])[col].sum().unstack('portid'); keep=(M.iloc[:12].sum()>0)&(M.mean()>=minmean)&((M==0).mean()<0.25); M=M.loc[:,keep]
    G=(2*(M-M.shift(12))/(M+M.shift(12))).iloc[12:].fillna(0)
    if end: G=G[G.index<=end]
    if start: G=G[G.index>=start]
    return M,G
def decomp(col,minmean,end=None,start=None,nnull=100):
    M,G=build(col,minmean,end,start); T,N=G.shape; ids=G.columns; w=(M.mean()/M.mean().sum()).values; Gv=G.values; Z=(Gv-Gv.mean(0))/np.where(Gv.std(0)>0,Gv.std(0),1)
    reg=P.continent.reindex(ids).values; iso=P.ISO3.reindex(ids).values
    U,S,Vt=np.linalg.svd(Z,full_matrices=False); F=U[:,0]*np.sqrt(T); agg=(Gv*w).sum(1)
    if np.corrcoef(F,agg)[0,1]<0: F=-F
    def proj(Y,f): Xc=np.column_stack([np.ones(len(f)),f]); B=np.linalg.lstsq(Xc,Y,rcond=None)[0]; fit=Xc@B; return fit-fit.mean(0),Y-fit
    Gc=Gv-Gv.mean(0); cg,E1=proj(Gc,F); cr=np.zeros_like(Gc); E2=E1.copy()
    for r in pd.unique(reg[pd.notna(reg)]):
        idx=np.where(reg==r)[0]
        if len(idx)<10: continue
        Zr=E1[:,idx]/np.where(E1[:,idx].std(0)>0,E1[:,idx].std(0),1); u=np.linalg.svd(Zr,full_matrices=False)[0][:,0]*np.sqrt(T); c,e=proj(E1[:,idx],u); cr[:,idx]=c; E2[:,idx]=e
    cc=np.zeros_like(Gc); E3=E2.copy()
    for c_ in pd.unique(iso[pd.notna(iso)]):
        idx=np.where(iso==c_)[0]
        if len(idx)<8: continue
        Zr=E2[:,idx]/np.where(E2[:,idx].std(0)>0,E2[:,idx].std(0),1); u=np.linalg.svd(Zr,full_matrices=False)[0][:,0]*np.sqrt(T); c,e=proj(E2[:,idx],u); cc[:,idx]=c; E3[:,idx]=e
    A={k:(v*w).sum(1) for k,v in(('glob',cg),('reg',cr),('cty',cc),('idio',E3))}; va=np.var(agg)
    out=dict(N=int(N),T=int(T),ev1=float(S[0]**2/(T*N)),sd_agg=float(agg.std()),**{'agg_'+k:float(np.var(v)/va) for k,v in A.items()},agg_cov=float(1-sum(np.var(v) for v in A.values())/va),
        port_glob_med=float(np.median(1-E1.var(0)/Gc.var(0))),port_idio_med=float(np.median(E3.var(0)/Gc.var(0))),indep_idio_sd=float(np.sqrt(((w*E3.std(0))**2).sum())),idio_sd=float(A['idio'].std()))
    top=np.argsort(-w)[:10]; out['agg_idio_top10']=float(np.var((E3[:,top]*w[top]).sum(1))/va); top=np.argsort(-w)[:100]; out['agg_idio_top100']=float(np.var((E3[:,top]*w[top]).sum(1))/va)
    rng=np.random.default_rng(5); ne=[]
    for k in range(nnull):
        Zs=np.column_stack([np.roll(Z[:,j],rng.integers(1,T)) for j in range(N)]); ne.append(np.linalg.svd(Zs,compute_uv=False)[0]**2/(T*N))
    out['null_ev1']=float(np.mean(ne)); out['null_ev1_p99']=float(np.percentile(ne,99)); out['p_ev1']=float((np.array(ne)>=out['ev1']).mean())
    return out,(G,F,E3,w,ids)
X['calls'],(G,F,E3,w,ids)=decomp('portcalls',10); X['thr'],_=decomp('thr',1000); X['container'],_=decomp('portcalls_container',5)
X['calls_to_2025m06'],_=decomp('portcalls',10,end='2025-06'); X['calls_from_2022'],_=decomp('portcalls',10,start='2022-01'); X['container_from_2022'],_=decomp('portcalls_container',5,start='2022-01'); X['calls_to_2025m12'],_=decomp('portcalls',10,end='2025-12')
# gostergeler: buyuk limanlar arasinda kuresel faktorle en yuksek korelasyon
S=pd.read_csv('w/out/shares_calls.csv',index_col=0); S['name']=P.fullname.reindex(S.index); big=S[S['mean']>=S['mean'].quantile(.8)]
X['bell']=[[r['name'],float(r['mean']),float(r.load),float(r.glob)] for _,r in big.sort_values('load',ascending=False).head(12).iterrows()]; X['anti']=[[r['name'],float(r['mean']),float(r.load)] for _,r in big.sort_values('load').head(6).iterrows()]
X['top20']=[[r['name'],float(r['mean']),float(r.glob),float(r.reg),float(r.cty),float(r.idio)] for _,r in S.sort_values('mean',ascending=False).head(20).iterrows()]
# kuresel faktorun dunya toplamini bir ay onceden ongorme gucu ve en buyuk N limanla faktorun yeniden kurulmasi
Z=(G.values-G.values.mean(0))/G.values.std(0); 
rec={}
for k in(10,25,50,100,250):
    top=np.argsort(-w)[:k]; f=Z[:,top].mean(1); rec[str(k)]=float(np.corrcoef(f,F)[0,1])
X['recon']=rec
rng=np.random.default_rng(1); X['recon_random']={str(k):float(np.mean([abs(np.corrcoef(Z[:,rng.choice(len(w),k,replace=False)].mean(1),F)[0,1]) for _ in range(200)])) for k in(10,25,50,100,250)}
R['x']=X; json.dump(R,open('w/out/results.json','w'),ensure_ascii=False,indent=1)
rd=lambda o:{k:(round(v,3) if isinstance(v,float) else v) for k,v in o.items()}
for k in('calls','thr','container','calls_to_2025m06','calls_to_2025m12','calls_from_2022','container_from_2022'): print(k,rd(X[k]))
print([[b[0][:28],round(b[1]),round(b[2],2),round(b[3],2)] for b in X['bell']]); print([[b[0][:28],round(b[2],2)] for b in X['anti']]); print(X['recon'],X['recon_random'])
print([[b[0][:22],round(b[1]),round(b[2],2),round(b[3],2),round(b[4],2),round(b[5],2)] for b in X['top20']])
'''

# ==============================================================================
# w/w_extra2.py
# ==============================================================================
MODULES['w/w_extra2.py'] = r'''
import pandas as pd, numpy as np, json
R=json.load(open('w/out/results.json')); S=pd.read_csv('w/out/shares_calls.csv',index_col=0); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid'); S['country']=P.country.reindex(S.index)
Y={}
g=S[S.has_cty].groupby('country'); c=pd.DataFrame({'n':g.size(),'cty':g.cty.median(),'glob':g.glob.median(),'idio':g.idio.median(),'w':g.w.sum()}).sort_values('cty',ascending=False)
Y['cty_top']=[[k,int(r.n),float(r.cty),float(r.glob),float(r.idio)] for k,r in c.head(8).iterrows()]; Y['cty_bottom']=[[k,int(r.n),float(r.cty),float(r.glob),float(r.idio)] for k,r in c.tail(5).iterrows()]
big=c.sort_values('w',ascending=False).head(10); Y['cty_big']=[[k,int(r.n),float(r.cty),float(r.glob),float(r.idio),float(r.w)] for k,r in big.iterrows()]
dom=S[['glob','reg','cty','idio']].idxmax(axis=1).value_counts(normalize=True); Y['dominant']={k:float(v) for k,v in dom.items()}
nonidio=S[['glob','reg','cty']]; d2=nonidio.idxmax(axis=1).value_counts(normalize=True); Y['dominant_common']={k:float(v) for k,v in d2.items()}
Y['common_gt50']=float(((1-S.idio)>0.5).mean()); Y['common_med']=float((1-S.idio).median()); Y['common_p90']=float((1-S.idio).quantile(.9)); Y['common_p10']=float((1-S.idio).quantile(.1))
f=pd.read_csv('w/out/factors_calls.csv'); rc=[c_ for c_ in f.columns if c_.startswith('reg_')]
Y['reg']={c_[4:]:dict(min=[f.m[f[c_].idxmin()],float(f[c_].min())],max=[f.m[f[c_].idxmax()],float(f[c_].max())],ac1=float(f[c_].autocorr())) for c_ in rc}
cm=f[rc].corr().values; Y['reg_factor_corr_absmax']=float(np.abs(cm[np.triu_indices(len(rc),1)]).max())
Y['wglob_cont']={k:v['wglob'] for k,v in R['calls']['by_cont'].items()}
R['y']=Y; json.dump(R,open('w/out/results.json','w'),ensure_ascii=False,indent=1)
for k,v in Y.items(): print(k,v)
'''

# ==============================================================================
# w/w_extra3.py
# ==============================================================================
MODULES['w/w_extra3.py'] = r'''
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
d=pd.read_parquet('data/pwd/daily.parquet'); d=d[d.date<='2026-10-31'].copy(); d['m']=d.date.dt.to_period('M'); R=json.load(open('w/out/results.json'))
def run(col,minmean,lab):
    M=d.groupby(['portid','m'])[col].sum().unstack('portid'); keep=(M.iloc[:12].sum()>0)&(M.mean()>=minmean)&((M==0).mean()<0.25); M=M.loc[:,keep]
    G=(2*(M-M.shift(12))/(M+M.shift(12))).iloc[12:].fillna(0); F=pd.read_csv(f'w/out/factors_{lab}.csv').Fg.values; Gv=G.values; T,N=Gv.shape; size=M.mean().values; out={}
    for h in(1,3,6,12):
        r0=[];r1=[]
        for j in range(N):
            y=Gv[h:,j]; x=Gv[:-h,j]; f=F[:-h]
            if y.std()==0: r0.append(np.nan); r1.append(np.nan); continue
            X0=np.column_stack([np.ones(len(y)),x]); e0=y-X0@np.linalg.lstsq(X0,y,rcond=None)[0]; X1=np.column_stack([X0,f]); e1=y-X1@np.linalg.lstsq(X1,y,rcond=None)[0]
            k0=1-e0.var()/y.var(); k1=1-e1.var()/y.var(); n=len(y); r0.append(1-(1-k0)*(n-1)/(n-2)); r1.append(1-(1-k1)*(n-1)/(n-3))
        r0=np.array(r0); r1=np.array(r1); dq=r1-r0; q=pd.qcut(size,5,labels=False)
        out[str(h)]=dict(ar=float(np.nanmedian(r0)),arf=float(np.nanmedian(r1)),gain_med=float(np.nanmedian(dq)),gain_mean=float(np.nanmean(dq)),share_pos=float(np.nanmean(dq>0)),share_gt5=float(np.nanmean(dq>0.05)),gain_top=float(np.nanmedian(dq[q==4])),gain_bottom=float(np.nanmedian(dq[q==0])))
    # toplam icin
    w=size/size.sum(); agg=(Gv*w).sum(1); o={}
    for h in(1,3,6):
        y=agg[h:]; X0=np.column_stack([np.ones(len(y)),agg[:-h]]); e0=y-X0@np.linalg.lstsq(X0,y,rcond=None)[0]; o[str(h)]=float(1-e0.var()/y.var())
    out['agg_ar']=o; return out
R['fc']={'calls':run('portcalls',10,'calls'),'container':run('portcalls_container',5,'container')}; json.dump(R,open('w/out/results.json','w'),ensure_ascii=False,indent=1)
for k,v in R['fc'].items():
    for h,o in v.items(): print(k,h,{a:round(b,3) for a,b in o.items()})
'''

# ==============================================================================
# w/w_cont.py
# ==============================================================================
MODULES['w/w_cont.py'] = r'''
# Konteyner ugrak dongusunun ayrismasi: tonaj, ugrak basina yuk ve tedarik zinciri baskisi (GSCPI)
import pandas as pd, numpy as np, json, warnings; warnings.filterwarnings('ignore')
d=pd.read_parquet('data/pwd/daily.parquet'); P=pd.read_csv('data/pw/ports_global.csv').set_index('portid'); d=d[d.date<='2026-10-31'].copy(); d['m']=d.date.dt.to_period('M')
d['tc']=d.import_container+d.export_container; d['cn']=d.portcalls-d.portcalls_container; d['tn']=d['import']+d['export']-d.tc
R=json.load(open('w/out/results.json')); Q={}
gs=pd.read_excel('data/gs/gscpi.xlsx',header=None,skiprows=5,names=['d','g']).dropna(); gs['m']=pd.to_datetime(gs.d,format='%d-%b-%Y').dt.to_period('M'); gs=gs.set_index('m').g.astype(float)
def dh(s): return (2*(s-s.shift(12))/(s+s.shift(12))).iloc[12:]
W=d.groupby('m')[['portcalls_container','tc','cn','tn','portcalls']].sum(); G=pd.DataFrame({'gC':dh(W.portcalls_container),'gT':dh(W.tc),'gCn':dh(W.cn),'gTn':dh(W.tn)}); G['gap']=G.gT-G.gC; G['gapn']=G.gTn-G.gCn
fc=pd.read_csv('w/out/factors_calls.csv'); fk=pd.read_csv('w/out/factors_container.csv'); idx=pd.PeriodIndex(fc.m,freq='M'); G=G.reindex(idx); G['Fc']=fc.Fg.values; G['Fk']=fk.Fg.values; G['gs']=gs.reindex(idx).values; G['aggk']=fk['agg'].values
# konteyner tonaji icin liman duzeyi kuresel etken
M=d.groupby(['portid','m']).tc.sum().unstack('portid'); keep=(M.iloc[:12].sum()>0)&(M.mean()>=1000)&((M==0).mean()<0.25); M=M.loc[:,keep]; Gp=dh(M).fillna(0); Z=((Gp-Gp.mean())/Gp.std(ddof=0)).fillna(0).values; T,N=Z.shape
U,S,Vt=np.linalg.svd(Z,full_matrices=False); Fkt=U[:,0]*np.sqrt(T); w=(M.mean()/M.mean().sum()).values; aggT=(Gp.values*w).sum(1)
if np.corrcoef(Fkt,aggT)[0,1]<0: Fkt=-Fkt
G['Fkt']=Fkt; rng=np.random.default_rng(4); ne=[np.linalg.svd(np.column_stack([np.roll(Z[:,j],rng.integers(1,T)) for j in range(N)]),compute_uv=False)[0]**2/(T*N) for _ in range(100)]
Q['tonnage_factor']=dict(N=int(N),ev1=float(S[0]**2/(T*N)),null=float(np.mean(ne)),null_p99=float(np.percentile(ne,99)))
cc=lambda a,b:float(np.corrcoef(G[a],G[b])[0,1])
Q['corr']={k:cc(*k.split('|')) for k in('Fk|Fkt','Fc|Fkt','Fc|Fk','gC|gT','gap|gs','Fk|gs','Fkt|gs','Fc|gs','gC|gs','gT|gs','gapn|gs','gC|Fc','gT|Fc','gap|Fc','gCn|gTn')}
def ll(a,b,L=range(-12,13,3)):  # b, a'nin l ay gerisinde
    o={}
    for l in L:
        x=G[a].values; y=G[b].values
        if l>=0: o[str(l)]=float(np.corrcoef(x[:len(x)-l] if l else x,y[l:])[0,1])
        else: o[str(l)]=float(np.corrcoef(x[-l:],y[:l])[0,1])
    return o
Q['ll_Fc_Fkt']=ll('Fc','Fkt'); Q['ll_Fc_Fk']=ll('Fc','Fk'); Q['ll_gs_Fk']=ll('gs','Fk',range(-6,13,3)); Q['ll_gs_gap']=ll('gs','gap',range(-6,13,3))
# baskidan arindirma: Fk ~ GSCPI (0 ve 6 ay gecikmeli); artigin genel etkenle iliskisi
import statsmodels.api as sm
X=pd.DataFrame({'gs':G.gs,'gs6':pd.Series(gs.reindex(idx-6).values,index=idx)}); m=sm.OLS(G.Fk,sm.add_constant(X)).fit(cov_type='HAC',cov_kwds={'maxlags':12}); Q['Fk_on_gs']=dict(b0=float(m.params['gs']),t0=float(m.tvalues['gs']),b6=float(m.params['gs6']),t6=float(m.tvalues['gs6']),r2=float(m.rsquared))
m1=sm.OLS(G.gC,sm.add_constant(G[['gT']])).fit(); m2=sm.OLS(G.gC,sm.add_constant(G[['gT','gs']])).fit(cov_type='HAC',cov_kwds={'maxlags':12}); Q['gC_on_gT']=dict(b=float(m1.params['gT']),r2=float(m1.rsquared)); Q['gC_on_gT_gs']=dict(bT=float(m2.params['gT']),bgs=float(m2.params['gs']),tgs=float(m2.tvalues['gs']),r2=float(m2.rsquared))
m3=sm.OLS(G.gap,sm.add_constant(G[['gs']])).fit(cov_type='HAC',cov_kwds={'maxlags':12}); Q['gap_on_gs']=dict(b=float(m3.params['gs']),t=float(m3.tvalues['gs']),r2=float(m3.rsquared))
m4=sm.OLS(G.gapn,sm.add_constant(G[['gs']])).fit(cov_type='HAC',cov_kwds={'maxlags':12}); Q['gapn_on_gs']=dict(b=float(m4.params['gs']),t=float(m4.tvalues['gs']),r2=float(m4.rsquared))
# yillik ortalamalar
G['y']=[p_.year for p_ in idx]; Q['by_year']={str(y):dict(gC=float(g.gC.mean()),gT=float(g.gT.mean()),gap=float(g.gap.mean()),gs=float(g.gs.mean()),gCn=float(g.gCn.mean()),gTn=float(g.gTn.mean()),Fk=float(g.Fk.mean()),Fkt=float(g.Fkt.mean()),Fc=float(g.Fc.mean())) for y,g in G.groupby('y')}
Q['gap_max']=[str(G.gap.idxmax()),float(G.gap.max())]; Q['gap_min']=[str(G.gap.idxmin()),float(G.gap.min())]; Q['gs_max']=[str(G.gs.idxmax()),float(G.gs.max())]; Q['Fkt_min']=[str(G.Fkt.idxmin()),float(G.Fkt.min())]; Q['Fkt_max']=[str(G.Fkt.idxmax()),float(G.Fkt.max())]
# bolgeler
d['cont']=P.continent.reindex(d.portid).values; Wr=d.groupby(['cont','m'])[['portcalls_container','tc']].sum(); reg={}
for c in Wr.index.levels[0]:
    a=dh(Wr.loc[c].portcalls_container).reindex(idx); b=dh(Wr.loc[c].tc).reindex(idx); gp=(b-a); reg[c]=dict(corr_gap_gs=float(np.corrcoef(gp,G.gs)[0,1]),corr_gC_gs=float(np.corrcoef(a,G.gs)[0,1]),corr_gC_gT=float(np.corrcoef(a,b)[0,1]),gap2021=float(gp[G.y==2021].mean()),gap2023=float(gp[G.y==2023].mean()),gC2021=float(a[G.y==2021].mean()),gC2023=float(a[G.y==2023].mean()))
Q['reg']=reg
# liman duzeyi: konteyner ugrak buyumesinin GSCPI ile korelasyonu
Mc=d.groupby(['portid','m']).portcalls_container.sum().unstack('portid'); keep=(Mc.iloc[:12].sum()>0)&(Mc.mean()>=5)&((Mc==0).mean()<0.25); Mc=Mc.loc[:,keep]; Gc=dh(Mc).fillna(0).reindex(idx)
r=np.array([np.corrcoef(Gc[c_].values,G.gs.values)[0,1] if Gc[c_].std()>0 else np.nan for c_ in Gc.columns]); size=Mc.mean().values; q=pd.qcut(size,5,labels=False)
Q['port_gs']=dict(N=int(len(r)),med=float(np.nanmedian(r)),share_neg=float(np.nanmean(r<0)),share_lt_m03=float(np.nanmean(r<-0.3)),by_size=[dict(q=int(k)+1,med=float(np.nanmedian(r[q==k])),neg=float(np.nanmean(r[q==k]<0))) for k in range(5)],wmean=float(np.nansum(r*size)/np.nansum(size*(~np.isnan(r)))))
rn=[]
for k in range(200):
    s_=np.roll(G.gs.values,rng.integers(6,T-6)); rn.append(np.nanmedian([np.corrcoef(Gc[c_].values,s_)[0,1] for c_ in Gc.columns[:300]]))
Q['port_gs']['null_med_sd']=float(np.std(rn)); Q['port_gs']['p']=float((np.array(rn)<=Q['port_gs']['med']).mean())
G.assign(m=idx.astype(str)).to_csv('w/out/container_divergence.csv',index=False)
R['q']=Q; json.dump(R,open('w/out/results.json','w'),ensure_ascii=False,indent=1,allow_nan=False)
rd=lambda o:{k:(round(v,3) if isinstance(v,float) else (rd(v) if isinstance(v,dict) else v)) for k,v in o.items()}
for k,v in Q.items(): print(k,rd(v) if isinstance(v,dict) else v)
'''

# ==============================================================================
# w/w_fig.py
# ==============================================================================
MODULES['w/w_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
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('w/out/results.json')); X=R['x']
B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'; Gn='#1e8449'; Gy='#909497'
L={'calls':'Tüm\nuğraklar','cont':'Konteyner','bulk':'Kuru\ndökme','tank':'Tanker','thr':'Yük\nhacmi'}
def ser(lab):
    f=pd.read_csv(f'w/out/factors_{lab}.csv'); f.index=pd.PeriodIndex(f.m,freq='M').to_timestamp(); return f
fc=ser('calls'); fk=ser('container'); ft=ser('thr')
# S1 kuresel faktor ve dunya toplam buyumesi
fig,ax=plt.subplots(2,1,figsize=(7.4,4.6),sharex=True)
a=ax[0]; a.plot(fc.index,fc.Fg,color=B,lw=1.4,label='Tüm uğraklar'); a.plot(fk.index,fk.Fg,color=O,lw=1.1,label='Konteyner uğrakları'); a.plot(ft.index,ft.Fg,color=Gn,lw=1.0,ls='--',label='Yük hacmi'); a.axhline(0,color='#999',lw=.6)
a.set_ylabel('Küresel faktör (standart birim)'); a.legend(frameon=False,fontsize=7.5,ncol=3,loc='lower right'); a.set_title('(a) Küresel liman faktörü',fontsize=9,loc='left')
a=ax[1]; a.plot(fc.index,100*fc['agg'],color=B,lw=1.4,label='Tüm uğraklar'); a.plot(fk.index,100*fk['agg'],color=O,lw=1.1,label='Konteyner uğrakları'); a.plot(ft.index,100*ft['agg'],color=Gn,lw=1.0,ls='--',label='Yük hacmi'); a.axhline(0,color='#999',lw=.6)
a.set_ylabel('Yıllık büyüme (%)'); a.set_title('(b) Örneklem toplamının ağırlıklı yıllık büyümesi',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('w/fig/Figure_1.png'); plt.close()
# S2 PC1 payi ve sifir dagilimi; liman duzeyi kuresel pay dagilimi
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; ks=['calls','cont','tank','bulk','thr']; x=np.arange(len(ks)); a.bar(x-.2,[100*R[k]['ev'][0] for k in ks],.4,color=B,label='Gözlenen'); a.bar(x+.2,[100*R[k]['null_ev1'] for k in ks],.4,color='#c8d0d6',label='Dairesel kaydırma (ortalama)')
a.scatter(x+.2,[100*R[k]['null_ev1_p99'] for k in ks],marker='_',color=Rd,s=120,zorder=3,label='Sıfır dağılımı, %99')
a.set_xticks(x); a.set_xticklabels([L[k] for k in ks],fontsize=7.5); a.set_ylabel('Birinci bileşenin varyans payı (%)'); a.legend(frameon=False,fontsize=7); a.set_title('(a) Ortak bileşenin gücü',fontsize=9,loc='left')
a=ax[1]; S=pd.read_csv('w/out/shares_calls.csv'); Sk=pd.read_csv('w/out/shares_container.csv'); bins=np.arange(0,0.71,0.025)
a.hist(100*S.glob,bins=100*bins,color=B,alpha=.75,label='Tüm uğraklar',density=True); a.hist(100*Sk.glob,bins=100*bins,histtype='step',color=O,lw=1.3,label='Konteyner',density=True)
a.axvline(100*R['calls']['null_glob_med'],color=Rd,lw=.9,ls=':'); a.text(100*R['calls']['null_glob_med']+1,0.085,'sıfır dağılımı\nmedyanı',fontsize=7,color=Rd)
a.set_xlabel('Küresel faktörün açıkladığı varyans (%)'); a.set_ylabel('Yoğunluk'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Limanlar arasında dağılım',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('w/fig/Figure_2.png'); plt.close()
# S3 varyans paylari: gemi tipi ve buyukluk
fig,ax=plt.subplots(1,2,figsize=(7.4,3.2)); cols=[B,O,Gn,'#d5d8dc']; nm=['Küresel','Bölgesel','Ülkesel','Limana özgü']
a=ax[0]; ks=['calls','cont','tank','bulk','thr']; bot=np.zeros(len(ks))
for j,k in enumerate(('glob','reg','cty','idio')):
    v=np.array([100*R[q]['mean'][k] for q in ks]); a.bar(range(len(ks)),v,bottom=bot,color=cols[j],label=nm[j],width=.65); bot+=v
a.set_xticks(range(len(ks))); a.set_xticklabels([L[k] for k in ks],fontsize=7.5); a.set_ylabel('Ortalama varyans payı (%)'); a.set_ylim(0,118); a.legend(frameon=False,fontsize=7,ncol=4,loc='upper center',columnspacing=.8,handlelength=1); a.set_title('(a) Seriye göre',fontsize=9,loc='left')
a=ax[1]; q=[1,2,3,4,5]; a.plot(q,[100*b['glob'] for b in R['calls']['by_size']],'o-',color=B,label='Tüm uğraklar'); a.plot(q,[100*b['glob'] for b in R['cont']['by_size']],'s-',color=O,label='Konteyner'); a.plot(q,[100*b['glob'] for b in R['thr']['by_size']],'^--',color=Gn,label='Yük hacmi')
a.set_xticks(q); a.set_xlabel('Liman büyüklüğü beşte birlik dilimi (1 = en küçük)'); a.set_ylabel('Küresel payın medyanı (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Büyüklüğe göre küresel pay',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('w/fig/Figure_3.png'); plt.close()
# S4 toplamin ayristirilmasi ve kayan pencere
fig,ax=plt.subplots(1,2,figsize=(7.4,3.2))
a=ax[0]; ks=['calls','container','thr']; lab=['Tüm uğraklar','Konteyner','Yük hacmi']; bot=np.zeros(3); c5=[B,O,Gn,'#d5d8dc','#7d6608']; n5=['Küresel','Bölgesel','Ülkesel','Limana özgü','Kovaryans']
for j,k in enumerate(('agg_glob','agg_reg','agg_cty','agg_idio','agg_cov')):
    v=np.array([100*max(X[q][k],0) for q in ks]); a.bar(range(3),v,bottom=bot,color=c5[j],label=n5[j],width=.6); bot+=v
a.set_xticks(range(3)); a.set_xticklabels(lab); a.set_ylabel('Toplam büyüme varyansındaki pay (%)'); a.legend(frameon=False,fontsize=7,loc='upper center',ncol=3,bbox_to_anchor=(.5,1.02)); a.set_ylim(0,125); a.set_title('(a)',fontsize=9,loc='left')
a=ax[1]
for k,c,l in(('calls',B,'Tüm uğraklar'),('cont',O,'Konteyner'),('tank',Rd,'Tanker'),('bulk',Gy,'Kuru dökme')):
    r=R[k]['roll']; a.plot(pd.PeriodIndex([z[0] for z in r],freq='M').to_timestamp(),[100*z[1] for z in r],color=c,lw=1.2,label=l)
a.set_ylabel('Birinci bileşenin payı (%), 36 aylık pencere'); a.set_xlabel('Pencerenin son ayı'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Zaman içinde eş zamanlılık',fontsize=9,loc='left'); import matplotlib.dates as md; a.xaxis.set_major_locator(md.YearLocator()); a.xaxis.set_major_formatter(md.DateFormatter('%Y'))
fig.tight_layout(); fig.savefig('w/fig/Figure_4.png'); plt.close()
# S5 one gecme-gecikme
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0),sharey=True); ll=R['leadlag']; k=list(range(-6,7))
a=ax[0]; a.plot(k,[ll['eu_import_vs_asia_export'][str(i)] for i in k],'o-',color=B,label='Avrupa'); a.plot(k,[ll['na_import_vs_asia_export'][str(i)] for i in k],'s-',color=O,label='Kuzey Amerika'); a.axvline(0,color='#999',lw=.6); a.axhline(0,color='#999',lw=.6)
a.set_xlabel('Asya ihracatının öne geçtiği ay sayısı'); a.set_ylabel('Korelasyon'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Konteyner ithalatı – Asya konteyner ihracatı',fontsize=9,loc='left')
a=ax[1]; a.plot(k,[ll['eu_calls_vs_asia_calls'][str(i)] for i in k],'o-',color=B,label='Avrupa'); a.plot(k,[ll['na_calls_vs_asia_calls'][str(i)] for i in k],'s-',color=O,label='Kuzey Amerika'); a.axvline(0,color='#999',lw=.6); a.axhline(0,color='#999',lw=.6)
a.set_xlabel('Asya uğraklarının öne geçtiği ay sayısı'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Toplam uğraklar',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('w/fig/Figure_5.png'); plt.close()
# S6-S7 konteyner ayrismasi
D=pd.read_csv('w/out/container_divergence.csv'); D.index=pd.PeriodIndex(D.m,freq='M').to_timestamp(); Q=R['q'] if 'q' in R else json.load(open('w/out/results.json'))['q']
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; a.plot(D.index,D.Fk,color=O,lw=1.5,label='Konteyner uğrakları'); a.plot(D.index,D.Fkt,color=Rd,lw=1.5,label='Konteyner tonajı'); a.plot(D.index,D.Fc,color=B,lw=1.0,ls='--',label='Tüm uğraklar'); a.axhline(0,color='#999',lw=.6)
a.set_ylabel('Küresel etken (standart birim)'); a.legend(frameon=False,fontsize=7,loc='lower right'); a.set_title('(a) Üç küresel etken',fontsize=9,loc='left')
a=ax[1]; a.plot(D.index,D.gs,color='#555',lw=1.4,label='Tedarik zinciri baskı endeksi'); a.plot(D.index,-D.Fk,color=O,lw=1.3,label='Konteyner uğrak etkeni (ters çevrilmiş)'); a.axhline(0,color='#999',lw=.6); a.set_ylabel('Standart birim'); a.legend(frameon=False,fontsize=7,loc='upper right'); a.set_ylim(-3,6.2); a.set_title('(b) Baskı endeksi ve uğrak etkeni',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('w/fig/Figure_6.png'); plt.close()
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; a.scatter(100*D.gCn,100*D.gTn,s=12,color=B,alpha=.6,lw=0,label='Konteyner dışı (r = %s)'%f"{Q['corr']['gCn|gTn']:.2f}".replace('.',',')); a.scatter(100*D.gC,100*D.gT,s=12,color=O,alpha=.75,lw=0,label='Konteyner (r = %s)'%f"{Q['corr']['gC|gT']:.2f}".replace('.',',')); a.axhline(0,color='#999',lw=.6); a.axvline(0,color='#999',lw=.6)
a.set_xlabel('Uğrak sayısı, yıllık büyüme (%)'); a.set_ylabel('Tonaj, yıllık büyüme (%)'); a.legend(frameon=False,fontsize=7.5,loc='upper left'); a.set_title('(a) Aylık dünya toplamları',fontsize=9,loc='left')
a=ax[1]; by=Q['by_year']; ys=sorted(by); x=np.arange(len(ys)); a.bar(x-.2,[100*by[y]['gC'] for y in ys],.4,color=O,label='Konteyner uğrakları'); a.bar(x+.2,[100*by[y]['gT'] for y in ys],.4,color=Rd,label='Konteyner tonajı'); a.axhline(0,color='#999',lw=.6)
a.set_xticks(x); a.set_xticklabels(ys); a.set_ylabel('Yıllık ortalama büyüme (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(b) Yıllara göre',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('w/fig/Figure_7.png'); plt.close()
'''

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