# -*- coding: utf-8 -*-

# ------------------------------------------------------------------------------
# # Limanlar için yük faktörü ve eşzamanlılık faktörü: altyapı verisi olmadan kapasite kullanımının ölçülmesi, rastgele varış sınaması ve Türkiye liman sisteminde havuzlama potansiyeli, 2019–2025
#
# *Load factor and coincidence factor for ports: measuring capacity utilisation without infrastructure data, a random-arrival test and the pooling potential of Türkiye's port system, 2019–2025*
#
# **GOMD Yayınlar · Deniz ekonomisi 04** — https://gomd.pirireis.edu.tr/publications.php?id=me04
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür. Veri, replikasyon paketinde (11 Kasım 2026 indirmesi) varsa oradan; yoksa özgün sağlayıcıdan (UNCTADstat, IMF PortWatch) indirilir. Bütün sayılar `out/results.json` dosyasına yazılır.
#
# This notebook reproduces **every number, table and figure** of the paper. Data are taken from the replication package (retrieved 11 November 2026) if present; otherwise they are downloaded from the original provider. All numbers are written to `out/results.json`.
#
# *Makale hazırlık aşamasındadır; dergi adı ve tam metin kabulden sonra eklenecektir. / The paper is in preparation; journal name and full text will be added after acceptance.*
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# ## 0. Kurulum / Setup
# ------------------------------------------------------------------------------
# [Colab Magic] !pip install -q pandas numpy statsmodels scipy matplotlib adjustText py7zr

# ------------------------------------------------------------------------------
# ## 1. Veri / Data
# Exact reproduction in Colab: upload the replication package (.zip) from the GOMD Publications page, unzip it, then run the cells.
# ------------------------------------------------------------------------------
VERI_KAYNAGI = 'kaynak'   # Colab'da 'zip' secenegi icin not defterine bakin / see the notebook for the 'zip' option in Colab
# Veri / Data. Paketteki veri varsa onu kullanır (11 Kasım 2026 indirmesi, birebir yeniden üretim);
# yoksa özgün kaynaktan indirir (güncel sürüm; sağlayıcı revizyonları nedeniyle sonuçlar küçük farklar gösterebilir).
# Uses the data shipped in the replication package if present (retrieved 11 Nov 2026, exact reproduction);
# otherwise downloads from the original provider (current release; results may differ slightly after provider revisions).
import os, glob, json, urllib.request
import pandas as pd
os.makedirs('data', exist_ok=True); os.makedirs('out', exist_ok=True); os.makedirs('fig', exist_ok=True)
def _unctad(code):
    d = 'data/US_' + code
    if glob.glob(d + '/*.csv'): return
    arc = d + '.7z'
    if not os.path.exists(arc):
        print('UNCTADstat indiriliyor / downloading:', code)
        req = urllib.request.Request('https://unctadstat-api.unctad.org/bulkdownload/US.%s/US_%s' % (code, code), headers={'User-Agent': 'gomd-replication/1.0'})
        with urllib.request.urlopen(req, timeout=300) as r, open(arc, 'wb') as f: f.write(r.read())
    import py7zr
    with py7zr.SevenZipFile(arc) as z: z.extractall(d)
for _c in ['PortCalls']: _unctad(_c)
def _pw_get(params, layer='Daily_Ports_Data'):
    import urllib.request, urllib.parse, time
    url = 'https://services9.arcgis.com/weJ1QsnbMYJlCHdG/arcgis/rest/services/' + layer + '/FeatureServer/0/query?' + urllib.parse.urlencode(params)
    for k in range(6):
        try:
            with urllib.request.urlopen(urllib.request.Request(url, headers={'User-Agent': 'gomd-replication/1.0'}), timeout=120) as r:
                d = json.load(r)
            if 'error' in d: raise RuntimeError(str(d['error'])[:200])
            return d
        except Exception as e:
            err = e; time.sleep(3 * (k + 1))
    raise err

def _pw_layer(layer, where):
    rows, o = [], 0
    while True:
        f = _pw_get(dict(where=where, outFields='*', orderByFields='ObjectId', resultOffset=o, resultRecordCount=1000, f='json'), layer)['features']
        rows += [x['attributes'] for x in f]; o += 1000
        if len(f) < 1000: return pd.DataFrame(rows)
_PW = 'data/portwatch/'; os.makedirs(_PW, exist_ok=True)
if not os.path.exists(_PW + 'ports.csv'):
    print('IMF PortWatch liman veritabani indiriliyor / downloading ports database')
    _pw_layer('PortWatch_ports_database', '1=1').to_csv(_PW + 'ports.csv', index=False)
for _iso in ['TUR', 'GRC', 'EGY', 'ISR', 'ROU', 'BGR', 'GEO', 'MLT', 'CYP', 'LBN', 'UKR', 'MAR', 'ITA', 'ESP', 'DZA', 'TUN', 'HRV', 'SVN']:
    if os.path.exists(_PW + 'daily_%s.csv' % _iso): continue
    print('IMF PortWatch gunluk liman verisi indiriliyor / downloading daily port data:', _iso)
    _pts = pd.read_csv(_PW + 'ports.csv'); _pts = _pts[_pts.ISO3 == _iso].portid.tolist()
    pd.concat([_pw_layer('Daily_Ports_Data', "portid='%s'" % _x) for _x in _pts]).to_csv(_PW + 'daily_%s.csv' % _iso, index=False)

# ------------------------------------------------------------------------------
# ## 2. Haftalık seri, yük faktörü ve rastgele varış modeli / Weekly series, load factor and random-arrival model
# ------------------------------------------------------------------------------
# Makale 4: ortak işlevler. PortWatch günlük liman verisinden haftalık seri, gözlenen kapasite, kullanım, düzenlilik
import pandas as pd, numpy as np, glob
PW='data/portwatch/'
def load(iso):
    d=pd.read_csv(PW+f'daily_{iso}.csv',parse_dates=['date']); d['q']=d['import']+d['export']; d['qc']=d.import_container+d.export_container; d['qt']=d.import_tanker+d.export_tanker
    d['w']=d.date.dt.to_period('W-SUN').dt.start_time; return d
def weekly(d):
    w=d.groupby(['portid','w']).agg(q=('q','sum'),calls=('portcalls','sum'),qc=('qc','sum'),qt=('qt','sum'),cc=('portcalls_container','sum')).reset_index()
    w['y']=w.w.dt.year; return w[(w.w>='2018-12-31')&(w.w<='2025-12-28')]          # 2019-2025 tam haftalar (ilk hafta 31 Aralık 2018)
def port_year(w, top=0.10):
    w=w.copy(); w.loc[w.w=='2018-12-31','y']=2019; R=[]
    for (p,y),g in w.groupby(['portid','y']):
        q=g.q.values; n=len(q); k=max(int(round(n*top)),1); pk=np.sort(q)[-k:].mean()
        R.append(dict(portid=p,y=y,n=n,qbar=q.mean(),peak=pk,calls=g.calls.mean(),disp=(g.calls.var(ddof=1)/g.calls.mean()) if g.calls.mean()>0 else np.nan,tpc=q.sum()/g.calls.sum() if g.calls.sum()>0 else np.nan,csh=g.qc.sum()/q.sum() if q.sum()>0 else np.nan,tsh=g.qt.sum()/q.sum() if q.sum()>0 else np.nan,cv=q.std(ddof=1)/q.mean() if q.mean()>0 else np.nan))
    A=pd.DataFrame(R).sort_values(['portid','y']); A['K']=A.groupby('portid').peak.cummax(); A['U']=A.qbar/A.K; A['Uin']=A.qbar/A.peak; return A
def null_U(d, rng, sims=400, top=0.10):
    # Kısıtsız rastgele varış (Poisson) altında beklenen yıl içi tepe oranı: liman-yıl bazında
    d=d[(d.date>='2018-12-31')&(d.date<='2025-12-28')].copy(); d['y']=d.w.dt.year; d.loc[d.w=='2018-12-31','y']=2019; out=[]
    for (p,y),g in d.groupby(['portid','y']):
        nw=g.w.nunique(); lam=g.portcalls.sum()/nw; gg=g[g.portcalls>0]
        if lam<=0 or len(gg)<5: out.append(dict(portid=p,y=y,U0=np.nan)); continue
        per=(gg.q/gg.portcalls).values; wt=gg.portcalls.values/gg.portcalls.sum(); k=max(int(round(nw*top)),1)
        N=rng.poisson(lam,size=(sims,nw)); tot=N.sum(); draws=rng.choice(per,size=tot,p=wt); idx=np.repeat(np.arange(sims*nw),N.ravel()); Q=np.bincount(idx,weights=draws,minlength=sims*nw).reshape(sims,nw)
        m=Q.mean(1); pk=np.sort(Q,axis=1)[:,-k:].mean(1); out.append(dict(portid=p,y=y,U0=float(np.mean(m/pk))))
    return pd.DataFrame(out)

# ------------------------------------------------------------------------------
# ## 3. Analizler ve şekiller / Analyses and figures
# ------------------------------------------------------------------------------
# Makale 4: limanlar için yük faktörü, rastgele varış (Poisson) karşılaştırması, eşzamanlılık faktörü ve havuzlama; 2023 depremi; Akdeniz-Karadeniz karşılaştırması
import json, os, warnings, statsmodels.api as sm, matplotlib
warnings.filterwarnings('ignore')
from scipy import stats
matplotlib.use('Agg'); import matplotlib.pyplot as plt
O_='out/'; F_='fig/'; R={}; rng=np.random.default_rng(20261111)
TRN={'Istanbul':'İstanbul (Ambarlı batısı)','Nemrut Bay':'Nemrut Körfezi','Yakacik':'Yakacık','Ambarli':'Ambarlı','Botas':'BOTAŞ Ceyhan','Mersin':'Mersin','Tekirdag':'Tekirdağ','Aliaga':'Aliağa','Yarimca':'Yarımca','Diliskelesi':'Diliskelesi','Botas Natural Gas Terminal':'BOTAŞ Marmara Ereğlisi','Gemlik':'Gemlik','Zonguldak':'Zonguldak','Samsun':'Samsun','Eregli':'Karadeniz Ereğli','Delta Terminal':'Delta Terminal','Toros Gubre':'Toros Gübre','Derince (Izmit-Kocaeli)':'Derince','Karabiga':'Karabiga','Gebze':'Gebze','Port Akdeniz':'Antalya (Port Akdeniz)','Borusan':'Borusan (Gemlik)','Izmir':'İzmir','Iskenderun':'İskenderun','Bandirma':'Bandırma','Golcuk':'Gölcük','Güllük':'Güllük','Trabzon':'Trabzon','Tuzla (Pendik)':'Tuzla (Pendik)','Yalova':'Yalova','Sakarya (Karasu)':'Karasu','Haydarpasa':'Haydarpaşa','Saraylar':'Saraylar','Canakkale':'Çanakkale','Gulf of Saros FSRU':'Saros FSRU','Hopa':'Hopa','Tasucu':'Taşucu','Cesme':'Çeşme'}
REG={'Ege':['Cesme','Nemrut Bay','Aliaga','Izmir','Güllük'],'İstanbul ve Batı Marmara':['Tekirdag','Botas Natural Gas Terminal','Ambarli','Istanbul','Haydarpasa','Tuzla (Pendik)','Gulf of Saros FSRU','Canakkale'],'Güney Marmara':['Karabiga','Saraylar','Bandirma','Borusan','Gemlik','Yalova'],'İzmit Körfezi':['Diliskelesi','Gebze','Yarimca','Derince (Izmit-Kocaeli)','Golcuk'],'Akdeniz':['Botas','Toros Gubre','Delta Terminal','Iskenderun','Yakacik','Mersin','Tasucu','Port Akdeniz'],'Karadeniz':['Sakarya (Karasu)','Eregli','Zonguldak','Samsun','Trabzon','Hopa']}
d=load('TUR'); nm=d.drop_duplicates('portid').set_index('portid').portname; reg={p:r for r,ps in REG.items() for p in ps}; assert set(nm.values)==set(reg), set(nm.values)^set(reg)
w=weekly(d); w['y2']=w.y.where(w.w!='2018-12-31',2019); A=port_year(w).merge(null_U(d,rng),on=['portid','y']); A['rho']=A.Uin/A.U0; A['en']=A.portid.map(nm); A['ad']=A.en.map(TRN); A['bolge']=A.en.map(reg)
A.round(4).to_csv(O_+'liman_yil.csv',index=False); R['n_port']=int(nm.size); R['date_max']=str(d.date.max().date())
YR=list(range(2019,2026)); d['y']=d.date.dt.year; dd=d[d.y.isin(YR)]
# A) Sistem: tonaj, çağrı, çağrı başına yük; segmentler
seg=['container','dry_bulk','general_cargo','roro','tanker']
S=dd.groupby('y').agg(q=('q','sum'),calls=('portcalls','sum'),imp=('import','sum'),exp=('export','sum'),**{f'q_{s}':(f'import_{s}','sum') for s in seg},**{f'x_{s}':(f'export_{s}','sum') for s in seg},**{f'c_{s}':(f'portcalls_{s}','sum') for s in seg})
for s in seg: S[f't_{s}']=S[f'q_{s}']+S[f'x_{s}']; S[f'tpc_{s}']=S[f't_{s}']/S[f'c_{s}']
S['tpc']=S.q/S.calls; R['sys']={int(y):{k:float(v) for k,v in r.items()} for y,r in S.iterrows()}
R['sys_growth']={k:round(float(np.log(S.loc[2025,k]/S.loc[2019,k])*100),1) for k in ['q','calls','tpc']+[f't_{s}' for s in seg]+[f'c_{s}' for s in seg]+[f'tpc_{s}' for s in seg]}
R['seg_share_2025']={s:round(float(S.loc[2025,f't_{s}']/S.loc[2025,'q']*100),1) for s in seg}
# B) Liman tablosu 2025 ve 2019
def prow(r): return dict(en=r.en,ad=r.ad,bolge=r.bolge,mt=round(r.qbar*r.n/1e6,2),calls=round(r.calls,1),tpc=round(r.tpc,0),csh=round(r.csh*100,1),tsh=round(r.tsh*100,1),U=round(r.Uin,3),U0=round(r.U0,3) if pd.notna(r.U0) else None,rho=round(r.rho,3) if pd.notna(r.rho) else None,disp=round(r.disp,2) if pd.notna(r.disp) else None,cv=round(r.cv,3) if pd.notna(r.cv) else None)
a25=A[A.y==2025].sort_values('qbar',ascending=False); a19=A[A.y==2019].set_index('en')
R['ports25']=[prow(r) for _,r in a25.iterrows()]; R['ports19']={r.en:prow(r) for _,r in A[A.y==2019].iterrows()}
big=A[(A.qbar*52>=1e6)&A.U0.notna()]; R['n_py_big']=int(len(big)); R['n_port_big25']=int(((a25.qbar*52)>=1e6).sum())
R['hhi']={int(y):round(float(((g.qbar/g.qbar.sum())**2).sum()*1e4),0) for y,g in A.groupby('y')}; R['top5']={int(y):round(float(g.qbar.nlargest(5).sum()/g.qbar.sum()*100),1) for y,g in A.groupby('y')}
R['reg_share']={int(y):{r:round(float(g[g.bolge==r].qbar.sum()/g.qbar.sum()*100),1) for r in REG} for y,g in A[A.y.isin([2019,2025])].groupby('y')}
# C) Ölçek yasası: yük faktörü ~ ln(haftalık çağrı); sıfır modelin açıklama gücü
m1=sm.OLS(big.Uin,sm.add_constant(np.log(big.calls))).fit(cov_type='cluster',cov_kwds={'groups':big.portid}); m2=sm.OLS(big.Uin,sm.add_constant(big.U0)).fit(cov_type='cluster',cov_kwds={'groups':big.portid})
m3=sm.OLS(big.Uin,sm.add_constant(pd.DataFrame({'U0':big.U0,'csh':big.csh,'tsh':big.tsh}))).fit(cov_type='cluster',cov_kwds={'groups':big.portid})
R['scale']=dict(b_lncalls=round(float(m1.params.iloc[1]),4),se=round(float(m1.bse.iloc[1]),4),r2=round(float(m1.rsquared),3),b_U0=round(float(m2.params.iloc[1]),3),se_U0=round(float(m2.bse.iloc[1]),3),a_U0=round(float(m2.params.iloc[0]),3),r2_U0=round(float(m2.rsquared),3),n=int(len(big)),nport=int(big.portid.nunique()),
  m3={k:[round(float(m3.params[k]),3),round(float(m3.bse[k]),3),round(float(m3.pvalues[k]),4)] for k in ['U0','csh','tsh']},r2_m3=round(float(m3.rsquared),3))
R['rho']=dict(med=round(float(big.rho.median()),3),q25=round(float(big.rho.quantile(.25)),3),q75=round(float(big.rho.quantile(.75)),3),share_lt1=round(float((big.rho<1).mean()*100),1),share_gt105=round(float((big.rho>1.05).mean()*100),1),wil_p=float(stats.wilcoxon(big.groupby('portid').rho.mean()-1).pvalue),n_ports=int(big.portid.nunique()),ports_gt1=int((big.groupby('portid').rho.mean()>1).sum()))
pm=big.groupby('en').agg(rho=('rho','mean'),U=('Uin','mean'),U0=('U0','mean'),calls=('calls','mean'),disp=('disp','mean'),n=('y','count')).sort_values('rho',ascending=False); pm['ad']=pm.index.map(TRN); R['rho_ports']=[dict(ad=r.ad,rho=round(r.rho,3),U=round(r.U,3),U0=round(r.U0,3),disp=round(r.disp,2),n=int(r.n)) for _,r in pm.iterrows()]
R['disp']=dict(med=round(float(big.disp.median()),2),share_lt1=round(float((big.disp<1).mean()*100),1))
R['U_by_year']={int(y):dict(med=round(float(g.Uin.median()),3),rho=round(float(g.rho.median()),3)) for y,g in big.groupby('y')}
# tavan sınaması: haftalık elleçlemenin çarpıklığı (kapasite tavanı varsa dağılım sola çarpık olur)
sk=w.merge(big[['portid','y']].rename(columns={'y':'y2'}),on=['portid','y2']).groupby(['portid','y2']).q.apply(lambda x: float(stats.skew(x))); R['skew']=dict(med=round(float(sk.median()),3),share_pos=round(float((sk>0).mean()*100),1),n=int(len(sk)))
# D) Eşzamanlılık: sistem yük faktörü, eşzamanlılık faktörü (eşzamanlı tepe / liman tepeleri toplamı), havuzlama kazancı
def pk(x,top=.10): k=max(int(round(len(x)*top)),1); return np.sort(np.asarray(x))[-k:].mean()
def coin(wk,col='q'):
    out={}
    for y,g in wk.groupby('y2'):
        P=g.pivot_table(index='w',columns='portid',values=col,aggfunc='sum').fillna(0); tot=P.sum(axis=1); sp=sum(pk(P[c]) for c in P.columns); cp=pk(tot)
        out[int(y)]=dict(mean=float(tot.mean()),cpeak=float(cp),speak=float(sp),LF=round(float(tot.mean()/cp),3),CF=round(float(cp/sp),3),save=float(sp-cp),LFavg=round(float(tot.mean()/sp),3),n=int((P.sum()>0).sum()))
    return out
R['coin_sys']=coin(w); w['en']=w.portid.map(nm); w['bolge']=w.en.map(reg)
R['coin_reg']={r:coin(w[w.bolge==r]) for r in REG}
cont_ports=a25[(a25.csh>=.2)].portid.tolist(); R['cont_ports']=[TRN[nm[p]] for p in cont_ports]; R['coin_cont']=coin(w[w.portid.isin(cont_ports)],'qc')
R['coin_cont_marmara']=coin(w[w.portid.isin(cont_ports)&w.bolge.isin(['İstanbul ve Batı Marmara','Güney Marmara','İzmit Körfezi'])],'qc')
R['coin_avg']={k:dict(LF=round(float(np.mean([v['LF'] for v in c.values()])),3),CF=round(float(np.mean([v['CF'] for v in c.values()])),3)) for k,c in [('sys',R['coin_sys'])]+[(r,R['coin_reg'][r]) for r in REG]+[('cont',R['coin_cont']),('cont_marmara',R['coin_cont_marmara'])]}
# bağımsızlık altında beklenen eşzamanlılık: liman serilerini yıl içinde rastgele kaydırarak (dairesel) 500 kez
def coin_null(wk,y,B=500):
    P=wk[wk.y2==y].pivot_table(index='w',columns='portid',values='q',aggfunc='sum').fillna(0).values; sp=sum(pk(P[:,j]) for j in range(P.shape[1])); n=P.shape[0]; cf=[]
    for b in range(B):
        Q=np.column_stack([np.roll(P[:,j],rng.integers(n)) for j in range(P.shape[1])]); cf.append(pk(Q.sum(1))/sp)
    return dict(mean=round(float(np.mean(cf)),3),p05=round(float(np.percentile(cf,5)),3),p95=round(float(np.percentile(cf,95)),3))
R['coin_null']={y:coin_null(w,y) for y in (2019,2022,2025)}
# haftalık liman serileri arası ortalama korelasyon (büyük limanlar)
bp=a25[(a25.qbar*52)>=2e6].portid; Pw=w[w.portid.isin(bp)].pivot_table(index='w',columns='portid',values='q').fillna(0); Cm=Pw.corr().values; R['corr_ports']=dict(mean=round(float(Cm[np.triu_indices_from(Cm,1)].mean()),3),n=int(len(bp)),share_neg=round(float((Cm[np.triu_indices_from(Cm,1)]<0).mean()*100),1))
# E) 2023 depremi (6 Şubat 2023): İskenderun, Yakacık; Akdeniz limanları ve ulusal toplam
q0=pd.Timestamp('2023-02-06'); W=w.pivot_table(index='w',columns='en',values='q').fillna(0); Wc=w.pivot_table(index='w',columns='en',values='calls').fillna(0)
base=W[(W.index<q0)&(W.index>=q0-pd.Timedelta(weeks=52))]; post=W[W.index>=q0]
def quake(p):
    b=base[p].mean(); x=post[p]; r=dict(base=float(b),base_calls=float(Wc[p][(Wc.index<q0)&(Wc.index>=q0-pd.Timedelta(weeks=52))].mean()))
    for k in (1,4,12,26,52): r[f'ratio_{k}']=round(float(x.iloc[:k].mean()/b),3); r[f'loss_{k}']=float((b-x.iloc[:k]).sum())
    rm=x.rolling(8).mean()/b; rec=rm[rm>=.9]; r['rec_weeks']=int((rec.index[0]-q0).days/7)+1 if len(rec) else None
    r['ratio_2024']=round(float(W[p][W.index.year==2024].mean()/b),3); r['ratio_2025']=round(float(W[p][W.index.year==2025].mean()/b),3); r['calls_4']=round(float(Wc[p][Wc.index>=q0].iloc[:4].mean()),1)
    # plasebo: aynı ölçünün (12 haftalık ortalama / önceki 52 hafta) diğer tüm tarihlerdeki dağılımı
    s=W[p]; pl=[(s.iloc[i:i+12].mean()/s.iloc[i-52:i].mean()) for i in range(52,len(s)-12) if not (q0-pd.Timedelta(weeks=12)<=s.index[i]<=q0+pd.Timedelta(weeks=60))]
    r['placebo_p']=round(float((np.array(pl)<=r['ratio_12']).mean()),4); r['placebo_p_hi']=round(float((np.array(pl)>=r['ratio_12']).mean()),4); r['peak_pre']=float(pk(base[p])); r['max_post12']=float(x.iloc[:12].max()); return r
R['quake']={TRN[p]:quake(p) for p in ['Iskenderun','Yakacik','Mersin','Botas','Toros Gubre','Delta Terminal','Port Akdeniz','Tasucu']}
oth=[c for c in W.columns if c not in ('Iskenderun','Yakacik')]; tot=W[oth].sum(axis=1); bt=tot[(tot.index<q0)&(tot.index>=q0-pd.Timedelta(weeks=52))].mean(); R['quake_others']={f'ratio_{k}':round(float(tot[tot.index>=q0].iloc[:k].mean()/bt),3) for k in (4,12,52)}
gulf=['Iskenderun','Yakacik']; g=W[gulf].sum(axis=1); bg=g[(g.index<q0)&(g.index>=q0-pd.Timedelta(weeks=52))].mean(); R['quake_gulf']=dict(base=float(bg),**{f'loss_{k}':float((bg-g[g.index>=q0].iloc[:k]).sum()) for k in (4,12,52)})
med=['Mersin','Port Akdeniz','Tasucu','Botas','Toros Gubre','Delta Terminal']; mm=W[med].sum(axis=1); bm=mm[(mm.index<q0)&(mm.index>=q0-pd.Timedelta(weeks=52))].mean(); R['quake_med']=dict(base=float(bm),**{f'excess_{k}':float((mm[mm.index>=q0].iloc[:k]-bm).sum()) for k in (4,12,52)})
# konteyner: İskenderun -> Mersin
Wq=w.pivot_table(index='w',columns='en',values='qc').fillna(0); bq=Wq[(Wq.index<q0)&(Wq.index>=q0-pd.Timedelta(weeks=52))]
R['quake_cont']={TRN[p]:dict(base=float(bq[p].mean()),**{f'ratio_{k}':round(float(Wq[p][Wq.index>=q0].iloc[:k].mean()/bq[p].mean()),3) for k in (4,12,52)},**{f'diff_{k}':float((Wq[p][Wq.index>=q0].iloc[:k]-bq[p].mean()).sum()) for k in (4,12,52)}) for p in ['Iskenderun','Mersin']}
# F) Karşılaştırma ülkeleri
isos=sorted(os.path.basename(f)[6:9] for f in glob.glob(PW+'daily_*.csv') if 'TUR' not in f); R['peer_isos']=isos; PR=[]; CC={}
ports=pd.read_csv(PW+'ports.csv').set_index('portid')
for iso in isos:
    dp=load(iso); wp=weekly(dp); wp['y2']=wp.y.where(wp.w!='2018-12-31',2019); Ap=port_year(wp).merge(null_U(dp,rng,sims=200),on=['portid','y']); Ap['iso']=iso; PR.append(Ap); CC[iso]=coin(wp)
if PR:
    PA=pd.concat(PR+[A.assign(iso='TUR')[PR[0].columns]]); PA['rho']=PA.Uin/PA.U0; PA['tr']=(PA.iso=='TUR').astype(float); PA['name']=PA.portid.map(ports.portname); PA.round(4).to_csv(O_+'karsilastirma_liman_yil.csv',index=False)
    pb=PA[(PA.qbar*52>=1e6)&PA.U0.notna()&(PA.calls>=2)]; R['peer_n']=dict(py=int(len(pb)),ports=int(pb.portid.nunique()),tr_ports=int(pb[pb.tr==1].portid.nunique()),countries=int(pb.iso.nunique()))
    X=pd.DataFrame({'U0':pb.U0,'csh':pb.csh,'tsh':pb.tsh,'tr':pb.tr}); X=pd.concat([X,pd.get_dummies(pb.y,prefix='y',drop_first=True).astype(float)],axis=1); mp=sm.OLS(pb.Uin,sm.add_constant(X)).fit(cov_type='cluster',cov_kwds={'groups':pb.portid})
    R['peer_reg']={k:[round(float(mp.params[k]),4),round(float(mp.bse[k]),4),round(float(mp.pvalues[k]),4)] for k in ['U0','csh','tsh','tr']}; R['peer_reg']['r2']=round(float(mp.rsquared),3); R['peer_reg']['n']=int(mp.nobs)
    mr=sm.OLS(pb.rho,sm.add_constant(pd.concat([pb[['csh','tsh','tr']],pd.get_dummies(pb.y,prefix='y',drop_first=True).astype(float)],axis=1))).fit(cov_type='cluster',cov_kwds={'groups':pb.portid}); R['peer_rho_reg']={k:[round(float(mr.params[k]),4),round(float(mr.bse[k]),4),round(float(mr.pvalues[k]),4)] for k in ['csh','tsh','tr']}
    R['peer_rho']={('TUR' if t else 'diger'):dict(med=round(float(g.rho.median()),3),U=round(float(g.Uin.median()),3),U0=round(float(g.U0.median()),3),n=int(len(g))) for t,g in pb.groupby('tr')}
    CC['TUR']=R['coin_sys']; R['peer_coin']={iso:dict(LF=round(float(np.mean([v['LF'] for v in c.values()])),3),CF=round(float(np.mean([v['CF'] for v in c.values()])),3),LFavg=round(float(np.mean([v['LFavg'] for v in c.values()])),3),n=int(c[2025]['n']),mt=round(float(c[2025]['mean']*52/1e6),1)) for iso,c in CC.items()}
    R['peer_scale']=dict(r2_U0=round(float(sm.OLS(pb.Uin,sm.add_constant(pb.U0)).fit().rsquared),3))
# G) UNCTAD: limanda geçen ortanca süre
u=pd.read_csv('data/US_PortCalls/US_PortCalls.csv',dtype={'Economy':str,'CommercialMarket':str}); u=u[u.Economy.str.len()==3]; tc='Median time in port (days)'
def tip(mk,y): x=u[(u.Year==y)&(u.CommercialMarket==mk)].set_index('Economy')[tc].dropna(); return dict(tr=round(float(x['792']),3),med=round(float(x.median()),3),rank=int(x.rank()['792']),n=int(len(x)))
R['tip']={mk:{y:tip(mk,y) for y in (2018,2023)} for mk in ['00','03','05','02','04']}
uw=pd.read_csv('data/US_PortCalls/US_PortCalls.csv',dtype={'Economy':str,'CommercialMarket':str}); R['tip_world']={mk:{int(y):float(v) for y,v in uw[(uw.Economy=='0000')&(uw.CommercialMarket==mk)].set_index('Year')[tc].items() if y in (2018,2023)} for mk in ['00','03']}
# H) Sağlamlık: tepe tanımı (en yüksek %5, %10, %20 hafta) ve toplama penceresi (1, 2, 4 hafta)
def robust(top,agg):
    dd2=d[(d.date>='2018-12-31')&(d.date<='2025-12-28')].copy(); k=((dd2.date-pd.Timestamp('2018-12-31')).dt.days//(7*agg)); dd2['w']=pd.Timestamp('2018-12-31')+pd.to_timedelta(k*7*agg,unit='D')
    ww=dd2.groupby(['portid','w']).agg(q=('q','sum'),calls=('portcalls','sum'),qc=('qc','sum'),qt=('qt','sum'),cc=('portcalls_container','sum')).reset_index(); ww['y']=ww.w.dt.year.where(ww.w!='2018-12-31',2019); ww.loc[ww.y==2018,'y']=2019
    Ax=port_year(ww,top).merge(null_U(dd2,rng,sims=200,top=top),on=['portid','y']); Ax['rho']=Ax.Uin/Ax.U0; b=Ax[(Ax.qbar*52/agg>=1e6)&Ax.U0.notna()]
    ww['y2']=ww.y; P=ww[ww.y2==2025].pivot_table(index='w',columns='portid',values='q',aggfunc='sum').fillna(0); cf=pk(P.sum(axis=1),top)/sum(pk(P[c],top) for c in P.columns)
    return dict(U=round(float(b.Uin.median()),3),rho=round(float(b.rho.median()),3),lt1=round(float((b.rho<1).mean()*100),1),r2=round(float(sm.OLS(b.Uin,sm.add_constant(b.U0)).fit().rsquared),3),CF=round(float(cf),3),n=int(len(b)))
R['robust']={f'{int(t*100)}_{g}':robust(t,g) for t,g in [(.05,1),(.10,1),(.20,1),(.10,2),(.10,4)]}
# konteyner limanlarında konteyner yükü için yük faktörü
wc=w[w.portid.isin(cont_ports)].copy(); R['cont_LF']={TRN[nm[p]]:{int(y):round(float(g.qc.mean()/pk(g.qc)),3) for y,g in wc[wc.portid==p].groupby('y2')} for p in cont_ports}
R['cont_mt25']={TRN[nm[p]]:round(float(wc[(wc.portid==p)&(wc.y2==2025)].qc.sum()/1e6),2) for p in cont_ports}
# liman yük faktöründe yıllar arası kalıcılık ve değişim
pv=big.pivot_table(index='portid',columns='y',values='Uin'); R['persist']=dict(corr_19_25=round(float(pv[[2019,2025]].dropna().corr().iloc[0,1]),3),n=int(len(pv[[2019,2025]].dropna())),icc=round(float(pv.dropna().mean(axis=1).var()/(pv.dropna().values.var())),3))
ch=(pv[2025]-pv[2019]).dropna().sort_values(); R['U_change']=dict(up=[dict(ad=TRN[nm[p]],d=round(float(v),3)) for p,v in ch.tail(4).items()][::-1],down=[dict(ad=TRN[nm[p]],d=round(float(v),3)) for p,v in ch.head(4).items()],med=round(float(ch.median()),3))
R['isk_annual']={int(y):round(float(W['Iskenderun'][W.index.year==y].sum()/1e6),2) for y in range(2019,2026)}
# ---------------- Şekiller
plt.rcParams.update({'font.size':9,'font.family':'DejaVu Sans','axes.spines.top':False,'axes.spines.right':False}); K1,K2,K3,K4,KG='#1f4e79','#c0392b','#2e8b57','#b8860b','#999999'
from adjustText import adjust_text
pm2=big.groupby('en').agg(U=('Uin','mean'),U0=('U0','mean'),calls=('calls','mean'),q=('qbar','mean')).reset_index(); pm2['ad']=pm2.en.map(TRN)
fig,ax=plt.subplots(figsize=(7.2,4.8)); ax.scatter(pm2.calls,pm2.U0,s=22,facecolors='none',edgecolors=KG,label='Rastgele varış modeli (beklenen)'); ax.scatter(pm2.calls,pm2.U,s=24,c=K1,label='Gözlenen yük faktörü',zorder=3)
for _,r in pm2.iterrows(): ax.plot([r.calls,r.calls],[r.U,r.U0],c=KG,lw=.5,zorder=1)
tx=[ax.text(r.calls,r.U,r.ad,fontsize=6.5) for _,r in pm2.iterrows() if r.q*52>=6e6 or r.calls<8]; ax.set_xscale('log'); adjust_text(tx,ax=ax,arrowprops=dict(arrowstyle='-',color='#bbbbbb',lw=.4))
ax.set_xlabel('Haftalık ortalama gemi çağrısı (logaritmik)'); ax.set_ylabel('Yük faktörü (ortalama / tepe)'); ax.legend(frameon=False,loc='lower right'); fig.tight_layout(); fig.savefig(F_+'Figure_1.png',dpi=300); plt.close()
if PR:
    g=pb.groupby(['portid','tr']).agg(U=('Uin','mean'),U0=('U0','mean')).reset_index(); fig,ax=plt.subplots(figsize=(5.6,5.2)); o=g[g.tr==0]; t=g[g.tr==1]
    ax.scatter(o.U0,o.U,s=14,c=KG,alpha=.7,label='Akdeniz–Karadeniz limanları (%d)'%len(o)); ax.scatter(t.U0,t.U,s=26,c=K2,label='Türkiye limanları (%d)'%len(t),zorder=3); ax.plot([.25,.95],[.25,.95],'--',c='k',lw=.7)
    ax.set_xlabel('Rastgele varış modelinin öngördüğü yük faktörü'); ax.set_ylabel('Gözlenen yük faktörü'); ax.set_xlim(.25,.95); ax.set_ylim(.25,.95); ax.legend(frameon=False,loc='upper left'); fig.tight_layout(); fig.savefig(F_+'Figure_2.png',dpi=300); plt.close()
rr=list(REG)+['Türkiye']; cr=[R['coin_reg'][r][2025] for r in REG]+[R['coin_sys'][2025]]; fig,axs=plt.subplots(1,2,figsize=(7.4,3.6),gridspec_kw={'width_ratios':[3,1.1]})
for ax,idx in ((axs[0],range(6)),(axs[1],[6])):
    x=np.arange(len(idx)); ax.bar(x-.2,[cr[i]['speak']/1e6 for i in idx],.4,color=KG,label='Liman tepelerinin toplamı'); ax.bar(x+.2,[cr[i]['cpeak']/1e6 for i in idx],.4,color=K1,label='Eşzamanlı tepe'); ax.scatter(x+.2,[cr[i]['mean']/1e6 for i in idx],marker='_',c='w',s=260,zorder=3,lw=1.6)
    for j,i in enumerate(idx): ax.text(j,max(cr[i]['speak'],cr[i]['cpeak'])/1e6*1.02,('EF = %.2f'%cr[i]['CF']).replace('.',','),ha='center',fontsize=7)
    ax.set_xticks(x); ax.set_xticklabels([rr[i].replace(' ve ','\nve ').replace('Güney ','Güney\n').replace('İzmit ','İzmit\n') for i in idx],fontsize=7)
axs[0].set_ylabel('Milyon ton / hafta'); axs[0].set_ylim(0,4.9); axs[0].legend(frameon=False,fontsize=7,loc='upper center',ncol=2); fig.tight_layout(); fig.savefig(F_+'Figure_3.png',dpi=300); plt.close()
fig,axs=plt.subplots(1,2,figsize=(7.4,3.4)); a=axs[0]; win=(W.index>=q0-pd.Timedelta(weeks=20))&(W.index<=q0+pd.Timedelta(weeks=60))
for p,c in (('Iskenderun',K2),('Yakacik',K4),('Mersin',K1)): a.plot((W.index[win]-q0).days/7,(W[p][win]/base[p].mean()*100).rolling(4,min_periods=1).mean(),c=c,label=TRN[p])
a.axvline(0,c='k',lw=.7,ls='--'); a.axhline(100,c=KG,lw=.6); a.set_xlabel('Depremden itibaren hafta'); a.set_ylabel('Deprem öncesi 52 hafta = 100'); a.set_title('(a) Haftalık elleçleme, 4 haftalık ortalama',fontsize=9); a.legend(frameon=False,fontsize=7)
a=axs[1]; yy=[2019,2020,2021,2022,2023,2024,2025]; isk=[W['Iskenderun'][W.index.year==y].sum()/1e6 for y in yy]; a.bar(yy,isk,color=[K1]*4+[K2]*3,width=.65); a.set_ylabel('Milyon ton'); a.set_title('(b) İskenderun, yıllık',fontsize=9); fig.tight_layout(); fig.savefig(F_+'Figure_4.png',dpi=300); plt.close()
fig,axs=plt.subplots(1,2,figsize=(7.4,3.3)); a=axs[0]
for k,c,l in (('q',K1,'Tonaj'),('calls',K2,'Gemi çağrısı'),('tpc',K3,'Çağrı başına yük')): a.plot(S.index,S[k]/S.loc[2019,k]*100,marker='o',ms=3,c=c,label=l)
a.axhline(100,c=KG,lw=.6); a.set_ylabel('2019 = 100'); a.legend(frameon=False,fontsize=7); a.set_title('(a) Türkiye liman sistemi',fontsize=9)
a=axs[1]; cs=R['coin_sys']; a.plot(YR,[cs[y]['LF'] for y in YR],marker='o',ms=3,c=K1,label='Sistem yük faktörü'); a.plot(YR,[cs[y]['CF'] for y in YR],marker='o',ms=3,c=K2,label='Eşzamanlılık faktörü'); a.plot(YR,[R['U_by_year'][y]['med'] for y in YR],marker='o',ms=3,c=K3,label='Liman yük faktörü (ortanca)'); a.set_ylim(.4,1); a.legend(frameon=False,fontsize=7); a.set_title('(b) Yük ve eşzamanlılık faktörleri',fontsize=9); fig.tight_layout(); fig.savefig(F_+'Figure_5.png',dpi=300); plt.close()
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,(np.integer,)): return int(o)
    if isinstance(o,(float,np.floating)): return None if (np.isnan(o) or np.isinf(o)) else float(o)
    return o
R=clean(R); json.dump(R,open(O_+'results.json','w'),ensure_ascii=False,indent=1)
s=open(O_+'results.json').read(); assert 'NaN' not in s
for k in ['sys_growth','seg_share_2025','hhi','top5','reg_share','scale','rho','disp','U_by_year','coin_avg','coin_null','corr_ports','quake_others','quake_gulf','quake_med','quake_cont','peer_n','peer_reg','peer_rho_reg','peer_rho','peer_coin','peer_scale','tip','tip_world']: print(k, json.dumps(R.get(k),ensure_ascii=False))
print(pd.DataFrame(R['ports25']).head(26).to_string()); print(pd.DataFrame(R['rho_ports']).to_string()); print(pd.DataFrame(R['quake']).T.to_string()); print(pd.DataFrame(R['coin_sys']).T.to_string()); print({r:R['coin_reg'][r]['2025'] for r in REG}); print(R['coin_cont']['2025'],R['coin_cont_marmara']['2025'],R['cont_ports'])

# ------------------------------------------------------------------------------
# ## Şekiller / Figures
# ------------------------------------------------------------------------------
print('Sekiller / figures:', sorted(glob.glob('fig/*.png')))