# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Dünya limanlarının takvimi: 2.065 limanda gemi uğraklarının haftalık ritmi, hafta sonu etkisi ve tatillerin küresel izi, 2019–2026
# The calendar of world ports: the weekly rhythm of ship calls, the weekend effect and the global footprint of holidays at 2,065 ports, 2019–2026
#
# GOMD Yayınlar · Makale 19 — https://gomd.pirireis.edu.tr/publications.php?id=rp19
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri `h/` klasörüne yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# veriyi `data/` ve `h/out/` klasörlerinden okur, sonuçları `h/out/` ve şekilleri `h/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to `h/` and runs them in order, each in its own process;
# it reads data from `data/` and `h/out/`, and writes results to `h/out/` and figures to `h/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib scipy statsmodels pyarrow
# Çalıştırma / Run:  python rp19_calendar_of_world_ports_weekend_holidays.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 = ['h/h_analysis.py', 'h/h_holiday.py', 'h/h_extra.py', 'h/h_extra2.py', 'h/h_robust.py', 'h/h_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

# ==============================================================================
# h/h_analysis.py
# ==============================================================================
MODULES['h/h_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
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['dow']=d.date.dt.dayofweek; d['iso']=d.portid.map(P.ISO3); d['cont']=d.portid.map(P.continent); d['yr']=d.date.dt.year
TY={'portcalls':'Tümü','portcalls_container':'Konteyner','portcalls_dry_bulk':'Kuru dökme','portcalls_tanker':'Tanker','portcalls_general_cargo':'Genel kargo','portcalls_roro':'Ro-Ro'}
DN=['Pzt','Sal','Çar','Per','Cum','Cmt','Paz']
R={'n_ports':int(d.portid.nunique()),'n_days':int(d.date.nunique()),'n_obs':int(len(d)),'first':str(d.date.min().date()),'last':str(d.date.max().date()),'calls_total':float(d.portcalls.sum())}
# ---- A. Kuresel haftalik profil
G=d.groupby('date')[list(TY)+['thr','import','export']].sum(); G['dow']=G.index.dayofweek; G['wk']=((G.index-G.index[0]).days//7)
def profile(s,dow,wk,B=1000,seed=1):
    # hafta icindeki paya gore: her hafta icin gun/hafta ortalamasi -> ortalama
    df=pd.DataFrame({'y':s.values,'dow':dow.values,'wk':wk.values}); w=df.pivot(index='wk',columns='dow',values='y').dropna(); r=w.div(w.mean(axis=1),axis=0)
    m=r.mean(); rng=np.random.default_rng(seed); bs=np.array([r.iloc[rng.integers(0,len(r),len(r))].mean().values for _ in range(B)])
    we=(r[[5,6]].mean(axis=1)/r[[0,1,2,3,4]].mean(axis=1)); web=np.array([we.iloc[rng.integers(0,len(we),len(we))].mean() for _ in range(B)])
    return dict(idx=[float(x) for x in m.values],lo=[float(x) for x in np.percentile(bs,2.5,axis=0)],hi=[float(x) for x in np.percentile(bs,97.5,axis=0)],wr=float(we.mean()),wr_lo=float(np.percentile(web,2.5)),wr_hi=float(np.percentile(web,97.5)),nweeks=int(len(r)))
R['global']={TY.get(k,k):profile(G[k],G.dow,G.wk) for k in list(TY)+['thr']}
# yillara gore kuresel hafta sonu orani
R['global_year']={int(y):profile(g.portcalls,g.dow,g.wk,B=300)['wr'] for y,g in G.groupby(G.index.year)}
R['global_year_thr']={int(y):profile(g.thr,g.dow,g.wk,B=300)['wr'] for y,g in G.groupby(G.index.year)}
# bolgelere gore
R['region']={}
for c,g in d.groupby('cont'):
    s=g.groupby('date')[['portcalls','portcalls_container','thr']].sum(); s['dow']=s.index.dayofweek; s['wk']=((s.index-s.index[0]).days//7)
    R['region'][c]={k:{a:b for a,b in profile(s[k],s.dow,s.wk,B=300).items() if a in('idx','wr','wr_lo','wr_hi')} for k in('portcalls','portcalls_container','thr')}; R['region'][c]['n']=int(g.portid.nunique()); R['region'][c]['share']=float(g.portcalls.sum()/d.portcalls.sum())
# ---- B. Liman duzeyi
pd_=d.groupby(['portid','dow']).portcalls.sum().unstack(); nd=d.groupby('dow').date.nunique(); rate=pd_/nd   # gunluk ortalama
tot=pd_.sum(axis=1); ok=tot>=1000
Q=pd.DataFrame({'tot':tot,'wr':rate[[5,6]].mean(axis=1)/rate[[0,1,2,3,4]].mean(axis=1)})
# ki-kare tekduzelik testi ve haftalik donemsellik gucu (Cramer V benzeri): chi2/(N)
exp=np.outer(tot,nd/nd.sum()); chi=((pd_.values-exp)**2/exp).sum(axis=1); Q['chi2']=chi; Q['p']=stats.chi2.sf(chi,6); Q['phi']=np.sqrt(chi/tot)
Q['peak']=rate.idxmax(axis=1); Q['trough']=rate.idxmin(axis=1); Q['amp']=(rate.max(axis=1)-rate.min(axis=1))/rate.mean(axis=1)
# hafta sonu etkisi icin Poisson z: hafta sonu payi beklenen 2/7 (gun sayisi duzeltmeli)
pwe=nd[[5,6]].sum()/nd.sum(); we=pd_[[5,6]].sum(axis=1); Q['z_we']=(we-tot*pwe)/np.sqrt(tot*pwe*(1-pwe))
for k in('ISO3','continent','country','portname','lat'): Q[k]=P[k]
ty=d.groupby('portid')[['portcalls_container','portcalls_dry_bulk','portcalls_tanker','portcalls_general_cargo','portcalls_roro']].sum(); Q['type']=ty.idxmax(axis=1).map(TY); sh=ty.div(ty.sum(axis=1),axis=0); Q['sh_cont']=sh.portcalls_container; Q['sh_roro']=sh.portcalls_roro; Q['sh_tank']=sh.portcalls_tanker; Q['sh_bulk']=sh.portcalls_dry_bulk
Q.to_csv('h/out/port_weekly.csv'); S=Q[ok].copy(); R['port']=dict(n=int(len(S)),share_calls=float(S.tot.sum()/tot.sum()),wr_med=float(S.wr.median()),wr_p10=float(S.wr.quantile(.1)),wr_p25=float(S.wr.quantile(.25)),wr_p75=float(S.wr.quantile(.75)),wr_p90=float(S.wr.quantile(.9)),
   sig_low=float(((S.z_we<-1.96)).mean()),sig_high=float(((S.z_we>1.96)).mean()),lt90=float((S.wr<0.9).mean()),lt80=float((S.wr<0.8).mean()),gt110=float((S.wr>1.1).mean()),
   nonuniform=float((S.p<0.01).mean()),amp_med=float(S.amp.median()),phi_med=float(S.phi.median()),
   peak_dist={DN[int(k)]:float(v) for k,v in S.peak.value_counts(normalize=True).sort_index().items()},trough_dist={DN[int(k)]:float(v) for k,v in S.trough.value_counts(normalize=True).sort_index().items()})
R['port_type']={k:dict(n=int(len(g)),wr_med=float(g.wr.median()),wr_w=float(np.average(g.wr,weights=g.tot)),lt90=float((g.wr<0.9).mean()),amp_med=float(g.amp.median()),phi_med=float(g.phi.median()),nonuniform=float((g.p<0.01).mean())) for k,g in S.groupby('type')}
S['sizeq']=pd.qcut(S.tot,5,labels=False); R['port_size']=[dict(q=int(k)+1,calls_med=float(g.tot.median()/R['n_days']),wr_med=float(g.wr.median()),amp_med=float(g.amp.median()),phi_med=float(g.phi.median())) for k,g in S.groupby('sizeq')]
# gelir grubu
try:
    wb=pd.read_csv('data/wb_country.csv'); col=[c for c in wb.columns if 'income' in c.lower()][0]; idc=[c for c in wb.columns if c.lower() in('iso3','id','iso3c','code')][0]; S['inc']=S.ISO3.map(wb.set_index(idc)[col])
    R['port_income']={str(k):dict(n=int(len(g)),wr_med=float(g.wr.median()),wr_w=float(np.average(g.wr,weights=g.tot))) for k,g in S.groupby('inc')}
except Exception as e: print('wb',e)
# regresyon: log wr ~ log size + paylar + bolge
S['lwr']=np.log(S.wr); S['ls']=np.log(S.tot/R['n_days']); S['abslat']=S.lat.abs()
SR=S.replace([np.inf,-np.inf],np.nan).dropna(subset=['lwr','ls','sh_cont','sh_roro','sh_tank','sh_bulk','continent','inc','ISO3']).copy(); SR['g']=pd.factorize(SR.ISO3)[0]
m=smf.ols('lwr ~ ls + sh_cont + sh_roro + sh_tank + sh_bulk + C(continent) + C(inc)',SR).fit(cov_type='cluster',cov_kwds={'groups':SR.g})
R['port_reg']={k:dict(b=float(m.params[k]),se=float(m.bse[k]),p=float(m.pvalues[k])) for k in m.params.index}; R['port_reg']['_r2']=float(m.rsquared); R['port_reg']['_n']=int(m.nobs)
# ---- C. Ulke duzeyi
cd=d.groupby(['iso','dow']).portcalls.sum().unstack(); ct=cd.sum(axis=1); cr=cd/nd; C=pd.DataFrame({'tot':ct,'wr':cr[[5,6]].mean(axis=1)/cr[[0,1,2,3,4]].mean(axis=1),'fri':cr[4]/cr.mean(axis=1),'sat':cr[5]/cr.mean(axis=1),'sun':cr[6]/cr.mean(axis=1),'mon':cr[0]/cr.mean(axis=1)}); C['name']=P.groupby('ISO3').country.first(); C=C[C.tot>=20000]
C.to_csv('h/out/country_weekly.csv'); R['country']=dict(n=int(len(C)),low=[[r['name'],float(r.wr),float(r.tot)] for _,r in C.sort_values('wr').head(12).iterrows()],high=[[r['name'],float(r.wr),float(r.tot)] for _,r in C.sort_values('wr').tail(8).iterrows()],
   big={r['name']:float(r.wr) for _,r in C.sort_values('tot',ascending=False).head(15).iterrows()},lt90=float((C.wr<0.9).mean()),med=float(C.wr.median()))
json.dump(R,open('h/out/results.json','w'),ensure_ascii=False,indent=1)
for k in('global_year','region','port','port_type','port_size','port_income','port_reg','country'): print(k,json.dumps(R.get(k),ensure_ascii=False)[:1400]); print()
print({k:(v['wr'],v['wr_lo'],v['wr_hi'],[round(x,3) for x in v['idx']]) for k,v in R['global'].items()})
'''

# ==============================================================================
# h/h_holiday.py
# ==============================================================================
MODULES['h/h_holiday.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['iso']=d.portid.map(P.ISO3); d['cont']=d.portid.map(P.continent)
R=json.load(open('h/out/results.json'))
H={'lny':['2019-02-05','2021-02-12','2022-02-01','2023-01-22','2024-02-10','2025-01-29','2026-02-17'],   # 2020 haric (salgin kapanmasi)
   'lny_all':['2019-02-05','2020-01-25','2021-02-12','2022-02-01','2023-01-22','2024-02-10','2025-01-29','2026-02-17'],
   'fitr':['2019-06-04','2020-05-24','2021-05-13','2022-05-02','2023-04-21','2024-04-10','2025-03-30','2026-03-20'],
   'adha':['2019-08-11','2020-07-31','2021-07-20','2022-07-09','2023-06-28','2024-06-16','2025-06-06','2026-05-27'],
   'xmas':['2019-12-25','2020-12-25','2021-12-25','2022-12-25','2023-12-25','2024-12-25','2025-12-25'],
   'golden':['2019-10-01','2020-10-01','2021-10-01','2022-10-01','2023-10-01','2024-10-01','2025-10-01','2026-10-01']}
H={k:pd.to_datetime(v) for k,v in H.items()}
GR={'Çin':['CHN'],'Doğu ve Güneydoğu Asya (Çin dışı, Ay Yeni Yılı kutlayan)':['VNM','KOR','SGP','MYS','TWN','HKG'],
    'Körfez ülkeleri':['SAU','ARE','QAT','KWT','BHR','OMN'],'Diğer Müslüman çoğunluklu ülkeler':['EGY','TUR','IDN','PAK','BGD','IRN','IRQ','JOR','DZA','MAR','TUN','LBY','MYS'],
    'Avrupa':None,'Kuzey Amerika':None,'Latin Amerika':None,'Japonya':['JPN'],'Dünya':'ALL'}
def series(g,cols=('portcalls','portcalls_container','thr','import_container','export_container','import','export')):
    if g=='Dünya': x=d
    elif g=='Avrupa': x=d[d.cont=='Europe']
    elif g=='Kuzey Amerika': x=d[d.iso.isin(['USA','CAN'])]
    elif g=='Latin Amerika': x=d[d.cont.isin(['South America'])|d.iso.isin(['MEX','PAN','COL','GTM','CRI','DOM','HND','SLV','NIC','CUB','JAM'])]
    else: x=d[d.iso.isin(GR[g])]
    return x.groupby('date')[list(cols)].sum()
def event(y,dates,lo=-21,hi=35,bins=None,excl=None):
    # Poisson GLM: gun-hafta + yil-ay sabit etkileri + goreli gun kuklalari (bins: [(ad,a,b)]) ; HAC
    df=pd.DataFrame({'y':y}); df['dow']=df.index.dayofweek; df['ym']=df.index.strftime('%Y')  # yil
    df['doy']=df.index.dayofyear
    rel=pd.Series(np.nan,index=df.index)
    for h in dates:
        k=(df.index-h).days; m=(k>=lo)&(k<=hi); rel[m]=k[m]
    X=pd.get_dummies(df.dow,prefix='d',drop_first=True).astype(float)
    X=X.join(pd.get_dummies(df.index.strftime('%Y'),prefix='y',drop_first=True).astype(float).set_index(df.index))
    # mevsimsellik: harmonikler (3) + yil icinde dogrusal egilim
    for j in(1,2,3):
        X['s%d'%j]=np.sin(2*np.pi*j*df.doy/365.25); X['c%d'%j]=np.cos(2*np.pi*j*df.doy/365.25)
    t=(df.index-df.index[0]).days/365.25; X['t']=t
    names=[]
    if bins is None: bins=[(str(k),k,k) for k in range(lo,hi+1)]
    for nm,a,b in bins:
        X['r_'+nm]=((rel>=a)&(rel<=b)).astype(float).values; names.append('r_'+nm)
    X=sm.add_constant(X); ok=np.ones(len(df),bool)
    if excl is not None: ok&=~((df.index>=excl[0])&(df.index<=excl[1]))
    m=sm.GLM(df.y[ok],X[ok],family=sm.families.Poisson()).fit(cov_type='HAC',cov_kwds={'maxlags':14})
    return {n[2:]:dict(b=float(m.params[n]),se=float(m.bse[n]),eff=float(np.exp(m.params[n])-1),lo=float(np.exp(m.params[n]-1.96*m.bse[n])-1),hi=float(np.exp(m.params[n]+1.96*m.bse[n])-1)) for n in names}, float(df.y[ok].mean())
COV=(pd.Timestamp('2020-01-20'),pd.Timestamp('2020-06-30'))
WB=[('h-3',-21,-15),('h-2',-14,-8),('h-1',-7,-1),('h0',0,6),('h+1',7,13),('h+2',14,20),('h+3',21,27),('h+4',28,34),('h+5',35,41),('h+6',42,48),('h+7',49,55),('h+8',56,62)]
S={g:series(g) for g in GR}
out={}
# Ay Yeni Yili
out['lny']={}
for g in('Çin','Doğu ve Güneydoğu Asya (Çin dışı, Ay Yeni Yılı kutlayan)','Japonya','Dünya'):
    out['lny'][g]={}
    for v in('portcalls','portcalls_container','thr'):
        e,mu=event(S[g][v],H['lny'],-21,35,excl=COV); w,_=event(S[g][v],H['lny'],-21,62,bins=WB,excl=COV)
        out['lny'][g][v]=dict(daily={k:[x['eff'],x['lo'],x['hi']] for k,x in e.items()},weekly={k:[x['eff'],x['lo'],x['hi']] for k,x in w.items()},mean=mu)
# Cin: yil yil (h0 ve h-1) 
out['lny_year']={}
for h in H['lny_all']:
    w,_=event(S['Çin']['portcalls'],pd.DatetimeIndex([h]),-21,34,bins=WB[:8],excl=None if h.year==2020 else COV); out['lny_year'][int(h.year)]={k:[x['eff'],x['lo'],x['hi']] for k,x in w.items()}
# kayip: Cin'de -7..+20 arasindaki beklenen ugraktan sapma (yilda)
e=out['lny']['Çin']['portcalls']; mu=e['mean']; out['lny_loss']=dict(calls=float(sum(e['daily'][str(k)][0]*mu for k in range(-7,21))),days_equiv=float(sum(e['daily'][str(k)][0] for k in range(-7,21))),min_day=min(e['daily'],key=lambda k:e['daily'][k][0]),min_eff=min(x[0] for x in e['daily'].values()))
e=out['lny']['Çin']['thr']; out['lny_loss']['days_equiv_thr']=float(sum(e['daily'][str(k)][0] for k in range(-7,21))); out['lny_loss']['mean_thr']=e['mean']
# yayilma: varis bolgelerinde konteyner ithalati (haftalik, +8 haftaya kadar)
out['spill']={}
for g in('Kuzey Amerika','Avrupa','Latin Amerika','Körfez ülkeleri'):
    out['spill'][g]={}
    for v in('import_container','portcalls_container','export_container'):
        w,_=event(S[g][v],H['lny'],-21,62,bins=WB,excl=COV); out['spill'][g][v]={k:[x['eff'],x['lo'],x['hi']] for k,x in w.items()}
# Cin ihracat/ithalat konteyner
out['lny_cn_dir']={v:{k:[x['eff'],x['lo'],x['hi']] for k,x in event(S['Çin'][v],H['lny'],-21,62,bins=WB,excl=COV)[0].items()} for v in('export_container','import_container','export','import')}
# Noel / Yilbasi
XB=[('-14..-8',-14,-8),('-7..-1',-7,-1),('25 Ara',0,0),('26-31 Ara',1,6),('1 Oca',7,7),('2-8 Oca',8,14),('9-15 Oca',15,21)]
out['xmas']={}
for g in('Avrupa','Kuzey Amerika','Latin Amerika','Çin','Japonya','Körfez ülkeleri','Dünya'):
    out['xmas'][g]={v:{k:[x['eff'],x['lo'],x['hi']] for k,x in event(S[g][v],H['xmas'],-14,21,bins=XB,excl=COV)[0].items()} for v in('portcalls','thr')}
e,mu=event(S['Avrupa']['portcalls'],H['xmas'],-14,21,excl=COV); out['xmas_daily_eu']={k:[x['eff'],x['lo'],x['hi']] for k,x in e.items()}
e,mu=event(S['Dünya']['portcalls'],H['xmas'],-14,21,excl=COV); out['xmas_daily_world']={k:[x['eff'],x['lo'],x['hi']] for k,x in e.items()}; out['xmas_world_days']=float(sum(x['eff'] for x in e.values())); out['world_mean']=mu
# Bayramlar
EB=[('-14..-8',-14,-8),('-7..-1',-7,-1),('0..3',0,3),('4..7',4,7),('8..14',8,14)]
out['eid']={}
for hn in('fitr','adha'):
    out['eid'][hn]={}
    for g in('Körfez ülkeleri','Diğer Müslüman çoğunluklu ülkeler','Avrupa','Dünya'):
        out['eid'][hn][g]={v:{k:[x['eff'],x['lo'],x['hi']] for k,x in event(S[g][v],H[hn],-14,14,bins=EB,excl=COV)[0].items()} for v in('portcalls','thr')}
# ulke bazinda bayram etkisi (0..3 gun)
out['eid_country']={}
for iso in ['SAU','ARE','EGY','TUR','IDN','PAK','BGD','MYS','MAR','QAT','OMN','KWT','IRN','DZA']:
    y=d[d.iso==iso].groupby('date').portcalls.sum()
    try:
        a=event(y,H['fitr'],-14,14,bins=EB,excl=COV)[0]['0..3']; b=event(y,H['adha'],-14,14,bins=EB,excl=COV)[0]['0..3']; out['eid_country'][iso]=dict(name=str(P[P.ISO3==iso].country.iloc[0]),fitr=[a['eff'],a['lo'],a['hi']],adha=[b['eff'],b['lo'],b['hi']],mean=float(y.mean()))
    except Exception as ex: print(iso,ex)
# Cin Ulusal Gunu
GB=[('-7..-1',-7,-1),('1-7 Eki',0,6),('8-14 Eki',7,13)]
out['golden']={v:{k:[x['eff'],x['lo'],x['hi']] for k,x in event(S['Çin'][v],H['golden'],-7,13,bins=GB,excl=COV)[0].items()} for v in('portcalls','portcalls_container','thr')}
# ---- BAE hafta sonu degisikligi (1 Oca 2022: Cum-Cmt -> Cmt-Paz)
def dowprof(y):
    df=pd.DataFrame({'y':y}); df['dow']=df.index.dayofweek; df['wk']=(df.index-pd.Timestamp('2018-12-31')).days//7; w=df.pivot(index='wk',columns='dow',values='y').dropna(); r=w.div(w.mean(axis=1),axis=0); return r
out['uae']={}
rng=np.random.default_rng(5)
for iso,nm in(('ARE','BAE'),('SAU','Suudi Arabistan'),('QAT','Katar'),('OMN','Umman'),('KWT','Kuveyt'),('EGY','Mısır')):
    y=d[d.iso==iso].groupby('date').portcalls.sum(); r=dowprof(y); wkdate=pd.Timestamp('2018-12-31')+pd.to_timedelta(r.index*7,'D'); pre=r[(wkdate<'2022-01-01')&~((wkdate>='2020-01-20')&(wkdate<='2020-06-30'))]; post=r[wkdate>='2022-01-01']
    dif=post.mean()-pre.mean(); bs=np.array([post.iloc[rng.integers(0,len(post),len(post))].mean().values-pre.iloc[rng.integers(0,len(pre),len(pre))].mean().values for _ in range(1000)])
    out['uae'][nm]=dict(pre=[float(x) for x in pre.mean()],post=[float(x) for x in post.mean()],diff=[float(x) for x in dif],lo=[float(x) for x in np.percentile(bs,2.5,axis=0)],hi=[float(x) for x in np.percentile(bs,97.5,axis=0)],n_pre=int(len(pre)),n_post=int(len(post)),mean=float(y.mean()))
# fark-in-farklar: BAE (Cum - Paz) degisimi eksi komsular ortalamasi
def fs(nm): o=out['uae'][nm]; return (o['post'][4]-o['post'][6])-(o['pre'][4]-o['pre'][6])
out['uae_did']=dict(uae=float(fs('BAE')),others={k:float(fs(k)) for k in out['uae'] if k!='BAE'})
R['hol']=out; json.dump(R,open('h/out/results.json','w'),ensure_ascii=False,indent=1)
f=lambda x:[round(100*v,1) for v in x]
print('LNY weekly'); [print(g,v,{k:f(x) for k,x in out['lny'][g][v]['weekly'].items()}) for g in out['lny'] for v in('portcalls','thr')]
print('loss',out['lny_loss']); print('year',{y:{k:f(x) for k,x in v.items() if k in('h-1','h0','h+1','h+2')} for y,v in out['lny_year'].items()})
print('spill'); [print(g,v,{k:f(x) for k,x in out['spill'][g][v].items()}) for g in out['spill'] for v in out['spill'][g]]
print('cn dir',{v:{k:f(x) for k,x in o.items()} for v,o in out['lny_cn_dir'].items()})
print('xmas',{g:{k:f(x) for k,x in o['portcalls'].items()} for g,o in out['xmas'].items()}); print(out['xmas_world_days'],out['world_mean'])
print('eid',{h:{g:{k:f(x) for k,x in o['portcalls'].items()} for g,o in v.items()} for h,v in out['eid'].items()})
print('eidc',{k:(f(v['fitr']),f(v['adha'])) for k,v in out['eid_country'].items()}); print('golden',{v:{k:f(x) for k,x in o.items()} for v,o in out['golden'].items()})
print('uae',{k:([round(x,3) for x in v['pre']],[round(x,3) for x in v['post']]) for k,v in out['uae'].items()}, out['uae_did'])
'''

# ==============================================================================
# h/h_extra.py
# ==============================================================================
MODULES['h/h_extra.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
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['iso']=d.portid.map(P.ISO3); d['cont']=d.portid.map(P.continent); d['dow']=d.date.dt.dayofweek
R=json.load(open('h/out/results.json')); X={}
G=d.groupby('date')[['portcalls','thr','portcalls_container','import_container']].sum()
COV=(pd.Timestamp('2020-01-20'),pd.Timestamp('2020-06-30'))
# ---- 1. hafta sonu oraninda egilim (haftalik seri)
def weekly_wr(s):
    df=pd.DataFrame({'y':s}); df['dow']=df.index.dayofweek; df['wk']=(df.index-pd.Timestamp('2018-12-31')).days//7; w=df.pivot(index='wk',columns='dow',values='y').dropna()
    wr=w[[5,6]].mean(axis=1)/w[[0,1,2,3,4]].mean(axis=1); wr.index=pd.Timestamp('2018-12-31')+pd.to_timedelta(wr.index*7,'D'); return wr
wr=weekly_wr(G.portcalls); t=(wr.index-wr.index[0]).days/365.25
m=sm.OLS(wr.values,sm.add_constant(t)).fit(cov_type='HAC',cov_kwds={'maxlags':8}); X['wr_trend']=dict(per_year=float(m.params[1]),se=float(m.bse[1]),p=float(m.pvalues[1]),n=int(len(wr)),first=float(wr[:52].mean()),last=float(wr[-52:].mean()))
wrt=weekly_wr(G.thr); m=sm.OLS(wrt.values,sm.add_constant(t)).fit(cov_type='HAC',cov_kwds={'maxlags':8}); X['wr_trend_thr']=dict(per_year=float(m.params[1]),se=float(m.bse[1]),p=float(m.pvalues[1]))
# gelir grubuna gore egilim
wb=pd.read_csv('data/wb_country.csv').set_index('iso3'); d['inc']=d.iso.map(wb.income)
X['wr_trend_inc']={}
for k,g in d.groupby('inc'):
    s=g.groupby('date').portcalls.sum(); w=weekly_wr(s); tt=(w.index-w.index[0]).days/365.25; m=sm.OLS(w.values,sm.add_constant(tt)).fit(cov_type='HAC',cov_kwds={'maxlags':8})
    X['wr_trend_inc'][k]=dict(mean=float(w.mean()),per_year=float(m.params[1]),se=float(m.bse[1]),p=float(m.pvalues[1]),share=float(g.portcalls.sum()/d.portcalls.sum()))
# ---- 2. ugrak basina yuk: hafta ici / hafta sonu
gd=d.groupby('dow')[['portcalls','thr']].sum(); cpc=gd.thr/gd.portcalls; X['cargo_per_call']=dict(idx=[float(x) for x in (cpc/cpc.mean()).values],weekend_vs_weekday=float(cpc[[5,6]].mean()/cpc[[0,1,2,3,4]].mean()),mean=float(d.thr.sum()/d.portcalls.sum()))
# ---- 3. cuma hafta sonu olan ulkeler ve Israil
FRI=['SAU','QAT','KWT','OMN','BHR','EGY','DZA','IRN','IRQ','JOR','LBY','BGD','YEM','SDN']
def prof(mask):
    s=d[mask].groupby('date').portcalls.sum(); df=pd.DataFrame({'y':s}); df['dow']=df.index.dayofweek; df['wk']=(df.index-pd.Timestamp('2018-12-31')).days//7; w=df.pivot(index='wk',columns='dow',values='y').dropna(); w=w[w.mean(axis=1)>0]; r=w.div(w.mean(axis=1),axis=0)
    rng=np.random.default_rng(2); bs=np.array([r.iloc[rng.integers(0,len(r),len(r))].mean().values for _ in range(500)])
    return dict(idx=[float(x) for x in r.mean()],lo=[float(x) for x in np.percentile(bs,2.5,axis=0)],hi=[float(x) for x in np.percentile(bs,97.5,axis=0)],calls=float(s.sum()),nports=int(d[mask].portid.nunique()))
X['fri_group']=prof(d.iso.isin(FRI)); X['israel']=prof(d.iso=='ISR'); X['sun_group']=prof(~d.iso.isin(FRI+['ISR','ARE']))
# ---- 4. liman turu: tepe gunu dagilimi ve profilin zaman icinde kaliciligi
Q=pd.read_csv('h/out/port_weekly.csv').set_index('portid'); S=Q[Q.tot>=1000]
X['peak_by_type']={k:{'Pzt':float((g.peak==0).mean()),'hafta_sonu':float(g.peak.isin([5,6]).mean()),'trough_we':float(g.trough.isin([5,6]).mean())} for k,g in S.groupby('type')}
a=d[d.date<'2023-01-01'].groupby(['portid','dow']).portcalls.sum().unstack(); b=d[d.date>='2023-01-01'].groupby(['portid','dow']).portcalls.sum().unstack()
a=a.div(a.mean(axis=1),axis=0).loc[S.index]; b=b.div(b.mean(axis=1),axis=0).loc[S.index]
cc=pd.Series([np.corrcoef(a.loc[i],b.loc[i])[0,1] for i in S.index],index=S.index)
X['persist']=dict(med=float(cc.median()),p25=float(cc.quantile(.25)),p75=float(cc.quantile(.75)),by_type={k:float(cc[g.index].median()) for k,g in S.groupby('type')},
   wr_corr=float(np.corrcoef((a[[5,6]].mean(axis=1)/a[[0,1,2,3,4]].mean(axis=1)),(b[[5,6]].mean(axis=1)/b[[0,1,2,3,4]].mean(axis=1)))[0,1]))
# ---- 5. varyans ayristirmasi: dunya gunluk ugraklari
H=R['hol']
df=pd.DataFrame({'y':G.portcalls}); df=df[~((df.index>=COV[0])&(df.index<=COV[1]))]
df['ly']=np.log(df.y); df['doy']=df.index.dayofyear; df['t']=(df.index-df.index[0]).days/365.25
B=pd.get_dummies(df.index.strftime('%Y'),prefix='y',drop_first=True).astype(float).set_index(df.index)
for j in(1,2,3): B['s%d'%j]=np.sin(2*np.pi*j*df.doy/365.25); B['c%d'%j]=np.cos(2*np.pi*j*df.doy/365.25)
B['t']=df.t; B=sm.add_constant(B)
DW=pd.get_dummies(df.index.dayofweek,prefix='d',drop_first=True).astype(float).set_index(df.index)
def hol_dummies(idx):
    Hd=pd.DataFrame(index=idx)
    lny=pd.to_datetime(['2019-02-05','2020-01-25','2021-02-12','2022-02-01','2023-01-22','2024-02-10','2025-01-29','2026-02-17'])
    fitr=pd.to_datetime(['2019-06-04','2020-05-24','2021-05-13','2022-05-02','2023-04-21','2024-04-10','2025-03-30','2026-03-20'])
    adha=pd.to_datetime(['2019-08-11','2020-07-31','2021-07-20','2022-07-09','2023-06-28','2024-06-16','2025-06-06','2026-05-27'])
    def win(dates,a,b):
        v=np.zeros(len(idx))
        for h in dates:
            k=(idx-h).days; v[(k>=a)&(k<=b)]=1
        return v
    Hd['xmas']=((idx.month==12)&(idx.day==25)).astype(float); Hd['xmas2']=((idx.month==12)&(idx.day>=26)).astype(float); Hd['xmas0']=((idx.month==12)&(idx.day==24)).astype(float); Hd['ny']=((idx.month==1)&(idx.day==1)).astype(float); Hd['ny2']=((idx.month==1)&(idx.day>=2)&(idx.day<=3)).astype(float)
    Hd['lny_m1']=win(lny,-3,-1); Hd['lny_0']=win(lny,0,6); Hd['lny_1']=win(lny,7,13); Hd['fitr']=win(fitr,0,3); Hd['adha']=win(adha,0,3); Hd['may1']=((idx.month==5)&(idx.day==1)).astype(float); Hd['gw']=((idx.month==10)&(idx.day<=7)).astype(float)
    return Hd
HD=hol_dummies(df.index)
def r2(Xm): return sm.OLS(df.ly,Xm).fit().rsquared
r_b=r2(B); r_bd=r2(B.join(DW)); r_bh=r2(B.join(HD)); r_all=r2(B.join(DW).join(HD))
res=sm.OLS(df.ly,B).fit().resid; tot=res.var()
X['variance']=dict(r2_base=float(r_b),r2_dow=float(r_bd),r2_hol=float(r_bh),r2_all=float(r_all),share_dow=float((r_bd-r_b)/(1-r_b)),share_hol=float((r_all-r_bd)/(1-r_b)),share_cal=float((r_all-r_b)/(1-r_b)),resid_sd=float(np.sqrt(tot)),n=int(len(df)))
mm=sm.OLS(df.ly,B.join(DW).join(HD)).fit(cov_type='HAC',cov_kwds={'maxlags':14}); X['variance']['hol_coef']={k:[float(np.exp(mm.params[k])-1),float(mm.bse[k]),float(mm.pvalues[k])] for k in HD.columns}
# kayip gun esdegeri (dunya, yilda): hafta sonu yok (yeniden dagilim), tatiller toplam
X['variance']['hol_days_per_year']=float(sum((np.exp(mm.params[k])-1)*HD[k].sum() for k in HD.columns)/ (len(df)/365.25))
# ---- 6. Ay Yeni Yili: yerel taban cizgisi (mevsim modelinden bagimsiz), yil yil
CN=d[d.iso=='CHN'].groupby('date')[['portcalls','export_container','thr']].sum()
def local(y,h,win=(0,13),base=((-49,-22),(36,63))):
    k=(y.index-h).days; s=pd.DataFrame({'y':y.values,'k':k,'dow':y.index.dayofweek})
    pre=s[(s.k>=base[0][0])&(s.k<=base[0][1])]; post=s[(s.k>=base[1][0])&(s.k<=base[1][1])]
    if len(post)<20 and len(pre)<20: return None
    # ayni gun-hafta icin taban: on ve arka pencere ortalamasinin ortalamasi (biri eksikse yalniz oteki)
    if len(pre)<20: bm=post.groupby('dow').y.mean()
    elif len(post)<20: bm=pre.groupby('dow').y.mean()
    else: bm=(pre.groupby('dow').y.mean()+post.groupby('dow').y.mean())/2
    w=s[(s.k>=win[0])&(s.k<=win[1])]; return float(w.y.sum()/w.dow.map(bm).sum()-1)
lny_all=pd.to_datetime(['2019-02-05','2020-01-25','2021-02-12','2022-02-01','2023-01-22','2024-02-10','2025-01-29','2026-02-17'])
X['lny_local']={int(h.year):dict(calls_w0=local(CN.portcalls,h,(0,6)),calls_w1=local(CN.portcalls,h,(7,13)),calls_2w=local(CN.portcalls,h,(0,13)),calls_pre=local(CN.portcalls,h,(-7,-1)),exp_cont_2w=local(CN.export_container,h,(7,20)),thr_2w=local(CN.thr,h,(0,13))) for h in lny_all}
v=[x['calls_2w'] for y,x in X['lny_local'].items() if y!=2020]; X['lny_local_mean']=dict(calls_2w=float(np.mean(v)),sd=float(np.std(v,ddof=1)),n=len(v),exp_cont=float(np.mean([x['exp_cont_2w'] for y,x in X['lny_local'].items() if y!=2020])),calls_w0=float(np.mean([x['calls_w0'] for y,x in X['lny_local'].items() if y!=2020])),calls_w1=float(np.mean([x['calls_w1'] for y,x in X['lny_local'].items() if y!=2020])))
# ham kontrol 2026 ve 2025: LNY cevresinde 7 gunluk ortalamalar
X['lny_raw']={int(h.year):[float(CN.portcalls[(CN.index>=h+pd.Timedelta(days=a))&(CN.index<h+pd.Timedelta(days=a+7))].mean()) for a in(-28,-21,-14,-7,0,7,14,21,28)] for h in lny_all}
# ---- 7. yayilma: ABD Bati Kiyisi konteyner ithalati, yerel taban cizgisi ile haftalar
USWC=P[(P.ISO3=='USA')&(P.lon<-115)&(P.lat>30)&(P.lat<50)].index; EU=P[P.continent=='Europe'].index
def spill(ids,col):
    y=d[d.portid.isin(ids)].groupby('date')[col].sum(); out={}
    for wk in range(-2,9):
        v=[local(y,h,(7*wk,7*wk+6),base=((-63,-29),(63,97))) for h in lny_all if h.year!=2020]; v=[x for x in v if x is not None]
        out[str(wk)]=[float(np.mean(v)),float(np.std(v,ddof=1)/np.sqrt(len(v))),len(v)]
    return out
X['spill_uswc']=dict(nports=int(len(USWC)),names=[str(P.portname[i]) for i in USWC][:12],imp=spill(USWC,'import_container'),calls=spill(USWC,'portcalls_container'))
X['spill_eu']=dict(imp=spill(EU,'import_container'))
X['spill_cn_exp']=spill(P[P.ISO3=='CHN'].index,'export_container')
# ---- 8. ek tatiller (Poisson GLM, h_holiday ile ayni model)
def series(mask,col='portcalls'): return d[mask].groupby('date')[col].sum()
def event(y,dates_fn,bins,excl=COV):
    df=pd.DataFrame({'y':y}); df['doy']=df.index.dayofyear
    Xm=pd.get_dummies(df.index.dayofweek,prefix='d',drop_first=True).astype(float).set_index(df.index)
    Xm=Xm.join(pd.get_dummies(df.index.strftime('%Y'),prefix='y',drop_first=True).astype(float).set_index(df.index))
    for j in(1,2,3): Xm['s%d'%j]=np.sin(2*np.pi*j*df.doy/365.25); Xm['c%d'%j]=np.cos(2*np.pi*j*df.doy/365.25)
    Xm['t']=(df.index-df.index[0]).days/365.25
    for nm,v in dates_fn(df.index).items(): Xm['r_'+nm]=v
    Xm=sm.add_constant(Xm); ok=~((df.index>=excl[0])&(df.index<=excl[1]))
    m=sm.GLM(df.y[ok],Xm[ok],family=sm.families.Poisson()).fit(cov_type='HAC',cov_kwds={'maxlags':14})
    return {n[2:]:[float(np.exp(m.params[n])-1),float(np.exp(m.params[n]-1.96*m.bse[n])-1),float(np.exp(m.params[n]+1.96*m.bse[n])-1)] for n in m.params.index if n.startswith('r_')}
def W(dates,a,b):
    def f(idx):
        v=np.zeros(len(idx))
        for h in dates:
            k=(idx-h).days; v[(k>=a)&(k<=b)]=1
        return v
    return f
easter=pd.to_datetime(['2019-04-21','2020-04-12','2021-04-04','2022-04-17','2023-04-09','2024-03-31','2025-04-20','2026-04-05'])
tg=pd.to_datetime(['2019-11-28','2020-11-26','2021-11-25','2022-11-24','2023-11-23','2024-11-28','2025-11-27'])
EUm=d.cont=='Europe'; NAm=d.iso.isin(['USA','CAN']); USm=d.iso=='USA'; JPm=d.iso=='JPN'; CNm=d.iso=='CHN'
X['more']={}
X['more']['easter_eu']=event(series(EUm),lambda idx:{'Kutsal Cuma':W(easter,-2,-2)(idx),'Paskalya Pazarı':W(easter,0,0)(idx),'Paskalya Pazartesisi':W(easter,1,1)(idx),'Cumartesi':W(easter,-1,-1)(idx)},None)
X['more']['easter_na']=event(series(NAm),lambda idx:{'Kutsal Cuma':W(easter,-2,-2)(idx),'Paskalya Pazarı':W(easter,0,0)(idx),'Paskalya Pazartesisi':W(easter,1,1)(idx)},None)
X['more']['thanksgiving_us']=event(series(USm),lambda idx:{'Şükran Günü':W(tg,0,0)(idx),'ertesi gün':W(tg,1,1)(idx)},None)
X['more']['july4_us']=event(series(USm),lambda idx:{'4 Temmuz':((idx.month==7)&(idx.day==4)).astype(float)},None)
X['more']['may1_eu']=event(series(EUm),lambda idx:{'1 Mayıs':((idx.month==5)&(idx.day==1)).astype(float)},None)
X['more']['may1_cn']=event(series(CNm),lambda idx:{'1-5 Mayıs':((idx.month==5)&(idx.day<=5)).astype(float)},None)
X['more']['gw_jp']=event(series(JPm),lambda idx:{'29 Nis-5 May':(((idx.month==4)&(idx.day>=29))|((idx.month==5)&(idx.day<=5))).astype(float),'3-5 May':((idx.month==5)&(idx.day>=3)&(idx.day<=5)).astype(float)},None)
X['more']['obon_jp']=event(series(JPm),lambda idx:{'13-16 Ağu':((idx.month==8)&(idx.day>=13)&(idx.day<=16)).astype(float)},None)
X['more']['ny_jp']=event(series(JPm),lambda idx:{'29-31 Ara':((idx.month==12)&(idx.day>=29)).astype(float),'1-3 Oca':((idx.month==1)&(idx.day<=3)).astype(float),'4-7 Oca':((idx.month==1)&(idx.day>=4)&(idx.day<=7)).astype(float)},None)
X['more']['israel_kippur']=None
R['x']=X; json.dump(R,open('h/out/results.json','w'),ensure_ascii=False,indent=1)
f=lambda v:[round(100*x,1) for x in v]
for k in('wr_trend','wr_trend_thr','wr_trend_inc','cargo_per_call','peak_by_type','persist','variance','lny_local','lny_local_mean','lny_raw'): print(k,json.dumps(X[k],ensure_ascii=False)[:1500]); print()
print('fri',[round(x,3) for x in X['fri_group']['idx']],X['fri_group']['calls'],'isr',[round(x,3) for x in X['israel']['idx']],X['israel']['calls'],'sun',[round(x,3) for x in X['sun_group']['idx']])
print('uswc',X['spill_uswc']['nports'],X['spill_uswc']['names'],{k:[round(100*v[0],1),round(100*v[1],1)] for k,v in X['spill_uswc']['imp'].items()}); print('uswc calls',{k:[round(100*v[0],1),round(100*v[1],1)] for k,v in X['spill_uswc']['calls'].items()})
print('eu',{k:[round(100*v[0],1),round(100*v[1],1)] for k,v in X['spill_eu']['imp'].items()}); print('cn exp',{k:[round(100*v[0],1),round(100*v[1],1)] for k,v in X['spill_cn_exp'].items()})
print('more',{k:({a:f(b) for a,b in v.items()} if v else None) for k,v in X['more'].items()})
'''

# ==============================================================================
# h/h_extra2.py
# ==============================================================================
MODULES['h/h_extra2.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['iso']=d.portid.map(P.ISO3)
R=json.load(open('h/out/results.json')); X=R['x']
Q=pd.read_csv('h/out/port_weekly.csv').set_index('portid'); S=Q[Q.tot>=1000]
a=d[d.date<'2023-01-01'].groupby(['portid','dow' if 'dow' in d else d.date.dt.dayofweek]).portcalls.sum().unstack() if False else None
d['dow']=d.date.dt.dayofweek
a=d[d.date<'2023-01-01'].groupby(['portid','dow']).portcalls.sum().unstack().loc[S.index]; b=d[d.date>='2023-01-01'].groupby(['portid','dow']).portcalls.sum().unstack().loc[S.index]
wa=a[[5,6]].mean(axis=1)/a[[0,1,2,3,4]].mean(axis=1); wb_=b[[5,6]].mean(axis=1)/b[[0,1,2,3,4]].mean(axis=1); ok=np.isfinite(wa)&np.isfinite(wb_)&(wa>0)&(wb_>0)
X['persist']['wr_corr']=float(np.corrcoef(np.log(wa[ok]),np.log(wb_[ok]))[0,1]); X['persist']['wr_spearman']=float(pd.Series(wa[ok]).corr(pd.Series(wb_[ok]),method='spearman')); X['persist']['n']=int(ok.sum())
# haftalik toplamlar: tatillerin acikladigi pay
G=d.groupby('date').portcalls.sum(); COV=(pd.Timestamp('2020-01-20'),pd.Timestamp('2020-06-30'))
idx=G.index
lny=pd.to_datetime(['2019-02-05','2020-01-25','2021-02-12','2022-02-01','2023-01-22','2024-02-10','2025-01-29','2026-02-17'])
fitr=pd.to_datetime(['2019-06-04','2020-05-24','2021-05-13','2022-05-02','2023-04-21','2024-04-10','2025-03-30','2026-03-20'])
adha=pd.to_datetime(['2019-08-11','2020-07-31','2021-07-20','2022-07-09','2023-06-28','2024-06-16','2025-06-06','2026-05-27'])
def win(dates,a_,b_):
    v=np.zeros(len(idx))
    for h in dates:
        k=(idx-h).days; v[(k>=a_)&(k<=b_)]=1
    return v
Hd=pd.DataFrame(index=idx); Hd['xmas']=(((idx.month==12)&(idx.day>=24))|((idx.month==1)&(idx.day<=3))).astype(float); Hd['lny']=win(lny,-3,13); Hd['eid']=win(fitr,0,3)+win(adha,0,3); Hd['may1']=((idx.month==5)&(idx.day==1)).astype(float)
wk=(idx-pd.Timestamp('2018-12-31')).days//7
Wy=G.groupby(wk).sum(); Wn=G.groupby(wk).size(); WH=Hd.groupby(wk).sum(); wdate=pd.Timestamp('2018-12-31')+pd.to_timedelta(Wy.index*7,'D')
ok=(Wn==7).values&~((wdate>=COV[0])&(wdate<=COV[1]))
df=pd.DataFrame({'ly':np.log(Wy.values)},index=wdate)[ok]; doy=df.index.dayofyear+3
B=pd.get_dummies(df.index.strftime('%Y'),prefix='y',drop_first=True).astype(float).set_index(df.index)
for j in(1,2,3): B['s%d'%j]=np.sin(2*np.pi*j*doy/365.25); B['c%d'%j]=np.cos(2*np.pi*j*doy/365.25)
B['t']=(df.index-df.index[0]).days/365.25; B=sm.add_constant(B); WHk=WH[ok]; WHk.index=df.index
r0=sm.OLS(df.ly,B).fit(); r1=sm.OLS(df.ly,B.join(WHk)).fit()
X['weekly_var']=dict(n=int(len(df)),r2_base=float(r0.rsquared),r2_hol=float(r1.rsquared),share_hol=float((r1.rsquared-r0.rsquared)/(1-r0.rsquared)),resid_sd_base=float(np.sqrt(r0.resid.var())),resid_sd_hol=float(np.sqrt(r1.resid.var())))
# Cin haftalik
GC=d[d.iso=='CHN'].groupby('date').portcalls.sum(); Wy=GC.groupby(wk).sum(); df=pd.DataFrame({'ly':np.log(Wy.values)},index=wdate)[ok]
r0=sm.OLS(df.ly,B).fit(); r1=sm.OLS(df.ly,B.join(WHk[['lny']])).fit(); X['weekly_var_cn']=dict(r2_base=float(r0.rsquared),r2_hol=float(r1.rsquared),share_hol=float((r1.rsquared-r0.rsquared)/(1-r0.rsquared)),resid_sd_base=float(np.sqrt(r0.resid.var())),resid_sd_hol=float(np.sqrt(r1.resid.var())))
R['x']=X; json.dump(R,open('h/out/results.json','w'),ensure_ascii=False,indent=1)
print(X['persist']); print(X['weekly_var'],X['weekly_var_cn']); print(X['lny_local']['2026'],X['lny_local']['2025'])
'''

# ==============================================================================
# h/h_robust.py
# ==============================================================================
MODULES['h/h_robust.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['iso']=d.portid.map(P.ISO3); d['cont']=d.portid.map(P.continent)
R=json.load(open('h/out/results.json')); RB={}
G=d.groupby('date').portcalls.sum()
def wr(s):
    df=pd.DataFrame({'y':s}); df['dow']=df.index.dayofweek; df['wk']=(df.index-pd.Timestamp('2018-12-31')).days//7; w=df.pivot(index='wk',columns='dow',values='y').dropna(); r=w.div(w.mean(axis=1),axis=0); return float((r[[5,6]].mean(axis=1)/r[[0,1,2,3,4]].mean(axis=1)).mean()),len(r)
RB['wr_all']=wr(G)[0]
RB['wr_ex2020']=wr(G[G.index.year!=2020])[0]
# tatil haftalari haric: Noel-yilbasi (20 Ara-7 Oca) ve LNY +-10 gun
lny=pd.to_datetime(['2019-02-05','2020-01-25','2021-02-12','2022-02-01','2023-01-22','2024-02-10','2025-01-29','2026-02-17'])
m=((G.index.month==12)&(G.index.day>=20))|((G.index.month==1)&(G.index.day<=7))
for h in lny: m|=(abs((G.index-h).days)<=10)
RB['wr_exhol']=wr(G[~m])[0]
# buyuk limanlar (>= 10 ugrak/gun) ve kucukler
tot=d.groupby('portid').portcalls.sum()/d.date.nunique(); big=tot[tot>=10].index; small=tot[tot<1].index
RB['wr_big']=wr(d[d.portid.isin(big)].groupby('date').portcalls.sum())[0]; RB['n_big']=int(len(big)); RB['wr_small']=wr(d[d.portid.isin(small)].groupby('date').portcalls.sum())[0]; RB['n_small']=int(len(small))
# plasebo: tatil tarihleri kaydirilmis
COV=(pd.Timestamp('2020-01-20'),pd.Timestamp('2020-06-30'))
def event(y,dates,a,b):
    df=pd.DataFrame({'y':y}); df['doy']=df.index.dayofyear
    X=pd.get_dummies(df.index.dayofweek,prefix='d',drop_first=True).astype(float).set_index(df.index)
    X=X.join(pd.get_dummies(df.index.strftime('%Y'),prefix='y',drop_first=True).astype(float).set_index(df.index))
    for j in(1,2,3): X['s%d'%j]=np.sin(2*np.pi*j*df.doy/365.25); X['c%d'%j]=np.cos(2*np.pi*j*df.doy/365.25)
    X['t']=(df.index-df.index[0]).days/365.25; v=np.zeros(len(df))
    for h in dates:
        k=(df.index-h).days; v[(k>=a)&(k<=b)]=1
    X['r']=v; X=sm.add_constant(X); ok=~((df.index>=COV[0])&(df.index<=COV[1]))
    mm=sm.GLM(df.y[ok],X[ok],family=sm.families.Poisson()).fit(cov_type='HAC',cov_kwds={'maxlags':14})
    return [float(np.exp(mm.params['r'])-1),float(np.exp(mm.params['r']-1.96*mm.bse['r'])-1),float(np.exp(mm.params['r']+1.96*mm.bse['r'])-1)]
CN=d[d.iso=='CHN'].groupby('date').portcalls.sum(); EU=d[d.cont=='Europe'].groupby('date').portcalls.sum()
l7=lny[lny.year!=2020]
RB['placebo_lny']={str(s):event(CN,l7+pd.Timedelta(days=s),0,13) for s in(-120,-60,60,120,180,240)}
RB['true_lny']=event(CN,l7,0,13)
xm=pd.to_datetime(['%d-12-25'%y for y in range(2019,2026)])
RB['placebo_xmas']={str(s):event(EU,xm+pd.Timedelta(days=s),0,0) for s in(-180,-120,-60,60,120,180)}
RB['true_xmas']=event(EU,xm,0,0)
# gunluk Poisson yerine haftalik profilde yuk agirlikli sonuc ve 2019 vs 2025 karsilastirmasi
RB['wr_year']={int(y):wr(g)[0] for y,g in G.groupby(G.index.year)}
R['rb']=RB; json.dump(R,open('h/out/results.json','w'),ensure_ascii=False,indent=1); print(json.dumps(RB,indent=0))
'''

# ==============================================================================
# h/h_fig.py
# ==============================================================================
MODULES['h/h_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('h/out/results.json')); H=R['hol']; Q=pd.read_csv('h/out/port_weekly.csv'); Q=Q[Q.tot>=1000]
B='#1b4f72'; O='#b9770e'; Rd='#b03a2e'; Gn='#1e8449'; DN=['Pzt','Sal','Çar','Per','Cum','Cmt','Paz']
# S1 kuresel haftalik profil
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1)); x=np.arange(7)
a=ax[0]
for k,c,l in(('Tümü',B,'Gemi uğrağı'),('thr',O,'Yük (ton)')):
    g=R['global'][k]; a.errorbar(x,100*np.array(g['idx']),yerr=[100*(np.array(g['idx'])-np.array(g['lo'])),100*(np.array(g['hi'])-np.array(g['idx']))],fmt='o-',color=c,ms=4,lw=1.2,capsize=2,label=l)
a.axhline(100,color='#999',lw=.6); a.set_xticks(x); a.set_xticklabels(DN); a.set_ylabel('Haftalık ortalamaya göre (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Dünya toplamı',fontsize=9,loc='left')
a=ax[1]
for k,c in(('Konteyner',B),('Kuru dökme',O),('Tanker',Rd),('Genel kargo',Gn),('Ro-Ro','#6c3483')):
    a.plot(x,100*np.array(R['global'][k]['idx']),'o-',color=c,ms=3,lw=1.1,label=k)
a.axhline(100,color='#999',lw=.6); a.set_xticks(x); a.set_xticklabels(DN); a.legend(frameon=False,fontsize=7,ncol=2); a.set_title('(b) Gemi türüne göre uğraklar',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('h/fig/Figure_1.png'); plt.close()
# S2 liman duzeyi dagilim + harita
e=pd.read_csv('data/etopo_30min.csv',skiprows=[1]); e=e[(e.altitude>=0)].iloc[::2]; P=pd.read_csv('data/pw/ports_global.csv').set_index('portid'); Q['lon']=Q.portid.map(P.lon)
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0),gridspec_kw={'width_ratios':[1,1.7]})
a=ax[0]; a.hist(Q.wr.clip(0.4,1.4),bins=40,color='#aed6f1',edgecolor='white',lw=.3); a.axvline(1,color='#999',lw=.7); a.axvline(Q.wr.median(),color=Rd,lw=1,ls='--'); a.set_xlabel('Hafta sonu / hafta içi uğrak oranı'); a.set_ylabel('Liman sayısı'); a.set_title('(a) Limanlar arası dağılım',fontsize=9,loc='left')
a=ax[1]; a.scatter(e.longitude,e.latitude,s=.25,color='#dcdcdc',marker='s',lw=0,rasterized=True); sc=a.scatter(Q.lon,Q.lat,c=Q.wr.clip(0.7,1.1),s=5+np.sqrt(Q.tot)/12,cmap='RdYlBu',vmin=0.7,vmax=1.1,lw=0,alpha=.85); cb=fig.colorbar(sc,ax=a,shrink=.75,pad=.01); cb.set_label('Hafta sonu / hafta içi',fontsize=7.5)
a.set_xlim(-180,180); a.set_ylim(-60,75); a.set_xticks([]); a.set_yticks([]); a.set_title('(b) Coğrafi dağılım',fontsize=9,loc='left')
for s in a.spines.values(): s.set_visible(False)
fig.tight_layout(); fig.savefig('h/fig/Figure_2.png'); plt.close()
# S3 Ay Yeni Yili: Cin gunluk + konteyner ihracati haftalik
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]; dly=H['lny']['Çin']['portcalls']['daily']; k=sorted(int(i) for i in dly); v=np.array([dly[str(i)] for i in k])*100; a.fill_between(k,v[:,1],v[:,2],color='#aed6f1',alpha=.6,lw=0); a.plot(k,v[:,0],color=B,lw=1.3,label='Çin')
dl2=H['lny']['Doğu ve Güneydoğu Asya (Çin dışı, Ay Yeni Yılı kutlayan)']['portcalls']['daily']; v2=np.array([dl2[str(i)] for i in k])*100; a.plot(k,v2[:,0],color=O,lw=1.1,label='Diğer kutlayan ülkeler')
a.axhline(0,color='#999',lw=.6); a.axvline(0,color='#999',lw=.6,ls=':'); a.set_xlabel('Ay Yeni Yılı gününe göre gün'); a.set_ylabel('Beklenen düzeyden sapma (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Gemi uğrakları',fontsize=9,loc='left')
a=ax[1]; X=R['x']; xx=np.arange(-2,9)
for key,c,l in((X['spill_cn_exp'],Rd,'Çin: konteyner ihracatı'),(X['spill_uswc']['imp'],B,'ABD Batı Kıyısı: konteyner ithalatı')):
    v=np.array([key[str(w)] for w in xx]); a.errorbar(xx,100*v[:,0],yerr=196*v[:,1],fmt='o-',color=c,ms=3.5,lw=1.1,capsize=2,label=l)
a.axhline(0,color='#999',lw=.6); a.set_xlabel('Ay Yeni Yılına göre hafta'); a.set_ylabel('Yerel taban düzeyinden sapma (%)'); a.set_ylim(-50,55); a.legend(frameon=False,fontsize=7.5,loc='upper right'); a.set_title('(b) Konteyner yükü ve yayılma',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('h/fig/Figure_3.png'); plt.close()
# S4 Noel/Yilbasi gunluk: Avrupa ve dunya; + bolgeler cubuk
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
a=ax[0]
for key,c,l in(('xmas_daily_eu',B,'Avrupa'),('xmas_daily_world',O,'Dünya')):
    dly=H[key]; k=sorted(int(i) for i in dly); v=np.array([dly[str(i)] for i in k])*100; a.plot(k,v[:,0],color=c,lw=1.2,label=l); a.fill_between(k,v[:,1],v[:,2],color=c,alpha=.15,lw=0)
a.axhline(0,color='#999',lw=.6); a.axvline(0,color='#999',lw=.6,ls=':'); a.axvline(7,color='#999',lw=.6,ls=':'); a.set_xlabel('25 Aralık gününe göre gün'); a.set_ylabel('Beklenen düzeyden sapma (%)'); a.legend(frameon=False,fontsize=7.5); a.set_title('(a) Noel ve yılbaşı çevresinde uğraklar',fontsize=9,loc='left')
a=ax[1]; gs=['Avrupa','Kuzey Amerika','Latin Amerika','Japonya','Çin','Körfez ülkeleri']; y=np.arange(len(gs))[::-1]
for i,(kk,c,l) in enumerate((('25 Ara',Rd,'25 Aralık'),('1 Oca',B,'1 Ocak'))):
    v=np.array([H['xmas'][g]['portcalls'][kk] for g in gs])*100; a.barh(y+0.2-i*0.4,v[:,0],height=.38,color=c,label=l,xerr=[v[:,0]-v[:,1],v[:,2]-v[:,0]],error_kw=dict(lw=.6,capsize=1.5))
a.set_yticks(y); a.set_yticklabels(gs,fontsize=7.5); a.axvline(0,color='#999',lw=.6); a.set_xlabel('Beklenen düzeyden sapma (%)'); a.legend(frameon=False,fontsize=7.5,loc='lower left'); a.set_title('(b) Bölgelere göre',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('h/fig/Figure_4.png'); plt.close()
# S5 Bayramlar ulke bazinda
C=H['eid_country']; ks=sorted(C,key=lambda k:C[k]['fitr'][0]); y=np.arange(len(ks))[::-1]
TR={'SAU':'Suudi Arabistan','ARE':'BAE','EGY':'Mısır','TUR':'Türkiye','IDN':'Endonezya','PAK':'Pakistan','BGD':'Bangladeş','MYS':'Malezya','MAR':'Fas','QAT':'Katar','OMN':'Umman','KWT':'Kuveyt','IRN':'İran','DZA':'Cezayir'}
fig,ax=plt.subplots(figsize=(6.2,3.6))
for i,(kk,c,l) in enumerate((('fitr',B,'Ramazan Bayramı'),('adha',O,'Kurban Bayramı'))):
    v=np.array([C[k][kk] for k in ks])*100; ax.barh(y+0.2-i*0.4,v[:,0],height=.38,color=c,label=l,xerr=[v[:,0]-v[:,1],v[:,2]-v[:,0]],error_kw=dict(lw=.6,capsize=1.5))
ax.set_yticks(y); ax.set_yticklabels([TR[k] for k in ks],fontsize=8); ax.axvline(0,color='#999',lw=.6); ax.set_xlabel('Bayramın ilk dört gününde beklenen düzeyden sapma (%)'); ax.legend(frameon=False,fontsize=8,loc='lower left')
fig.tight_layout(); fig.savefig('h/fig/Figure_5.png'); plt.close()

# S6 yerel hafta sonu
X=R['x']; fig,ax=plt.subplots(1,2,figsize=(7.4,3.1)); x=np.arange(7)
a=ax[0]
for key,c,l in(('sun_group',B,'Cumartesi–pazar tatili olan ülkeler'),('fri_group',O,'Cuma tatili olan ülkeler'),('israel',Rd,'İsrail')):
    g=X[key]; v=100*np.array(g['idx']); a.errorbar(x,v,yerr=[v-100*np.array(g['lo']),100*np.array(g['hi'])-v],fmt='o-',color=c,ms=3.5,lw=1.1,capsize=2,label=l)
a.axhline(100,color='#999',lw=.6); a.set_xticks(x); a.set_xticklabels(DN); a.set_ylabel('Haftalık ortalamaya göre (%)'); a.legend(frameon=False,fontsize=7,loc='lower left'); a.set_title('(a) Yerel hafta sonuna göre',fontsize=9,loc='left')
a=ax[1]; u=H['uae']['BAE']; a.plot(x,100*np.array(u['pre']),'o-',color=B,ms=3.5,lw=1.1,label='2019–2021 (cuma–cumartesi)'); a.plot(x,100*np.array(u['post']),'s-',color=Rd,ms=3.5,lw=1.1,label='2022–2026 (cumartesi–pazar)')
a.axhline(100,color='#999',lw=.6); a.set_xticks(x); a.set_xticklabels(DN); a.legend(frameon=False,fontsize=7,loc='lower left'); a.set_title('(b) Birleşik Arap Emirlikleri',fontsize=9,loc='left')
fig.tight_layout(); fig.savefig('h/fig/Figure_6.png'); plt.close()
'''

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