# -*- coding: utf-8 -*-
# ------------------------------------------------------------------------------
# Küresel su ürünleri ticaretinde ihracatçı çöküşlerinden sonra kaybedilen arzın yerine konması
# Replacement of lost supply after exporter collapses in the global seafood trade
#
# GOMD Yayınlar · Makale 29 — https://gomd.pirireis.edu.tr/publications.php?id=rp29
#
# Bu dosya makalenin bütün modüllerini içerir. Replikasyon paketi klasöründe çalıştırıldığında modülleri yerlerine yazar ve sırayla, her birini ayrı süreçte çalıştırır;
# sonuçları `paper/out/` ve şekilleri `paper/fig/` klasörüne yazar.
# This file contains all modules of the paper. Run in the replication-package folder, it writes the modules to their places and runs them in order, each in its own process;
# it writes results to `paper/out/` and figures to `paper/fig/`.
#
# Kurulum / Setup:  pip install pandas numpy matplotlib pyarrow pyfixest
# Çalıştırma / Run:  python rp29_seafood_supplier_collapse_replacement.py
#
# Trade data: BACI, CEPII (HS 1992, version 202601), distributed under the Etalab Open Licence 2.0; Gaulier, G. and Zignago, S. (2010) BACI: International Trade Database at the Product-Level, CEPII Working Paper 2010-23, https://www.cepii.fr/CEPII/en/bdd_modele/bdd_modele_item.asp?id=37 . The package contains the seafood subset (105 product codes, 2,724,451 flows) extracted from the full database (269,894,500 flows). GDP per capita: World Bank, World Development Indicators (NY.GDP.PCAP.KD), CC BY 4.0.
# 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 = ['baci_extract.py', 'b1_events.py', 'b4_main.py', 'run_alt.py', 'b5_desc.py', 'b6_fig.py']
SKIP_IF_EXISTS = {'baci_extract.py': 'baci_out/seafood.parquet'}   # hazırlık modülü: çıktısı pakette varsa atlanır / preparation module: skipped if its output is in the package
DIRS = ['baci_out', 'paper', 'paper/out', 'paper/fig']

# ==============================================================================
# baci_extract.py
# ==============================================================================
MODULES['baci_extract.py'] = r'''
# BACI HS92 zip dosyasindan su urunleri kodlarini ayiklar (ham dosya: https://www.cepii.fr/DATA_DOWNLOAD/baci/data/BACI_HS92_V202601.zip, 2,4 GB)
# Kullanim: python baci_extract.py [ZIP] [CIKTI_KLASORU]
import zipfile,sys,io,csv,gzip,os
ZIP=sys.argv[1] if len(sys.argv)>1 else 'baci/BACI_HS92_V202601.zip'; out=sys.argv[2] if len(sys.argv)>2 else 'baci_out'; os.makedirs(out,exist_ok=True)
import contextlib
LOG=open(out+'/extract.log','w')
z=zipfile.ZipFile(ZIP)
pref=("03","1604","1605","1504","230120")
tot={}
with gzip.open(out+"/seafood.csv.gz","wt") as g:
    g.write("t,i,j,k,v,q\n")
    for inf in z.infolist():
        if not inf.filename.startswith("BACI_HS92_Y"):
            open(out+"/"+inf.filename,"wb").write(z.read(inf)); continue
        n=0;m=0;allv=0.0
        with z.open(inf) as f:
            tf=io.TextIOWrapper(f,encoding="utf-8"); hdr=tf.readline()
            for line in tf:
                n+=1
                p=line.rstrip("\n").split(",")
                k=p[3].strip().zfill(6)
                try: allv+=float(p[4])
                except: pass
                if k.startswith(pref):
                    m+=1; g.write(",".join(x.strip() for x in p[:3])+","+k+","+p[4].strip()+","+p[5].strip()+"\n")
        print(inf.filename,hdr.strip(),n,m,allv,flush=True); LOG.write('%s %s %d %d %s\n'%(inf.filename,hdr.strip(),n,m,allv)); LOG.flush()
LOG.close()
import pandas as pd
d=pd.read_csv(out+'/seafood.csv.gz',dtype={'k':str},na_values=['NA','']); d.to_parquet(out+'/seafood.parquet'); print('su urunleri akimi',len(d))
'''

# ==============================================================================
# b1_events.py
# ==============================================================================
MODULES['b1_events.py'] = r'''
# Tedarikci cokusu olaylari ve ithalatci-urun paneli
import pandas as pd, numpy as np, json, sys
D='baci_out/'
d=pd.read_parquet(D+'seafood.parquet'); d=d[d.q.notna()&(d.q>0)]
MAP={56:58,442:58,710:711,72:711,516:711,426:711,748:711,688:891,499:891,729:736,728:736,531:530,534:530,535:530}
d['i']=d.i.replace(MAP); d['j']=d.j.replace(MAP)
d=d[(d.i!=d.j)&(d.t<=2023)&(~d.i.isin([849,175]))&(~d.j.isin([849,175]))]
DIG=int(sys.argv[4]) if len(sys.argv)>4 else 6
d['k']=d.k.str[:DIG]
d=d.groupby(['t','i','j','k'],as_index=False).agg(v=('v','sum'),q=('q','sum'))
DROP=float(sys.argv[1]) if len(sys.argv)>1 else 0.5; TAG=sys.argv[2] if len(sys.argv)>2 else ''
# ihracatci-urun-yil
X=d.groupby(['i','k','t']).q.sum().rename('X').reset_index()
full=pd.MultiIndex.from_product([X.i.unique(),X.k.unique(),range(1995,2024)],names=['i','k','t'])
Xf=X.set_index(['i','k','t']).reindex(full,fill_value=0).reset_index().sort_values(['i','k','t'])
g=Xf.groupby(['i','k']).X
Xf['B']=(g.shift(1)+g.shift(2)+g.shift(3))/3
Xf['minpre']=np.minimum(np.minimum(g.shift(1),g.shift(2)),g.shift(3))
W=d.groupby(['k','t']).q.sum().rename('W').reset_index().sort_values(['k','t'])
gw=W.groupby('k').W; W['WB']=(gw.shift(1)+gw.shift(2)+gw.shift(3))/3
Xf=Xf.merge(W,on=['k','t'],how='left'); Xf['wshare']=Xf.B/Xf.WB
Xf['ratio']=Xf.X/Xf.B
Xf['col_any']=((Xf.ratio<=1-DROP)&(Xf.B>=500)&(Xf.minpre>0)).astype(int)
gc=Xf.groupby(['i','k']).col_any
Xf['prior']=(gc.shift(1).fillna(0)+gc.shift(2).fillna(0)+gc.shift(3).fillna(0))
for h in range(1,6): Xf['Xr%d'%h]=Xf.groupby(['i','k']).X.shift(-h)/Xf.B
Xf.to_parquet(D+'X_ikt%s.parquet'%TAG)
# ikili akim paneli: 3 yillik taban
b=d.groupby(['i','j','k','t']).agg(q=('q','sum'),v=('v','sum')).reset_index()
# ithalatci-urun toplamlar
M=b.groupby(['j','k','t']).agg(M=('q','sum'),V=('v','sum'),ns=('i','nunique')).reset_index()
fullM=pd.MultiIndex.from_product([M.j.unique(),M.k.unique(),range(1995,2024)],names=['j','k','t'])
Mf=M.set_index(['j','k','t']).reindex(fullM,fill_value=0).reset_index().sort_values(['j','k','t'])
gm=Mf.groupby(['j','k'])
Mf['B']=(gm.M.shift(1)+gm.M.shift(2)+gm.M.shift(3))/3
Mf['BV']=(gm.V.shift(1)+gm.V.shift(2)+gm.V.shift(3))/3
Mf['minpre']=np.minimum(np.minimum(gm.M.shift(1),gm.M.shift(2)),gm.M.shift(3))
for h in range(0,6): Mf['M%d'%h]=gm.M.shift(-h); Mf['V%d'%h]=gm.V.shift(-h)
Mf['Mm4']=gm.M.shift(4); Mf['Mm5']=gm.M.shift(5); Mf['Mm6']=gm.M.shift(6)
Mf=Mf[(Mf.t>=1998)]
# taban donemi ikili paylar: (i,j,k,t) icin onceki 3 yil ortalama akim
keys=[]
for lag in (1,2,3):
    x=b[['i','j','k','t','q']].copy(); x['t']=x.t+lag; keys.append(x)
pb=pd.concat(keys).groupby(['i','j','k','t']).q.sum().div(3).rename('bq').reset_index()
pb=pb[pb.t<=2023]
pb=pb.merge(Mf[['j','k','t','B']],on=['j','k','t']); pb['s']=pb.bq/pb.B
# LOO ihracatci cokusu: i'nin j disindaki ihracati
pb=pb.merge(Xf[['i','k','t','X','B','wshare','minpre','prior']].rename(columns={'X':'Xi','B':'XB','minpre':'Xminpre'}),on=['i','k','t'],how='left')
cur=b[['i','j','k','t','q']].rename(columns={'q':'cq'})
pb=pb.merge(cur,on=['i','j','k','t'],how='left'); pb['cq']=pb.cq.fillna(0)
pb['Xloo']=pb.Xi-pb.cq; pb['XBloo']=pb.XB-pb.bq
pb['ratio_loo']=np.where(pb.XBloo>0,pb.Xloo/pb.XBloo,np.nan)
# olay tanimi: tedarikcinin toplam (j haric) ihracati tabana gore >=DROP dusmus; tedarikci anlamli (taban>=500 t, j haric taban>=250 t, onceki 3 yilin hepsi >0)
ONSET=int(sys.argv[3]) if len(sys.argv)>3 else 1
pb['collapse']=((pb.ratio_loo<=1-DROP)&(pb.Xi<=(1-DROP)*pb.XB)&(pb.XB>=500)&(pb.XBloo>=250)&(pb.Xminpre>0)&((pb.prior==0)|(ONSET==0))).astype(int)
pb['dloo']=np.where(pb.collapse==1,1-pb.ratio_loo.clip(0,1),0.0)
pb['E']=pb.s*pb.dloo            # beklenen kayip payi
pb['Es']=pb.s*pb.collapse       # cokmus tedarikcilerin payi
# ileri donem akimlar: tedarikciden gelen
bi=b.set_index(['i','j','k','t']).q
for h in range(0,6):
    idx=pd.MultiIndex.from_arrays([pb.i,pb.j,pb.k,pb.t+h]); pb['q%d'%h]=bi.reindex(idx).fillna(0).values
pb.to_parquet(D+'pairs%s.parquet'%TAG)
# ithalatci-urun duzeyi
pb['s2']=pb.s**2
A=pb.groupby(['j','k','t']).agg(E=('E','sum'),Es=('Es','sum'),ncol=('collapse','sum'),hhi=('s2','sum'),nsup=('s','size'),smax=('s','max')).reset_index()
# cokmus tedarikcilerden gelen akim (h)
c=pb[pb.collapse==1].groupby(['j','k','t'])[['bq']+['q%d'%h for h in range(6)]].sum().rename(columns={'bq':'cb',**{'q%d'%h:'c%d'%h for h in range(6)}}).reset_index()
# cokmeyenler icinde hhi ve sayi
nc=pb[pb.collapse==0].copy(); tot=nc.groupby(['j','k','t']).bq.transform('sum'); nc['so2']=(nc.bq/tot)**2
A2=nc.groupby(['j','k','t']).agg(hhi_o=('so2','sum'),nsup_o=('s','size'),sum_o=('bq','sum')).reset_index()
# etkin sayi icin: anlamli tedarikci (pay>=5%)
nc5=nc[nc.s>=0.05].groupby(['j','k','t']).size().rename('n5_o').reset_index()
P=Mf.merge(A,on=['j','k','t'],how='left').merge(c,on=['j','k','t'],how='left').merge(A2,on=['j','k','t'],how='left').merge(nc5,on=['j','k','t'],how='left')
for col in ['E','Es','ncol','cb']+['c%d'%h for h in range(6)]+['n5_o']: P[col]=P[col].fillna(0)
P=P[(P.B>=100)&(P.minpre>0)]
P.to_parquet(D+'panel%s.parquet'%TAG)
ev=pb[pb.collapse==1].drop_duplicates(['i','k','t'])
print('DROP',DROP,'panel',len(P),'ithalatci',P.j.nunique(),'urun',P.k.nunique(),'maruz(E>0)',(P.E>0).sum(),'Es>=0.1',(P.Es>=0.1).sum(),'Es>=0.25',(P.Es>=0.25).sum(),'olay(i,k,t)',len(ev))
'''

# ==============================================================================
# b4_main.py
# ==============================================================================
MODULES['b4_main.py'] = r'''
# Ana sonuclar: maruziyet kutularina gore kayip, ikame, net; ikame orani; heterojenlik; saglamlik
import pandas as pd, numpy as np, json, sys, warnings; warnings.filterwarnings('ignore')
import pyfixest as pf
TAG=sys.argv[1] if len(sys.argv)>1 else ''; NB=int(sys.argv[2]) if len(sys.argv)>2 else 200
P=pd.read_parquet('baci_out/panel%s.parquet'%TAG); pb=pd.read_parquet('baci_out/pairs%s.parquet'%TAG)
# kuresel sikisiklik: cokmus tedarikcilerin dunya payi (agirlikli)
c=pb[pb.collapse==1].copy(); c['w']=c.bq
gw=c.groupby(['j','k','t']).apply(lambda x: np.average(x.wshare.fillna(0),weights=x.w)).rename('cw').reset_index()
P=P.merge(gw,on=['j','k','t'],how='left'); P['cw']=P.cw.fillna(0)
P['jt']=P.j.astype(str)+'_'+P.t.astype(str); P['kt']=P.k+'_'+P.t.astype(str); P['jk']=P.j.astype(str)+'_'+P.k
H=range(6)
for h in H:
    P['yC%d'%h]=((P['c%d'%h]-P.cb)/P.B).clip(upper=1.0); yO=((P['M%d'%h]-P['c%d'%h])-(P.B-P.cb))/P.B; P['yO%d'%h]=yO.clip(upper=yO.quantile(0.99)); P['yT%d'%h]=P['yC%d'%h]+P['yO%d'%h]
    P['uv%d'%h]=(np.log((P['V%d'%h]/P['M%d'%h]).where(P['M%d'%h]>0))-np.log(P.BV/P.B)); P['uv%d'%h]=P['uv%d'%h].clip(P['uv%d'%h].quantile(.01),P['uv%d'%h].quantile(.99))
for h,col in ((4,'Mm4'),(5,'Mm5'),(6,'Mm6')):
    y=(P[col]-P.B)/P.B; P['yPm%d'%h]=y.clip(upper=y.quantile(0.99))
pre=(P.Mm4+P.Mm5+P.Mm6)/3; g=np.log(P.B/pre.where(pre>0)); P['gpre']=g.clip(g.quantile(.01),g.quantile(.99))
edges=[-1,0,0.1,0.25,0.5,0.75,1.01]; labs=['0','b1','b2','b3','b4','b5']; names={'b1':'0-10','b2':'10-25','b3':'25-50','b4':'50-75','b5':'75-100'}
P['Sb']=pd.cut(P.Es,edges,labels=labs).astype(str)
R={'n':len(P),'importers':int(P.j.nunique()),'products':int(P.k.nunique()),'years':[int(P.t.min()),int(P.t.max())]}
ev=pb[pb.collapse==1]
R['events']=int(len(ev.drop_duplicates(['i','k','t']))); R['exporters_with_event']=int(ev.i.nunique())
R['bins']={names[b]:{'n':int((P.Sb==b).sum()),'Es_mean':float(P[P.Sb==b].Es.mean()),'E_mean':float(P[P.Sb==b].E.mean()),'B_median':float(P[P.Sb==b].B.median()),'value_share':float(P[P.Sb==b].BV.sum()/P.BV.sum())} for b in labs[1:]}
R['bins']['none']={'n':int((P.Sb=='0').sum())}
R['exposed_any']=int((P.Es>0).sum()); R['exposed25']=int((P.Es>=0.25).sum())
def est(df,y,rhs='i(Sb, ref="0")',fe='jt + kt',w=None,cl='j + k'):
    d=df.dropna(subset=[y])
    m=pf.feols('%s ~ %s | %s'%(y,rhs,fe),data=d,vcov={'CRV1':cl},weights=w); t=m.tidy()
    return {str(k).replace('C(Sb, contr.treatment(base=0))','').replace('Sb::','').replace('[T.','').replace(']',''):[float(t.loc[k,'Estimate']),float(t.loc[k,'Std. Error'])] for k in t.index},int(m._N)
def binres(df,hs=H,fe='jt + kt',w=None,extra=''):
    out={}
    for h in hs:
        o={}
        for comp in ('T','C','O'):
            e,n=est(df,'y%s%d'%(comp,h),'i(Sb, ref="0")'+extra,fe,w); o[comp]={names[k]:v for k,v in e.items() if k in names}; o['n']=n
        out[h]=o
    return out
ex=pb.groupby(['j','k','t'])[['q%d'%h for h in H]].sum().rename(columns={'q%d'%h:'ex%d'%h for h in H}).reset_index()
P=P.merge(ex,on=['j','k','t'],how='left')
for h in H:
    P['yN%d'%h]=((P['M%d'%h]-P['ex%d'%h].fillna(0))/P.B); P['yN%d'%h]=P['yN%d'%h].clip(upper=P['yN%d'%h].quantile(0.99)); P['yX%d'%h]=P['yO%d'%h]-P['yN%d'%h]
R['source']={h:{comp:{names[k]:v for k,v in est(P,'y%s%d'%(comp,h))[0].items() if k in names} for comp in 'NX'} for h in (0,1,3)}
print('source h1',{b:[round(R['source'][1][c][b][0],3) for c in 'NX'] for b in names.values()},flush=True)
main=binres(P); R['main']=main
for h in H: print(h,{b:[round(main[h][c][b][0],3) for c in 'CTO'] for b in names.values()},flush=True)
if TAG!='':
    json.dump(R,open('baci_out/res%s.json'%TAG,'w'),indent=1); print('alt bitti'); sys.exit(0)
# onceki donem (plasebo)
R['pre']={h:{names[k]:v for k,v in est(P,'yPm%d'%h)[0].items() if k in names} for h in (4,5,6)}; print('pre',{h:{b:round(v[0],3) for b,v in R['pre'][h].items()} for h in R['pre']})
R['uv']={h:{names[k]:v for k,v in est(P,'uv%d'%h)[0].items() if k in names} for h in H}; print('uv',{h:{b:round(v[0],3) for b,v in R['uv'][h].items()} for h in (0,1,3)})
# ikame orani + importer kume bootstrap
def rates(df):
    out={}
    for h in (0,1,2,3,5):
        d=df.dropna(subset=['yT%d'%h])
        mT=pf.feols('yT%d ~ i(Sb, ref="0") | jt + kt'%h,data=d,vcov='iid').coef(); mC=pf.feols('yC%d ~ i(Sb, ref="0") | jt + kt'%h,data=d,vcov='iid').coef()
        for k in mT.index:
            b=[x for x in names if x in str(k)][0]; out[(h,names[b])]=1-mT[k]/mC[k]
    return out
r0=rates(P); rng=np.random.default_rng(1); js=P.j.unique(); grp={j:df for j,df in P.groupby('j')}; bs=[]
for it in range(NB):
    samp=rng.choice(js,len(js)); parts=[]
    for n_,j in enumerate(samp):
        x=grp[j].copy(); x['jt']=str(n_)+'_'+x.t.astype(str); parts.append(x)
    try: bs.append(rates(pd.concat(parts)))
    except Exception as e: pass
R['rate']={'%d|%s'%k:[float(v),float(np.percentile([b[k] for b in bs],2.5)),float(np.percentile([b[k] for b in bs],97.5))] for k,v in r0.items()}; R['rate_nboot']=len(bs)
print('rate',{k:[round(x,2) for x in v] for k,v in R['rate'].items()},flush=True)
json.dump(R,open('baci_out/res_partial.json','w'),indent=1)
# saglamlik
rev=[2002,2007,2012,2017,2022]
S={'excl_HS_revision_years':P[~P.t.isin(rev+[r+1 for r in rev])],'baseline_ge_1000t':P[P.B>=1000],'1998_2009':P[P.t<2010],'2010_2023':P[P.t>=2010],'excl_2020_2021':P[~P.t.isin([2019,2020,2021])]}
NES=['030199','030219','030229','030239','030269','030329','030339','030349','030379','030490','030549','030559','030569','030619','030629','030791','030799','160419','160420','160540','160590']
S['excl_nes_codes']=P[~P.k.isin(NES)]
R['nes_share_obs']=float(P.k.isin(NES).mean())
rob={nm:binres(df,hs=(0,1,3)) for nm,df in S.items()}
rob['importer_product_FE']=binres(P,hs=(0,1,3),fe='jt + kt + jk'); rob['value_weighted']=binres(P,hs=(0,1,3),w='BV'); rob['control_pregrowth']=binres(P.dropna(subset=['gpre']),hs=(0,1,3),extra=' + gpre')
R['rob']=rob
for nm in rob: print(nm,{h:{b:[round(rob[nm][h][c][b][0],2) for c in 'CTO'] for b in ('25-50','50-75','75-100')} for h in (1,)},flush=True)
# heterojenlik: (i) alternatif tedarikci sayisi, (ii) kuresel sikisiklik, (iii) gelir, (iv) urun grubu -- Es>=0.25 olaylari, dogrusal E ile
P['hi']=(P.Es>=0.25).astype(float)
def lin(df,h,z):
    o={}
    for comp in 'TCO':
        e,n=est(df,'y%s%d'%(comp,h),'Es + Es:%s + %s'%(z,z)); o[comp]=e; o['n']=n
    return o
P['alt']=P.n5_o.clip(0,4); P['alt2']=(P.n5_o>=2).astype(float); P['nsup_o_ln']=np.log1p(P.nsup_o.fillna(0)); P['cw10']=(P.cw>=0.10).astype(float)
het={}
for z in ('alt2','nsup_o_ln','cw10','hhi_o'):
    het[z]={h:lin(P,h,z) for h in (0,1,3)}; print(z,{h:{c:{k:[round(x,3) for x in v] for k,v in het[z][h][c].items()} for c in 'O'} for h in (1,)},flush=True)
wb=json.load(open('baci_out/wb_gdppc.json'))[1]; cc=pd.read_csv('baci_out/country_codes_V202601.csv')
g=pd.DataFrame([(x['countryiso3code'],int(x['date']),x['value']) for x in wb if x['value'] is not None],columns=['iso','t','gdp']).merge(cc[['country_code','country_iso3']],left_on='iso',right_on='country_iso3').rename(columns={'country_code':'j'})[['j','t','gdp']]
P=P.merge(g,on=['j','t'],how='left'); P['inc']=P.groupby('t').gdp.transform(lambda s: pd.qcut(s,3,labels=False,duplicates='drop'))
R['gdp_missing_share']=float(P.gdp.isna().mean())
R['het_inc']={['low','middle','high'][int(i)]:binres(df,hs=(0,1,3)) for i,df in P.dropna(subset=['inc']).groupby('inc')}
for i in R['het_inc']: print('inc',i,{b:[round(R['het_inc'][i][1][c][b][0],2) for c in 'CTO'] for b in ('25-50','50-75','75-100')})
def grp(k):
    if len(k)<4: return 'all'
    return {'0301':'live','0302':'fresh','0303':'frozen','0304':'fillets','0305':'cured','0306':'crustaceans','0307':'molluscs','1604':'prepared','1605':'prepared','2301':'meal/oil','1504':'meal/oil'}[k[:4]]
P['grp']=P.k.map(grp); P['E25']=(P.Es>=0.25).astype(float)
hg={}
for gname,df in P.groupby('grp'):
    o={}
    for h in (0,1,3):
        o[h]={comp:est(df,'y%s%d'%(comp,h),'Es',cl='j')[0]['Es'] for comp in 'TCO'}
    hg[gname]=o|{'n':len(df),'n25':int((df.Es>=0.25).sum())}; print(gname,len(df),{h:[round(o[h][c][0],2) for c in 'CTO'] for h in (0,1,3)},flush=True)
R['het_grp']=hg; R['het']=het
# ikamenin kaynagi: mevcut tedarikciler mi yeni tedarikciler mi (Es>=0.25, h=1)
json.dump(R,open(sys.argv[3] if len(sys.argv)>3 else 'baci_out/res%s.json'%TAG,'w'),indent=1)
'''

# ==============================================================================
# run_alt.py
# ==============================================================================
MODULES['run_alt.py'] = r'''
# Alternatif tanimlar: cokus esigi %30 ve %70, dort haneli urun duzeyi, baslangic kisiti olmadan butun cokus yillari
import subprocess, sys
for a in (['0.3','_d30'],['0.7','_d70'],['0.5','_hs4','1','4'],['0.5','_all','0']):
    print('== alternatif',a,flush=True)
    subprocess.run([sys.executable,'b1_events.py']+a,check=True)
    subprocess.run([sys.executable,'-W','ignore','b4_main.py',a[1],'2'],check=True)
'''

# ==============================================================================
# b5_desc.py
# ==============================================================================
MODULES['b5_desc.py'] = r'''
import pandas as pd, numpy as np, json
X=pd.read_parquet('baci_out/X_ikt.parquet'); pb=pd.read_parquet('baci_out/pairs.parquet'); P=pd.read_parquet('baci_out/panel.parquet')
cc=pd.read_csv('baci_out/country_codes_V202601.csv').set_index('country_code').country_name
d=pd.read_parquet('baci_out/seafood.parquet')
D={}
D['raw_rows_all_products']=sum(int(l.split()[2]) for l in open('baci_out/extract.log') if l.startswith('BACI_HS92_Y'))
D['seafood_rows']=int(len(d)); D['value_2023_busd']=float(d[d.t==2023].v.sum()/1e6); D['qty_2023_mt']=float(d[d.t==2023].q.sum()/1e6); D['value_1995_busd']=float(d[d.t==1995].v.sum()/1e6)
ev=pb[pb.collapse==1].groupby(['i','k','t']).agg(nimp=('j','nunique'),n25=('s',lambda s:int((s>=0.25).sum()))).reset_index().merge(X[['i','k','t','B','ratio','wshare','Xr1','Xr3','Xr5']],on=['i','k','t'])
D['events']=len(ev); D['event_exporters']=int(ev.i.nunique()); D['event_products']=int(ev.k.nunique())
D['events_per_year']={int(k):int(v) for k,v in ev.t.value_counts().sort_index().items()}
D['median_ratio']=float(ev.ratio.median()); D['median_Xr1']=float(ev.Xr1.median()); D['median_Xr3']=float(ev.Xr3.median()); D['median_Xr5']=float(ev.Xr5.median())
D['share_below_half_after3']=float((ev.Xr3.dropna()<0.5).mean()); D['share_recovered90_after3']=float((ev.Xr3.dropna()>=0.9).mean()); D['share_recovered90_after5']=float((ev.Xr5.dropna()>=0.9).mean())
D['median_wshare']=float(ev.wshare.median()); D['share_events_wshare_ge10']=float((ev.wshare>=0.10).mean()); D['median_nimp']=float(ev.nimp.median()); D['events_with_dep25']=int((ev.n25>0).sum())
# yillik: cokmus tedarikcilerin taban ihracatinin dunya ticaretine orani
W=d[d.q.notna()].groupby('t').q.sum()
D['collapsed_baseline_share_by_year']={int(t):float(g.B.sum()/W.get(t-1,np.nan)) for t,g in ev.groupby('t')}
# panel ozet
P['Sb']=pd.cut(P.Es,[-1,0,0.1,0.25,0.5,0.75,1.01],labels=['none','0-10','10-25','25-50','50-75','75-100'])
D['panel']={str(b):{'n':int(len(g)),'B_median':float(g.B.median()),'nsup_median':float(g.nsup.median()),'hhi_median':float(g.hhi.median()),'importers':int(g.j.nunique())} for b,g in P.groupby('Sb')}
D['panel_n']=len(P); D['share_obs_dep_ge50']=float((P.Es>=0.5).mean()); D['share_obs_dep_ge25']=float((P.Es>=0.25).mean())
# en buyuk tedarikci payi dagilimi (maruziyetten bagimsiz): tek tedarikciye bagimlilik
D['smax_ge50_share']=float((P.smax>=0.5).mean()); D['smax_ge75_share']=float((P.smax>=0.75).mean()); D['smax_median']=float(P.smax.median())
D['smax_ge50_by_period']={'1998-2009':float((P[P.t<2010].smax>=0.5).mean()),'2010-2023':float((P[P.t>=2010].smax>=0.5).mean())}
# adlandirilmis olaylar
cases=[('Peru','230120',1998,604,'Fish meal'),('Chile','230120',1998,152,'Fish meal'),('Thailand','030613',2013,764,'Frozen shrimp'),('Chile','030322',2010,152,'Frozen Atlantic salmon'),('Chile','030410',2010,152,'Fresh fish fillets')]
out=[]
for nm,k,t,i,lab in cases:
    e=ev[(ev.i==i)&(ev.k==k)&(ev.t==t)]
    if len(e)==0: out.append({'exporter':nm,'product':lab,'year':t,'found':False}); continue
    e=e.iloc[0]; x=pb[(pb.i==i)&(pb.k==k)&(pb.t==t)&(pb.s>=0.25)].merge(P[['j','k','t','B','M0','M1','c0','c1','cb']],on=['j','k','t'],suffixes=('','_p'))
    out.append({'exporter':nm,'product':lab,'hs':k,'year':t,'found':True,'baseline_t':float(e.B),'ratio':float(e.ratio),'wshare':float(e.wshare),'Xr1':float(e.Xr1),'Xr3':float(e.Xr3),'n_importers':int(e.nimp),'n_dep25':int(len(x)),
        'dep25_median_net0':float(((x.M0-x.B_p)/x.B_p).median()) if len(x) else None,'dep25_median_net1':float(((x.M1-x.B_p)/x.B_p).median()) if len(x) else None})
D['cases']=out
top=ev.sort_values('B',ascending=False).head(12); D['top']=[{'exporter':cc.get(r.i,str(r.i)),'hs':r.k,'year':int(r.t),'B':float(r.B),'ratio':float(r.ratio),'wshare':float(r.wshare)} for r in top.itertuples()]
json.dump(D,open('baci_out/desc.json','w'),indent=1); print(json.dumps({k:v for k,v in D.items() if k not in('events_per_year','collapsed_baseline_share_by_year','top')},indent=0)[:3500])
'''

# ==============================================================================
# b6_fig.py
# ==============================================================================
MODULES['b6_fig.py'] = r'''
from figstyle import *
import json, numpy as np
R=json.load(open('baci_out/res.json')); D=json.load(open('baci_out/desc.json'))
for t in ('_d30','_d70','_hs4','_all'):
    try: R['alt'+t]=json.load(open('baci_out/res%s.json'%t))
    except Exception as e: print('yok',t)
R['desc']=D; json.dump(R,open('paper/out/results.json','w'))
bins=['10-25','25-50','50-75','75-100']; col={'10-25':C['lb'],'25-50':C['b'],'50-75':C['o'],'75-100':C['r']}
# Sekil 1
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0))
y=sorted(int(k) for k in D['events_per_year']); ax[0].bar(y,[D['events_per_year'][str(k)] for k in y],color=C['lb'],edgecolor=C['b'],lw=.4); ax[0].set_ylabel('Collapse events (exporter–product)'); ax[0].set_title('(a) Exporter collapses per year')
a2=ax[0].twinx(); a2.plot(y,[100*D['collapsed_baseline_share_by_year'][str(k)] for k in y],color=C['r'],lw=1.2,marker='o',ms=2.5); a2.set_ylabel('Baseline exports of collapsed suppliers\n(% of world seafood trade, tonnes)',color=C['r'],fontsize=7.5); a2.grid(False); a2.spines['right'].set_visible(True)
P=D['panel']; lab=['0-10','10-25','25-50','50-75','75-100']; ax[1].bar(range(5),[P[b]['n'] for b in lab],color=[C['gr']]+[col[b] for b in bins]); ax[1].set_yscale('log'); ax[1].set_xticks(range(5)); ax[1].set_xticklabels(lab); ax[1].set_xlabel('Pre-shock import share from collapsed suppliers (%)'); ax[1].set_ylabel('Importer–product–year observations'); ax[1].set_title('(b) Exposure of import markets')
fig.tight_layout(); fig.savefig('paper/fig/Figure_1.png'); plt.close()
# Sekil 2: olay calismasi
fig,ax=plt.subplots(1,3,figsize=(7.6,2.9),sharey=True); H=[0,1,2,3,4,5]
for a,(comp,tt) in zip(ax,(('C','(a) From collapsed suppliers'),('O','(b) From all other suppliers'),('T','(c) Total imports'))):
    for b in bins:
        e=np.array([R['main'][str(h)][comp][b][0] for h in H]); s=np.array([R['main'][str(h)][comp][b][1] for h in H])
        a.plot(H,100*e,color=col[b],lw=1.3,marker='o',ms=3,label=b+'%'); a.fill_between(H,100*(e-1.96*s),100*(e+1.96*s),color=col[b],alpha=.13,lw=0)
        if comp=='T':
            hp=[-6,-5,-4]; ep=np.array([R['pre'][str(-h)][b][0] for h in hp]); sp=np.array([R['pre'][str(-h)][b][1] for h in hp]); a.plot(hp,100*ep,color=col[b],lw=1,marker='o',ms=2.5,ls='--'); a.fill_between(hp,100*(ep-1.96*sp),100*(ep+1.96*sp),color=col[b],alpha=.10,lw=0)
    a.axhline(0,color='k',lw=.6); a.set_title(tt,fontsize=8.5); a.set_xlabel('Years from collapse')
ax[0].set_ylabel('Change, % of baseline imports'); ax[0].legend(title='Pre-shock share',fontsize=6.5,title_fontsize=6.5,loc='lower left')
fig.tight_layout(); fig.savefig('paper/fig/Figure_2.png'); plt.close()
# Sekil 3: ikame orani
fig,ax=plt.subplots(1,2,figsize=(7.4,3.0)); w=.25
for n_,h in enumerate((0,1,3)):
    v=np.array([R['rate']['%d|%s'%(h,b)] for b in bins]); x=np.arange(4)+(n_-1)*w
    ax[0].bar(x,100*v[:,0],w,color=[C['lb'],C['b'],C['k']][n_],label='Year %d'%h if h else 'Collapse year'); ax[0].errorbar(x,100*v[:,0],yerr=[100*(v[:,0]-v[:,1]),100*(v[:,2]-v[:,0])],fmt='none',ecolor='k',lw=.6,capsize=1.5)
ax[0].axhline(100,color='k',lw=.6,ls=':'); ax[0].axhline(0,color='k',lw=.6); ax[0].set_xticks(range(4)); ax[0].set_xticklabels(bins); ax[0].set_ylim(-40,260); ax[0].set_xlabel('Pre-shock import share from collapsed suppliers (%)'); ax[0].set_ylabel('Replacement rate (% of lost supply)'); ax[0].set_title('(a) Share of lost supply replaced'); ax[0].legend(fontsize=7)
for b in bins:
    e=np.array([R['uv'][str(h)][b][0] for h in H]); s=np.array([R['uv'][str(h)][b][1] for h in H]); ax[1].plot(H,100*(np.exp(e)-1),color=col[b],lw=1.3,marker='o',ms=3,label=b+'%'); ax[1].fill_between(H,100*(np.exp(e-1.96*s)-1),100*(np.exp(e+1.96*s)-1),color=col[b],alpha=.13,lw=0)
ax[1].axhline(0,color='k',lw=.6); ax[1].set_xlabel('Years from collapse'); ax[1].set_ylabel('Change in import unit value (%)'); ax[1].set_title('(b) Import unit values'); ax[1].legend(title='Pre-shock share',fontsize=7,title_fontsize=7)
fig.tight_layout(); fig.savefig('paper/fig/Figure_3.png'); plt.close()
# Sekil 4: heterojenlik
fig,ax=plt.subplots(1,2,figsize=(7.4,3.1))
G=R['het_grp']; ks=sorted(G,key=lambda k: -G[k]['1']['O'][0]/G[k]['1']['C'][0]); yy=np.arange(len(ks))
rr=[-G[k]['1']['O'][0]/G[k]['1']['C'][0] for k in ks]
ax[0].barh(yy,[100*r for r in rr],color=C['b']); ax[0].set_yticks(yy); ax[0].set_yticklabels([k.capitalize() for k in ks]); ax[0].axvline(0,color='k',lw=.6); ax[0].set_xlabel('Replacement per unit of supply lost, year 1 (%)'); ax[0].set_title('(a) By product group')
I=R['het_inc']; x=np.arange(3); bb=['25-50','50-75','75-100']
for n_,(k,lab_) in enumerate((('low','Low income'),('middle','Middle income'),('high','High income'))):
    e=[I[k]['1']['T'][b][0] for b in bb]; s=[I[k]['1']['T'][b][1] for b in bb]; ax[1].bar(x+(n_-1)*w,[100*v for v in e],w,color=[C['r'],C['o'],C['g']][n_],label=lab_); ax[1].errorbar(x+(n_-1)*w,[100*v for v in e],yerr=[196*v for v in s],fmt='none',ecolor='k',lw=.6,capsize=1.5)
ax[1].axhline(0,color='k',lw=.6); ax[1].set_xticks(x); ax[1].set_xticklabels(bb); ax[1].set_xlabel('Pre-shock import share from collapsed suppliers (%)'); ax[1].set_ylabel('Net change in imports, year 1 (% of baseline)'); ax[1].set_title('(b) By importer income'); ax[1].legend(fontsize=7)
fig.tight_layout(); fig.savefig('paper/fig/Figure_4.png'); plt.close(); print('ok')
'''

# ==============================================================================
# figstyle.py
# ==============================================================================
MODULES['figstyle.py'] = r'''
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 9, 'axes.spines.top': False, 'axes.spines.right': False, 'axes.grid': True, 'grid.alpha': .25, 'grid.linewidth': .5,
                     'axes.titlesize': 9.5, 'axes.titleweight': 'bold', 'axes.titlelocation': 'left', 'legend.frameon': False, 'legend.fontsize': 8, 'savefig.dpi': 300, 'savefig.bbox': 'tight', 'figure.dpi': 110})
C = {'b': '#1f5a94', 'r': '#c0392b', 'g': '#2e7d5b', 'o': '#d98200', 'k': '#333333', 'p': '#6a4c93', 'lb': '#8fb8de', 'gr': '#9aa0a6'}
'''

if __name__ == "__main__":
    for _d in DIRS: os.makedirs(_d, exist_ok=True)
    for _m in MODULES:
        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, '-W', 'ignore', _m], check=True, env=dict(os.environ, MPLBACKEND='Agg'))   # her modul ayri surecte calisir / each module runs in its own process
    print('Bitti / done. Sonuclar / results: paper/out/results.json')