# -*- coding: utf-8 -*-
"""Weather and daily port activity: Evidence from 434 ports, 2019–2026

GOMD arastirma makalesi 11 - makaledeki tum sayilari, tablolari ve sekilleri ureten kod (Google Colab surumu).
Colab: kodu tek bir hucreye yapistirip calistirin ya da .ipynb dosyasini Dosya > Not defteri yukle ile acin.
Kaynak: GOMD acik veri platformu (https://gomd.pirireis.edu.tr) - Piri Reis Universitesi ve
T.C. Cevre, Sehircilik ve Iklim Degisikligi Bakanligi Cevre Yonetimi Genel Mudurlugu.
"""

# # Weather and daily port activity: Evidence from 434 ports, 2019–2026
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür (replikasyon paketindeki betiklerle aynı kod, hücrelere bölünmüş hâlde).
#
# **Çalıştırma:** Çalışma zamanı → Tümünü çalıştır (*Runtime → Run all*). Tahmini süre: 50–70 dakika ve yaklaşık 8 GB bellek (1,24 milyon gözlemli Poisson modelleri). Kesilirse son kod hücresini `sys.argv.append('--resume')` ekleyerek yeniden çalıştırın.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: IMF PortWatch (AIS tabanlı günlük liman verisi) ve NASA POWER günlük hava verisi; canlı indirme yaklaşık 900 istek ve 30–60 dakika sürer). Makaledeki sayıları birebir yeniden üretmek için replikasyon paketi ZIP'ini yükleyin (`VERI_KAYNAGI = "zip"`); kaynak kurumlar verilerini revize ederse canlı indirmede küçük farklar olabilir. Verilerin çekim tarihi 8 Kasım 2026'dır.
#
# Şekiller hem hücre altında görünür hem `fig/` klasörüne, sayısal sonuçlar `out/` klasörüne yazılır. Son hücre hepsini ZIP olarak indirir.

# ## 0. Kurulum

# Installs the packages used in the paper. If Colab asks to restart the session, restart and run all cells again.
# [Colab Magic] !pip install -q pyfixest statsmodels

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp11", exist_ok=True); os.chdir("/content/rp11")
for d in ("data", "out", "fig"): os.makedirs(d, exist_ok=True)
if VERI_KAYNAGI == "zip":
    import zipfile
    from google.colab import files
    up = files.upload()                      # choose *_replication_package.zip
    with zipfile.ZipFile(next(iter(up))) as z:
        for n in z.namelist():
            if n.startswith("data/") or n == "entities.json": z.extract(n, ".")
    print(sorted(os.listdir("data")))
print("Calisma klasoru:", os.getcwd())

# ## 1. Veri

# --- common.py dosyasi yaziliyor (notebook'taki %%writefile hucresi) ---
open('common.py', 'w', encoding='utf-8').write('"""Shared download helpers (public APIs; no keys needed)."""\nimport io, json, time, zipfile\nfrom pathlib import Path\nimport requests, pandas as pd\nS = requests.Session(); S.headers["User-Agent"] = "research-replication-client"\nDATA = Path("data"); DATA.mkdir(exist_ok=True)\n\ndef get(url, tries=6, **kw):\n    for k in range(tries):\n        try:\n            r = S.get(url, timeout=300, **kw)\n            if r.status_code == 200: return r\n            if r.status_code in (429, 500, 502, 503, 504): time.sleep(5 * (k + 1)); continue\n            r.raise_for_status()\n        except requests.RequestException:\n            if k == tries - 1: raise\n            time.sleep(5 * (k + 1))\n    raise RuntimeError(url)\n\ndef power_daily(lat, lon, params, start="19810101", end="20261031"):\n    url = (f"https://power.larc.nasa.gov/api/temporal/daily/point?parameters={\',\'.join(params)}&community=RE"\n           f"&longitude={lon}&latitude={lat}&start={start}&end={end}&format=JSON")\n    d = get(url).json()["properties"]["parameter"]; time.sleep(0.3)\n    return {p: {k: v for k, v in ser.items() if v != -999} for p, ser in d.items()}\n\ndef power_monthly(lat, lon, params, start=1981, end=2026):\n    url = (f"https://power.larc.nasa.gov/api/temporal/monthly/point?parameters={\',\'.join(params)}&community=RE"\n           f"&longitude={lon}&latitude={lat}&start={start}&end={end}&format=JSON")\n    d = get(url).json()["properties"]["parameter"]; time.sleep(0.3)\n    return {p: {k: v for k, v in ser.items() if v != -999 and k[4:] != "13"} for p, ser in d.items()}\n\ndef jsonstat_frame(js):\n    dims = js["id"]; size = js["size"]\n    cats = [sorted(js["dimension"][d]["category"]["index"].items(), key=lambda kv: kv[1]) for d in dims]\n    rows = []\n    for k, v in js["value"].items():\n        k = int(k); idx = []\n        for s in reversed(size):\n            idx.append(k % s); k //= s\n        idx = list(reversed(idx)); rows.append([cats[i][idx[i]][0] for i in range(len(dims))] + [v])\n    return pd.DataFrame(rows, columns=dims + ["value"])\n\ndef eurostat(dataset, query):\n    url = f"https://ec.europa.eu/eurostat/api/dissemination/statistics/1.0/data/{dataset}?format=JSON&lang=EN&{query}"\n    return jsonstat_frame(get(url).json())\n\ndef edgar(fn):\n    """EDGAR 2026 GHG release (JRC/IEA): downloads the zip and extracts the xlsx into data/."""\n    base = "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/EDGAR/datasets/EDGAR_2026_GHG/"\n    with zipfile.ZipFile(io.BytesIO(get(base + fn).content)) as z:\n        for n in z.namelist():\n            if n.endswith(".xlsx"):\n                (DATA / Path(n).name).write_bytes(z.read(n)); print("EDGAR", Path(n).name)\n\ndef worldbank(ind, years="1985:2025"):\n    d = get(f"https://api.worldbank.org/v2/country/all/indicator/{ind}?format=json&per_page=20000&date={years}").json()\n    rows = [(r["countryiso3code"], r["country"]["value"], int(r["date"]), r["value"]) for r in d[1] if r["value"] is not None and r["countryiso3code"]]\n    df = pd.DataFrame(rows, columns=["iso3", "name", "year", "value"]); df.to_csv(DATA / f"wb_{ind}.csv", index=False); print("WB", ind, len(df))\n\ndef gomd_panel(indicator, entities, out_csv):\n    """Series retrieved through the public GOMD open data API (https://gomd.pirireis.edu.tr) with the gomd Python client."""\n    from gomd import GOMDClient\n    client = GOMDClient("https://gomd.pirireis.edu.tr"); recs = []\n    for code, name, etype in entities:\n        for k in range(4):\n            try:\n                rows = client.get_series(indicator=indicator, entity=code, limit=50000, as_df=False); break\n            except Exception:\n                if k == 3: raise\n                time.sleep(2 ** k)\n        recs += [(code, name, etype, str(r["period"]), float(r["value"])) for r in rows if str(r.get("entity_code") or "") == code and r.get("value") is not None]\n    df = pd.DataFrame(recs, columns=["entity", "name", "etype", "period", "value"]).drop_duplicates(["entity", "period"])\n    df.to_csv(out_csv, index=False); print("GOMD", indicator, len(df)); return df\n')

# Downloads the inputs from the original sources (only when VERI_KAYNAGI == 'api').
if VERI_KAYNAGI == 'api':
    """Downloads the inputs of the port-weather paper into data/:
    IMF PortWatch port metadata and daily port activity (AIS-based port calls and cargo estimates, 1 Jan 2019 - 31 Oct 2026)
    for ports with at least 1,000 vessel calls a year, and NASA POWER daily weather at each port (MERRA-2)."""
    import sys, time
    import pandas as pd
    from common import get, power_daily, DATA
    PW = "https://services9.arcgis.com/weJ1QsnbMYJlCHdG/arcgis/rest/services"
    END = "2026-10-31"

    def arcgis(layer, where, fields, order="ObjectId"):
        rows, off = [], 0
        while True:
            r = get(f"{PW}/{layer}/FeatureServer/0/query", params={"where": where, "outFields": fields, "orderByFields": order,
                    "resultOffset": off, "resultRecordCount": 1000, "f": "json"}).json()
            rows += [f["attributes"] for f in r["features"]]
            if not r.get("exceededTransferLimit"): break
            off += 1000
        return pd.DataFrame(rows)

    if __name__ == "__main__":
        ports = arcgis("PortWatch_ports_database", "1=1", "*")
        ports.to_csv(DATA / "ports.csv", index=False)
        sel = ports[ports.vessel_count_total >= 1000].sort_values("vessel_count_total", ascending=False)
        print("ports:", len(ports), "selected:", len(sel), flush=True)
        F = "date,portid,portcalls,portcalls_container,portcalls_dry_bulk,portcalls_general_cargo,portcalls_roro,portcalls_tanker,import,export"
        out = DATA / "port_daily.csv.gz"
        parts = []
        for i, pid in enumerate(sel.portid):
            parts.append(arcgis("Daily_Ports_Data", f"portid='{pid}' AND date<='{END}'", F, order="date"))
            if i % 50 == 0: print("ports", i, flush=True)
        pd.concat(parts).to_csv(out, index=False, compression="gzip")
        rows = []
        for i, r in enumerate(sel.itertuples()):
            d = power_daily(round(r.lat, 3), round(r.lon, 3), ["WS10M_MAX", "T2M_MAX", "PRECTOTCORR", "WS50M_MAX"], start="20180101", end=END.replace("-", ""))
            for p, ser in d.items():
                rows += [(r.portid, p, k, v) for k, v in ser.items()]
            if i % 50 == 0: print("weather", i, flush=True)
        w = pd.DataFrame(rows, columns=["portid", "param", "date", "value"])
        w.pivot_table(index=["portid", "date"], columns="param", values="value").reset_index().to_csv(DATA / "weather_daily.csv.gz", index=False, compression="gzip")
        print("done", flush=True)
print(sorted(os.listdir('data')))

# ## 2. Hava koşulları ve liman faaliyeti: Poisson modelleri, dinamikler, gemi tipleri, uyum ve ampirik Bayes, kayıp liman çağrıları, şekiller

# Paper 11: weather and the daily operations of global ports. Daily AIS-based port calls and cargo volumes at 434 ports (IMF PortWatch, 1 Jan 2019 - 30 Oct 2026) and daily weather at each port (NASA POWER, MERRA-2): Poisson pseudo-maximum-likelihood with high-dimensional fixed effects, distributed leads and lags (displacement versus loss), vessel types, adaptation (port-specific effects with empirical Bayes shrinkage) and attributable losses.
import json, sys, gc, warnings
import numpy as np, pandas as pd
import pyfixest as pf
from pyfixest.estimation import demean
import matplotlib
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
plt.rcParams.update({"font.family": "DejaVu Serif", "font.size": 9, "axes.titlesize": 10, "legend.fontsize": 8,
                     "savefig.dpi": 300, "axes.spines.top": False, "axes.spines.right": False})
QUICK = "--quick" in sys.argv
R = {}

# ---------------- data ----------------
P = pd.read_csv("data/port_daily.csv.gz", parse_dates=["date"])
W = pd.read_csv("data/weather_daily.csv.gz"); W["date"] = pd.to_datetime(W.date.astype(str), format="%Y%m%d")
meta = pd.read_csv("data/ports.csv")
W = W.sort_values(["portid", "date"]).reset_index(drop=True)
W["yr"] = W.date.dt.year
# data-quality screen for precipitation: a port is excluded from the precipitation analysis if any annual mean from 2021 onwards
# exceeds twice the largest annual mean of 2018-2020 by more than 1 mm/day (drifting reanalysis values; see Section 3)
a = W.groupby(["portid", "yr"]).PRECTOTCORR.mean().unstack()
mx = a[[2018, 2019, 2020]].max(1)
bad_p = a.index[((a[[2021, 2022, 2023, 2024, 2025, 2026]].gt(2 * mx, axis=0)) & (a[[2021, 2022, 2023, 2024, 2025, 2026]].sub(mx, axis=0) > 1)).any(axis=1)]
R["precip_excluded_ports"] = meta.set_index("portid").loc[bad_p, "portname"].tolist()
W.loc[W.portid.isin(bad_p), "PRECTOTCORR"] = np.nan
# wind: port-years whose Jan-Oct mean deviates from the 2018-2022 mean by more than 6 standard deviations are set to missing
# (temperature is not screened: recent record heat is a genuine signal)
W["m"] = W.date.dt.month
flag = {}
for v in ["WS10M_MAX"]:
    am = W[W.m <= 10].groupby(["portid", "yr"])[v].mean().unstack()
    z = am.sub(am[[2018, 2019, 2020, 2021, 2022]].mean(1), axis=0).div(am[[2018, 2019, 2020, 2021, 2022]].std(1), axis=0)
    fl = [(p, y) for p in z.index for y in [2023, 2024, 2025, 2026] if abs(z.loc[p, y]) > 6]
    for p, y in fl: W.loc[(W.portid == p) & (W.yr == y), v] = np.nan
    flag[v] = [(meta.set_index("portid").loc[p, "portname"], int(y)) for p, y in fl]
R["wind_flagged_port_years"] = flag["WS10M_MAX"]

W["sw"] = (W.WS10M_MAX >= 12).astype("float32").where(W.WS10M_MAX.notna())
W["hot"] = (W.T2M_MAX >= 35).astype("float32").where(W.T2M_MAX.notna())
W["rain"] = (W.PRECTOTCORR >= 50).astype("float32").where(W.PRECTOTCORR.notna())
g = W.groupby("portid")
for v in ["sw", "hot", "rain"]:
    for L in range(1, 11): W[f"{v}_l{L}"] = g[v].shift(L).astype("float32")
    for L in (1, 2): W[f"{v}_f{L}"] = g[v].shift(-L).astype("float32")
P = P.drop(columns=[c for c in P.columns if c.startswith("import_") or c.startswith("export_")], errors="ignore")
D = P.merge(W.drop(columns=["yr", "m", "WS50M_MAX"], errors="ignore"), on=["portid", "date"], how="left")
del W, P; gc.collect()
for c in D.columns:
    if D[c].dtype == np.float64 and c not in ("import", "export"): D[c] = D[c].astype("float32")
D["cargo"] = (D["import"] + D["export"]).astype("float64")
D["port"] = D.portid.astype("category").cat.codes
D["d"] = D.date.dt.strftime("%Y%m%d").astype(int)
D["pym"] = D.port.astype(np.int64) * 1000000 + D.date.dt.year * 100 + D.date.dt.month
D["pdow"] = D.port.astype(np.int64) * 10 + D.date.dt.dayofweek
for c in ["port", "d", "pym", "pdow"]: D[c] = D[c].astype("category").cat.codes.astype(np.int64)
R["n_ports"] = int(D.portid.nunique()); R["n_obs"] = int(len(D)); R["period"] = [str(D.date.min().date()), str(D.date.max().date())]
R["mean_daily_calls"] = round(float(D.portcalls.mean()), 2)
R["share_days"] = {"wind_ge_12": round(float(D.sw.mean()), 4), "tmax_ge_35": round(float(D.hot.mean()), 4), "precip_ge_50": round(float(D.rain.mean()), 4)}
FE = "pym + d + pdow"
VC = {"CRV1": "port"}

def fit(y, X, data):
    sub_ = data[[y] + list(X) + ["pym", "d", "pdow", "port"]].dropna()
    m_ = pf.fepois(f"{y} ~ {' + '.join(X)} | {FE}", data=sub_, vcov=VC)
    del sub_; gc.collect()
    return m_

def coef_table(m, names):
    t = m.tidy()
    return {n: {"b": round(float(t.loc[n, "Estimate"]), 4), "se": round(float(t.loc[n, "Std. Error"]), 4)} for n in names}

def lincomb(m, names):
    b = m.coef(); V = m._vcov; ix = [list(b.index).index(n) for n in names]
    return round(float(b.iloc[ix].sum()), 4), round(float(np.sqrt(V[np.ix_(ix, ix)].sum())), 4)

# ---------------- A. contemporaneous non-linear responses ----------------
D["w8"] = ((D.WS10M_MAX >= 8) & (D.WS10M_MAX < 10)).astype("float32"); D["w10"] = ((D.WS10M_MAX >= 10) & (D.WS10M_MAX < 12)).astype("float32")
D["w12"] = ((D.WS10M_MAX >= 12) & (D.WS10M_MAX < 14)).astype("float32"); D["w14"] = (D.WS10M_MAX >= 14).astype("float32")
D["t0"] = (D.T2M_MAX < 0).astype("float32"); D["t30"] = ((D.T2M_MAX >= 30) & (D.T2M_MAX < 35)).astype("float32")
D["t35"] = ((D.T2M_MAX >= 35) & (D.T2M_MAX < 40)).astype("float32"); D["t40"] = (D.T2M_MAX >= 40).astype("float32")
D["p10"] = ((D.PRECTOTCORR >= 10) & (D.PRECTOTCORR < 25)).astype("float32"); D["p25"] = ((D.PRECTOTCORR >= 25) & (D.PRECTOTCORR < 50)).astype("float32")
D["p50"] = (D.PRECTOTCORR >= 50).astype("float32")
BINS = ["w8", "w10", "w12", "w14", "t0", "t30", "t35", "t40", "p10", "p25", "p50"]
RESUME = "--resume" in sys.argv
if RESUME:
    _old = json.load(open("out/results_partial.json"))
if RESUME and "contemporaneous_wind_all_ports" in _old:
    for k_ in ["sample_main", "contemporaneous", "contemporaneous_wind_all_ports"]: R[k_] = _old[k_]
    main = R["contemporaneous"]
else:
    S = D[["WS10M_MAX", "T2M_MAX", "PRECTOTCORR", "portcalls", "cargo", "pym", "d", "pdow", "port", "portid"] + BINS].dropna(subset=["WS10M_MAX", "T2M_MAX", "PRECTOTCORR"])
    R["sample_main"] = {"obs": int(len(S)), "ports": int(S.portid.nunique())}
    main = {}
    for y in ["portcalls", "cargo"]:
        m = fit(y, BINS, S); main[y] = coef_table(m, BINS); main[y]["_nobs"] = int(m._N)
    R["contemporaneous"] = main
    # wind only, all ports
    m = fit("portcalls", ["w8", "w10", "w12", "w14"], D[D.WS10M_MAX.notna()]); R["contemporaneous_wind_all_ports"] = coef_table(m, ["w8", "w10", "w12", "w14"])
    del S, m; gc.collect()
    json.dump(R, open("out/results_partial.json", "w"), indent=1, default=str)

# ---------------- B. dynamics: leads and lags ----------------
LAGNAMES = lambda v: [f"{v}_f2", f"{v}_f1", v] + [f"{v}_l{L}" for L in range(1, 11)]
if RESUME and "dynamics" in _old:
    dyn = _old["dynamics"]; R["dynamics"] = dyn
    for v in ["hot", "rain"]: D.drop(columns=[c for c in LAGNAMES(v) if c != v], inplace=True)
else:
    dyn = {}
    for v in ["sw", "hot", "rain"]:
        for y in ["portcalls", "cargo"]:
            m = fit(y, LAGNAMES(v), D)
            ct = coef_table(m, LAGNAMES(v))
            dyn[f"{v}|{y}"] = {"coefs": ct, "same_day": ct[v], "cum_0_10": lincomb(m, [v] + [f"{v}_l{L}" for L in range(1, 11)]),
                               "cum_m2_10": lincomb(m, LAGNAMES(v)), "nobs": int(m._N)}
            print(v, y, dyn[f"{v}|{y}"]["same_day"], dyn[f"{v}|{y}"]["cum_m2_10"], flush=True)
        if v != "sw": D.drop(columns=[c for c in LAGNAMES(v) if c != v], inplace=True)
        gc.collect()
    R["dynamics"] = dyn
    json.dump(R, open("out/results_partial.json", "w"), indent=1, default=str)

# ---------------- C. vessel types ----------------
TYPES = {"Container": "portcalls_container", "Dry bulk": "portcalls_dry_bulk", "Tanker": "portcalls_tanker", "General cargo": "portcalls_general_cargo", "Ro-ro": "portcalls_roro"}
vt = {}
for k, y in TYPES.items():
    m = fit(y, LAGNAMES("sw"), D)
    vt[k] = {"same_day": coef_table(m, ["sw"])["sw"], "lead_1": coef_table(m, ["sw_f1"])["sw_f1"], "cum_m2_10": lincomb(m, LAGNAMES("sw")),
             "mean_daily_calls": round(float(D[y].mean()), 2)}
    print(k, vt[k], flush=True)
R["vessel_types"] = vt
gc.collect()
json.dump(R, open("out/results_partial.json", "w"), indent=1, default=str)

# ---------------- D. adaptation: interactions with port windiness and size ----------------
base_days = D[D.date < "2023-01-01"]
pc = base_days.groupby("portid").agg(windiness=("sw", "mean"), calls=("portcalls", "mean"), cont=("portcalls_container", "sum"), tot=("portcalls", "sum"))
pc["cont_share"] = pc.cont / pc.tot
pc = pc.join(meta.set_index("portid")[["portname", "country", "continent", "lat", "lon"]])
pc["wind_tercile"] = pd.qcut(pc.windiness.rank(method="first"), 3, labels=["low", "middle", "high"])
pc["size_tercile"] = pd.qcut(pc.calls.rank(method="first"), 3, labels=["small", "medium", "large"])
D = D.merge(pc[["wind_tercile", "size_tercile"]].reset_index(), on="portid", how="left")
het = {}
sub = D[["sw", "portcalls", "pym", "d", "pdow", "port", "wind_tercile", "size_tercile"]].dropna(subset=["sw"]).copy()
for col in ["wind_tercile", "size_tercile"]:
    for lev in pc[col].cat.categories:
        sub[f"sw_{lev}"] = (sub.sw * (sub[col] == lev)).astype("float32")
    names = [f"sw_{lev}" for lev in pc[col].cat.categories]
    m = fit("portcalls", names, sub); het[col] = coef_table(m, names)
    het[col]["_threshold_days_share"] = {lev: round(float(pc.loc[pc[col] == lev, "windiness"].mean()), 4) for lev in pc[col].cat.categories}
R["heterogeneity"] = het
del sub; gc.collect()

# port-specific effects (Frisch-Waugh-Lovell after removing the fixed effects), empirical Bayes shrinkage and meta-regression
sub = D[["portid", "sw", "portcalls", "pym", "d", "pdow"]].dropna(subset=["sw"]).reset_index(drop=True)
Y = np.log1p(sub.portcalls.values).astype(np.float64); X = sub.sw.values.astype(np.float64)
fl = sub[["pym", "d", "pdow"]].values.astype(np.uint64)
res, ok = demean(np.column_stack([Y, X]), fl, np.ones(len(sub)))
sub["ry"], sub["rx"] = res[:, 0], res[:, 1]
pe = sub.groupby("portid").apply(lambda g_: pd.Series({"b": (g_.rx * g_.ry).sum() / (g_.rx ** 2).sum() if (g_.rx ** 2).sum() > 0 else np.nan,
    "se": np.sqrt(((g_.ry - g_.rx * ((g_.rx * g_.ry).sum() / (g_.rx ** 2).sum())) ** 2).mean() / (g_.rx ** 2).sum()) if (g_.rx ** 2).sum() > 0 else np.nan,
    "n_sw": float(g_.sw.sum())}))
pe = pe[(pe.n_sw >= 10)].dropna()
mu = np.average(pe.b, weights=1 / pe.se ** 2); tau2 = max(0.0, float(np.var(pe.b, ddof=1) - np.mean(pe.se ** 2)))
pe["b_eb"] = mu + tau2 / (tau2 + pe.se ** 2) * (pe.b - mu)
pe = pe.join(pc)
import statsmodels.api as sm
Xm = pd.DataFrame({"windiness": pe.windiness * 100, "log_calls": np.log(pe.calls), "abs_lat": pe.lat.abs() / 10, "container_share": pe.cont_share * 10})
ols = sm.WLS(pe.b, sm.add_constant(Xm), weights=1 / (tau2 + pe.se ** 2)).fit(cov_type="HC1")
R["port_specific"] = {"n_ports": int(len(pe)), "precision_weighted_mean": round(float(mu), 4), "tau": round(float(np.sqrt(tau2)), 4),
                      "share_negative_eb": round(float((pe.b_eb < 0).mean()), 3),
                      "eb_quantiles": {k: round(float(v), 4) for k, v in pe.b_eb.quantile([0.1, 0.25, 0.5, 0.75, 0.9]).items()},
                      "meta_regression": {k: {"b": round(float(ols.params[k]), 4), "se": round(float(ols.bse[k]), 4), "p": round(float(ols.pvalues[k]), 4)} for k in ols.params.index},
                      "by_continent": {k: round(float(v), 4) for k, v in pe.groupby("continent").b_eb.mean().items()},
                      "most_sensitive": pe.sort_values("b_eb").head(10)[["portname", "country", "b_eb", "windiness"]].round(4).values.tolist()}
pe.to_csv("out/port_effects.csv")
del sub, res, fl, Y, X; gc.collect()

# ---------------- E. attributable port calls ----------------
# Same-day effects of a strong-wind day by the port's windiness tercile (adaptation) and of a heavy-rain day, scaled to the full
# event window (days -2 to +10) with the ratio of window to same-day effects; applied to each port's mean daily calls.
ratio_w = dyn["sw|portcalls"]["cum_m2_10"][0] / dyn["sw|portcalls"]["same_day"]["b"]
ratio_r = dyn["rain|portcalls"]["cum_m2_10"][0] / dyn["rain|portcalls"]["same_day"]["b"]
bw = {lev: het["wind_tercile"][f"sw_{lev}"]["b"] for lev in ["low", "middle", "high"]}
br = dyn["rain|portcalls"]["same_day"]["b"]
port_w = pc.wind_tercile.astype(str).map(lambda l: 1 - np.exp(bw[l] * ratio_w))
D["year"] = D.date.dt.year
days = D.groupby(["year", "portid"]).agg(sw_days=("sw", "sum"), rain_days=("rain", "sum"), calls=("portcalls", "mean"), total=("portcalls", "sum"))
days = days.join(port_w.rename("ew"), on="portid")
days["lost_wind"] = days.sw_days * days.calls * days.ew
days["lost_rain"] = days.rain_days * days.calls * (1 - np.exp(br * ratio_r))
yr = days.groupby(level=0)[["lost_wind", "lost_rain", "total"]].sum()
R["attributable"] = {"window_to_same_day_ratio_wind": round(float(ratio_w), 3), "window_to_same_day_ratio_rain": round(float(ratio_r), 3),
                     "effect_per_strong_wind_day_by_tercile": {k: round(float(1 - np.exp(v * ratio_w)), 4) for k, v in bw.items()},
                     "effect_per_heavy_rain_day": round(float(1 - np.exp(br * ratio_r)), 4),
                     "by_year": {int(k): {"wind_calls": int(v.lost_wind), "rain_calls": int(v.lost_rain),
                                          "wind_share_pct": round(100 * v.lost_wind / v.total, 3), "rain_share_pct": round(100 * v.lost_rain / v.total, 3)} for k, v in yr.iterrows()}}
ap = days.groupby(level=1)[["lost_wind", "lost_rain", "total"]].sum()
ap["share_pct"] = 100 * (ap.lost_wind + ap.lost_rain.fillna(0)) / ap.total
ap = ap.join(meta.set_index("portid")[["portname", "country", "continent"]])
R["attributable_top_ports"] = ap.sort_values("share_pct", ascending=False).head(15)[["portname", "country", "share_pct"]].round(2).values.tolist()
cont = ap.groupby("continent")[["lost_wind", "lost_rain", "total"]].sum()
R["attributable_by_continent"] = {k: {"wind_pct": round(float(100 * v.lost_wind / v.total), 3), "rain_pct": round(float(100 * v.lost_rain / v.total), 3)} for k, v in cont.iterrows()}
R["attributable_total_pct"] = {"wind": round(float(100 * ap.lost_wind.sum() / ap.total.sum()), 3), "rain": round(float(100 * ap.lost_rain.sum() / ap.total.sum()), 3)}
ap.to_csv("out/attributable_by_port.csv")
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
# Fig 1: non-linear responses (wind, temperature, precipitation) for calls and cargo
fig, axs = plt.subplots(1, 3, figsize=(10, 3.4))
groups = [("(a) Daily maximum wind speed at 10 m", ["w8", "w10", "w12", "w14"], ["8–10", "10–12", "12–14", "≥14"], "m/s (reference < 8)"),
          ("(b) Daily maximum temperature", ["t0", "t30", "t35", "t40"], ["<0", "30–35", "35–40", "≥40"], "°C (reference 0–30)"),
          ("(c) Daily precipitation", ["p10", "p25", "p50"], ["10–25", "25–50", "≥50"], "mm/day (reference < 10)")]
for a_, (ttl, keys, labs, xl) in zip(axs, groups):
    for j, (y, c_) in enumerate([("portcalls", "#1f77b4"), ("cargo", "#d62728")]):
        b = np.array([main[y][k]["b"] for k in keys]); se = np.array([main[y][k]["se"] for k in keys])
        x = np.arange(len(keys)) + (j - 0.5) * 0.18
        a_.errorbar(x, 100 * (np.exp(b) - 1), yerr=[100 * (np.exp(b) - np.exp(b - 1.96 * se)), 100 * (np.exp(b + 1.96 * se) - np.exp(b))], fmt="o", color=c_, capsize=2, ms=4,
                    label="Port calls" if y == "portcalls" else "Cargo volume")
    a_.axhline(0, color="grey", lw=0.6); a_.set_xticks(range(len(keys))); a_.set_xticklabels(labs); a_.set_xlabel(xl); a_.set_title(ttl)
axs[0].set_ylabel("Change relative to reference (%)"); axs[0].legend(frameon=False)
fig.tight_layout(); fig.savefig("fig/fig1_response.png"); plt.show()
# Fig 2: leads and lags of a strong-wind day (calls and cargo), heat and heavy rain (calls)
fig, axs = plt.subplots(1, 3, figsize=(10, 3.4), sharey=True)
xs = np.arange(-2, 11)
for a_, (v, ttl) in zip(axs, [("sw", "(a) Wind ≥ 12 m/s"), ("hot", "(b) Maximum temperature ≥ 35 °C"), ("rain", "(c) Precipitation ≥ 50 mm")]):
    for y, c_ in [("portcalls", "#1f77b4"), ("cargo", "#d62728")]:
        ct = dyn[f"{v}|{y}"]["coefs"]; b = np.array([ct[n]["b"] for n in LAGNAMES(v)]); se = np.array([ct[n]["se"] for n in LAGNAMES(v)])
        a_.plot(xs + (0.12 if y == "cargo" else 0), 100 * b, "o-", color=c_, ms=3, lw=1, label="Port calls" if y == "portcalls" else "Cargo volume")
        a_.fill_between(xs + (0.12 if y == "cargo" else 0), 100 * (b - 1.96 * se), 100 * (b + 1.96 * se), color=c_, alpha=0.15, lw=0)
    a_.axhline(0, color="grey", lw=0.6); a_.axvline(0, color="grey", lw=0.6, ls=":"); a_.set_title(ttl); a_.set_xlabel("Days relative to the weather event")
axs[0].set_ylabel("Effect (log points × 100)"); axs[0].legend(frameon=False)
fig.tight_layout(); fig.savefig("fig/fig2_dynamics.png"); plt.show()
# Fig 3: vessel types
fig, ax = plt.subplots(figsize=(6.4, 3.4))
ks = list(vt); x = np.arange(len(ks))
sd = np.array([vt[k]["same_day"]["b"] for k in ks]); sds = np.array([vt[k]["same_day"]["se"] for k in ks])
cm = np.array([vt[k]["cum_m2_10"][0] for k in ks]); cms = np.array([vt[k]["cum_m2_10"][1] for k in ks])
ax.bar(x - 0.2, 100 * sd, 0.4, yerr=196 * sds, color="#9ecae1", capsize=2, label="Same day")
ax.bar(x + 0.2, 100 * cm, 0.4, yerr=196 * cms, color="#08519c", capsize=2, label="Net, days −2 to +10")
ax.axhline(0, color="grey", lw=0.6); ax.set_xticks(x); ax.set_xticklabels(ks); ax.set_ylabel("Effect of a day with wind ≥ 12 m/s\n(log points × 100)"); ax.legend(frameon=False)
fig.tight_layout(); fig.savefig("fig/fig3_vessel_types.png"); plt.show()
# Fig 4: port-specific effects vs windiness and size
fig, axs = plt.subplots(1, 2, figsize=(9, 3.6), sharey=True)
axs[0].scatter(100 * pe.windiness, 100 * pe.b_eb, s=10, alpha=0.6, color="k"); axs[0].set_xscale("symlog", linthresh=1)
axs[0].set_xlabel("Share of days with wind ≥ 12 m/s, 2019–2022 (%)"); axs[0].set_ylabel("Port-specific effect, empirical Bayes\n(log points × 100)")
axs[1].scatter(pe.calls, 100 * pe.b_eb, s=10, alpha=0.6, color="k"); axs[1].set_xscale("log"); axs[1].set_xlabel("Mean daily port calls, 2019–2022")
for a_ in axs: a_.axhline(0, color="grey", lw=0.6); a_.axhline(100 * mu, color="#d62728", lw=0.8, ls="--")
axs[0].set_title("(a) Exposure"); axs[1].set_title("(b) Port size")
fig.tight_layout(); fig.savefig("fig/fig4_port_effects.png"); plt.show()
print(json.dumps({k: R[k] for k in ["n_ports", "n_obs", "period", "share_days", "precip_excluded_ports", "wind_flagged_port_years", "contemporaneous",
                                    "heterogeneity", "port_specific", "attributable", "attributable_by_continent", "attributable_total_pct", "attributable_top_ports"]}, indent=1, default=str)[:8000])

# ## Sonuçlar ve çıktılar

# Packs figures and numerical results for download.
import shutil
from google.colab import files
shutil.make_archive("/content/rp11_outputs", "zip", ".", ".")
files.download("/content/rp11_outputs.zip")