# -*- coding: utf-8 -*-
"""Which shocks bite? Source-market loss, extreme heat and the resilience of European tourism demand, 2015–2026

GOMD arastirma makalesi 08 - 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.
"""

# # Which shocks bite? Source-market loss, extreme heat and the resilience of European tourism demand, 2015–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: 5–10 dakika.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: Eurostat turizm istatistikleri ve NASA POWER günlük sıcaklık). 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 scipy

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp08", exist_ok=True); os.chdir("/content/rp08")
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')

# --- sites.py dosyasi yaziliyor (notebook'taki %%writefile hucresi) ---
open('sites.py', 'w', encoding='utf-8').write('"""Point locations used to query NASA POWER (approximate coordinates, decimal degrees)."""\n# Paper 7: offshore wind sites in European seas (existing, planned or studied offshore wind zones)\nOFFSHORE = [\n ("Dogger Bank", "GBR", "North Sea", 54.75, 1.90), ("Hornsea", "GBR", "North Sea", 53.90, 1.80), ("East Anglia", "GBR", "North Sea", 52.65, 2.60),\n ("Moray Firth", "GBR", "North Sea", 58.20, -2.80), ("Morecambe Bay", "GBR", "Irish Sea", 53.90, -3.70), ("Celtic Sea", "IRL", "Celtic Sea", 51.10, -6.30),\n ("Borssele", "NLD", "North Sea", 51.70, 3.00), ("IJmuiden Ver", "NLD", "North Sea", 52.85, 3.80), ("German Bight", "DEU", "North Sea", 54.30, 6.50),\n ("Horns Rev", "DNK", "North Sea", 55.50, 7.70), ("Sorlige Nordsjo II", "NOR", "North Sea", 56.80, 4.90), ("Utsira Nord", "NOR", "North Sea", 59.30, 4.50),\n ("Kriegers Flak", "DNK", "Baltic Sea", 55.00, 12.90), ("Baltic Power", "POL", "Baltic Sea", 54.90, 17.90), ("Gulf of Riga", "EST", "Baltic Sea", 58.30, 22.00),\n ("Gulf of Bothnia", "SWE", "Baltic Sea", 62.50, 20.50), ("South Brittany", "FRA", "Atlantic", 47.20, -3.00), ("Galicia", "ESP", "Atlantic", 43.40, -8.80),\n ("Viana do Castelo", "PRT", "Atlantic", 41.70, -9.10), ("Canary Islands", "ESP", "Atlantic", 27.90, -15.30), ("Gulf of Lion", "FRA", "Mediterranean", 43.00, 3.80),\n ("Sicily Channel", "ITA", "Mediterranean", 37.30, 12.20), ("Northern Adriatic", "ITA", "Mediterranean", 44.70, 12.90), ("Saros Gulf", "TUR", "Aegean", 40.50, 26.40),\n ("Gokceada-Bozcaada", "TUR", "Aegean", 39.90, 25.70), ("Cyclades", "GRC", "Aegean", 36.80, 25.40), ("Crete", "GRC", "Mediterranean", 35.60, 24.00),\n ("Kiyikoy", "TUR", "Black Sea", 41.75, 28.40), ("Romanian Black Sea", "ROU", "Black Sea", 44.30, 29.40), ("Cyprus South", "CYP", "Mediterranean", 34.50, 33.50),\n]\n# Paper 8: tourism destinations (two to four points per country, spanning main tourism regions)\nTOURISM = {\n "AT": [(48.21, 16.37), (47.26, 11.39), (47.80, 13.04)], "BE": [(50.85, 4.35), (51.23, 2.92)], "BG": [(43.21, 27.91), (42.70, 23.32), (42.50, 27.47)],\n "HR": [(43.51, 16.44), (42.65, 18.09), (45.81, 15.98)], "CY": [(34.99, 34.00), (34.77, 32.42)], "CZ": [(50.08, 14.44), (49.20, 16.61)],\n "DK": [(55.68, 12.57), (56.16, 10.20)], "EE": [(59.44, 24.75), (58.39, 24.50)], "FI": [(60.17, 24.94), (66.50, 25.73)],\n "FR": [(48.86, 2.35), (43.70, 7.27), (44.84, -0.58)], "DE": [(52.52, 13.40), (48.14, 11.58), (53.55, 9.99)], "EL": [(37.98, 23.73), (36.43, 28.22), (35.34, 25.14)],\n "HU": [(47.50, 19.04), (46.90, 18.05)], "IE": [(53.35, -6.26), (53.27, -9.05)], "IT": [(41.90, 12.50), (45.44, 12.32), (40.85, 14.27)],\n "LV": [(56.95, 24.11), (56.97, 23.77)], "LT": [(54.69, 25.28), (55.71, 21.13)], "LU": [(49.61, 6.13)], "MT": [(35.90, 14.51)],\n "NL": [(52.37, 4.90), (51.50, 3.60)], "PL": [(50.06, 19.94), (54.35, 18.65), (49.30, 19.95)], "PT": [(38.72, -9.14), (37.02, -7.93), (41.15, -8.61)],\n "RO": [(44.43, 26.10), (44.17, 28.64), (45.65, 25.60)], "SK": [(48.15, 17.11), (49.13, 20.22)], "SI": [(46.06, 14.51), (45.51, 13.59)],\n "ES": [(40.42, -3.70), (41.39, 2.17), (36.72, -4.42), (39.57, 2.65)], "SE": [(59.33, 18.07), (57.71, 11.97)], "IS": [(64.15, -21.94)],\n "NO": [(59.91, 10.75), (60.39, 5.32)], "CH": [(47.37, 8.54), (46.20, 6.14), (46.69, 7.86)], "TR": [(36.90, 30.70), (41.01, 28.98), (37.04, 27.43)],\n "ME": [(42.29, 18.84)], "MK": [(41.12, 20.80), (42.00, 21.43)], "AL": [(41.33, 19.82), (41.32, 19.45), (39.87, 20.00)], "RS": [(44.79, 20.45)],\n}\n# Paper 9: major container and Mediterranean ports\nPORTS = [\n ("Shanghai", 31.23, 121.47), ("Singapore", 1.26, 103.82), ("Ningbo-Zhoushan", 29.87, 121.55), ("Shenzhen", 22.54, 114.05), ("Guangzhou Nansha", 22.75, 113.60),\n ("Qingdao", 36.07, 120.38), ("Busan", 35.10, 129.04), ("Tianjin", 38.99, 117.70), ("Hong Kong", 22.30, 114.17), ("Port Klang", 3.00, 101.39),\n ("Rotterdam", 51.95, 4.14), ("Jebel Ali", 25.01, 55.06), ("Antwerp", 51.27, 4.35), ("Xiamen", 24.48, 118.08), ("Tanjung Pelepas", 1.36, 103.55),\n ("Kaohsiung", 22.61, 120.28), ("Los Angeles", 33.74, -118.27), ("Hamburg", 53.54, 9.97), ("New York-New Jersey", 40.67, -74.04), ("Laem Chabang", 13.08, 100.88),\n ("Colombo", 6.95, 79.84), ("Cai Mep", 10.53, 107.03), ("Tanjung Priok", -6.10, 106.88), ("Jeddah", 21.48, 39.17), ("Valencia", 39.45, -0.32),\n ("Piraeus", 37.94, 23.63), ("Algeciras", 36.13, -5.44), ("Nhava Sheva", 18.95, 72.95), ("Santos", -23.98, -46.30), ("Tanger Med", 35.89, -5.50),\n ("Savannah", 32.08, -81.09), ("Felixstowe", 51.95, 1.32), ("Bremerhaven", 53.55, 8.58), ("Le Havre", 49.48, 0.11), ("Barcelona", 41.35, 2.17),\n ("Gioia Tauro", 38.44, 15.90), ("Marsaxlokk", 35.83, 14.54), ("Ambarli", 40.97, 28.68), ("Mersin", 36.79, 34.64), ("Izmir", 38.44, 27.14),\n ("Port Said", 31.26, 32.30), ("Durban", -29.87, 31.03), ("Lagos Apapa", 6.44, 3.37), ("Mombasa", -4.06, 39.67), ("Djibouti", 11.60, 43.13),\n ("Salalah", 16.94, 54.00), ("Houston", 29.73, -95.27), ("Vancouver", 49.29, -123.11), ("Manzanillo", 19.05, -104.31), ("Cartagena", 10.39, -75.52),\n ("Callao", -12.05, -77.15), ("Melbourne", -37.84, 144.92), ("Sydney Botany", -33.97, 151.22), ("Tokyo", 35.62, 139.78), ("Novorossiysk", 44.72, 37.78),\n ("Constanta", 44.17, 28.66), ("Gdansk", 54.40, 18.67), ("Gothenburg", 57.69, 11.86), ("St Petersburg", 59.88, 30.21), ("Balboa", 8.95, -79.57),\n]\n')

# Downloads the inputs from the original sources (only when VERI_KAYNAGI == 'api').
if VERI_KAYNAGI == 'api':
    """Re-downloads the inputs of the tourism paper into data/: Eurostat tour_occ_nim (monthly nights spent by residence) and tour_occ_ninraw
    (annual nights by country of origin), and daily maximum temperature at the 80 destination points (NASA POWER), 1 Jan 1995 - 31 Oct 2026."""
    import pandas as pd
    from common import eurostat, power_daily, DATA
    from sites import TOURISM
    frames = [eurostat("tour_occ_nim", f"unit=NR&nace_r2=I551-I553&c_resid={r}&sinceTimePeriod=2000-01") for r in ("TOTAL", "DOM", "FOR")]
    pd.concat(frames).to_csv(DATA / "eurostat_tour_occ_nim.csv", index=False)
    eurostat("tour_occ_ninraw", "unit=NR&nace_r2=I551-I553&sinceTimePeriod=2015").to_csv(DATA / "eurostat_tour_occ_ninraw.csv", index=False)
    cols = {}
    for geo, pts in TOURISM.items():
        for j, (lat, lon) in enumerate(pts):
            cols[f"{geo}_{j}"] = pd.Series(power_daily(lat, lon, ["T2M_MAX"], start="19950101")["T2M_MAX"]); print(geo, j)
    df = pd.DataFrame(cols); df.index.name = "date"; df[sorted(df.columns)].to_csv(DATA / "tmax_daily_points.csv.gz", compression="gzip")
print(sorted(os.listdir('data')))

# ## 2. Kaynak-pazar şoku, yabani küme önyüklemesi, sentetik fark-farklar, sıcaklık ve şekiller

# Paper 8: which shocks bite? The loss of Russian source markets, extreme heat and the resilience of European tourism demand, 2015-2026. Data: Eurostat monthly nights spent (tour_occ_nim, Jan 2015 - Aug/Sep 2026) by residence (total, domestic, foreign); Eurostat annual nights by country of origin (tour_occ_ninraw, 2015-2025); NASA POWER daily maximum temperature at 80 destination points.
import json, sys, time, warnings
import numpy as np, pandas as pd
import pyfixest as pf
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
plt.rcParams.update({"font.family": "DejaVu Serif", "font.size": 9, "axes.titlesize": 9.5, "legend.fontsize": 7.5,
                     "savefig.dpi": 300, "axes.spines.top": False, "axes.spines.right": False})
QUICK = "--quick" in sys.argv; BW = 99 if QUICK else 1999
T0 = time.time(); R = {}; rng = np.random.default_rng(42)
def log(*a): print(f"[{time.time() - T0:6.0f}s]", *a, flush=True)
AGG = ("EU", "MD", "XK")
NONCLOSERS = {"TR", "RS"}            # destinations that kept air links with Russia after February 2022 (diversion cases)
EXCL_COVID = ("2020-01", "2021-12")
BREAK = {"ME"}                       # Montenegro: break in origin statistics in 2017 (Russian nights fall from 2.9 to 0.6 million in one year)

# =============================== data ===============================
m = pd.read_csv("data/eurostat_tour_occ_nim.csv"); m = m[~m.geo.str.startswith(AGG)]
m["year"] = m.time.str[:4].astype(int); m["month"] = m.time.str[5:].astype(int); m = m[m.year >= 2015]
a = pd.read_csv("data/eurostat_tour_occ_ninraw.csv"); a = a[~a.geo.str.startswith(AGG)]
A = a.pivot_table(index=["geo", "time"], columns="c_resid", values="value")
A["FORX"] = A["WORLD"] - A["DOM"]
A19 = A.xs(2019, level="time")
EXPO = pd.DataFrame({"ru": A19.RU / A19.FORX, "ua": A19.UA / A19.FORX, "cn": A19.CN / A19.FORX, "us": A19.US / A19.FORX,
                     "dom_share": A19.DOM / A19.WORLD}).dropna(subset=["ru"])
# Montenegro and Serbia lack 2019 origin detail: use their latest pre-2020 year with origin data
for g in ("ME", "RS", "NO"):
    if g in A.index.get_level_values(0) and g not in EXPO.index:
        sub = A.xs(g, level="geo"); yrs = [y for y in sub.index if y <= 2019 and np.isfinite(sub.loc[y, "RU"]) and np.isfinite(sub.loc[y, "FORX"]) and sub.loc[y, "FORX"] > 0]
        if yrs:
            y = max(yrs); r = sub.loc[y]; EXPO.loc[g] = [r.RU / r.FORX, r.UA / r.FORX, r.CN / r.FORX, r.US / r.FORX, r.DOM / r.WORLD]
EXPO = EXPO[(EXPO.ru >= 0) & (EXPO.ru < 1)]
EXPO.to_csv("out/exposure_2019.csv")
R["exposure_top"] = EXPO.ru.sort_values(ascending=False).head(12).round(4).to_dict()
ru_idx = A.RU.unstack("time")
R["ru_nights_index_2025_vs_2019"] = (ru_idx[2025] / ru_idx[2019]).dropna().round(3).to_dict()
R["ru_nights_eu27_2019_2025"] = [float(a[(a.geo == "EU27_2020")].shape[0])]
eu = pd.read_csv("data/eurostat_tour_occ_ninraw.csv"); eu = eu[(eu.geo == "EU27_2020") & (eu.c_resid == "RU")].set_index("time").value
R["ru_nights_eu27"] = {int(k): float(v) for k, v in eu.items()}

M = m.pivot_table(index=["geo", "time", "year", "month"], columns="c_resid", values="value").reset_index()
M = M.merge(EXPO, left_on="geo", right_index=True, how="inner")
for k in ("TOTAL", "DOM", "FOR"): M["ln_" + k] = np.log(M[k].where(M[k] > 0))
M["gm"] = M.geo + "_" + M.month.astype(str); M["ym"] = M.time
covid = (M.time >= EXCL_COVID[0]) & (M.time <= EXCL_COVID[1])
# heat exposure: monthly number of days with Tmax >= 32 / 35 C, averaged over each country's points
tw = pd.read_csv("data/tmax_daily_points.csv.gz", index_col="date")
tx = tw.stack().rename("value").reset_index(); tx.columns = ["date", "pt", "value"]
tx["geo"] = tx.pt.str.rsplit("_", n=1).str[0]; tx["point"] = tx.pt.str.rsplit("_", n=1).str[1].astype(int)
tx["time"] = pd.to_datetime(tx.date.astype(str), format="%Y%m%d").dt.strftime("%Y-%m")
H = tx.assign(h32=(tx.value >= 32), h35=(tx.value >= 35)).groupby(["geo", "point", "time"])[["h32", "h35"]].sum().groupby(["geo", "time"]).mean().reset_index()
M = M.merge(H, on=["geo", "time"], how="left")
M.to_csv("out/panel_monthly.csv", index=False)
S_ALL = M[~covid & ~M.geo.isin(NONCLOSERS)].copy()
S = S_ALL[~S_ALL.geo.isin(BREAK)].copy()
R["sample"] = {"countries": int(S.geo.nunique()), "months": [S.time.min(), S.time.max()], "obs_for": int(S.ln_FOR.notna().sum()),
               "geos": sorted(S.geo.unique().tolist())}
log("sample", R["sample"])

# =============================== estimation helpers ===============================
def demean2(X, g1, g2, iters=200, tol=1e-12):
    X = X.astype(float).copy()
    for _ in range(iters):
        X1 = X - pd.DataFrame(X).groupby(g1).transform("mean").values
        X2 = X1 - pd.DataFrame(X1).groupby(g2).transform("mean").values
        if np.max(np.abs(X2 - X)) < tol: X = X2; break
        X = X2
    return X

def crv1(Xd, u, cl):
    XtX = np.linalg.inv(Xd.T @ Xd); G = np.unique(cl); meat = np.zeros((Xd.shape[1],) * 2)
    for g in G:
        s = Xd[cl == g].T @ u[cl == g]; meat += np.outer(s, s)
    n, k = Xd.shape; c = len(G) / (len(G) - 1) * (n - 1) / (n - k)
    return c * XtX @ meat @ XtX

def wild_cluster(yd, Xd, cl, j, B=BW):
    """Wild cluster restricted (WCR) bootstrap-t for H0: beta_j = 0 with Webb six-point weights (MacKinnon & Webb, 2017; Roodman et al., 2019)."""
    b = np.linalg.lstsq(Xd, yd, rcond=None)[0]; u = yd - Xd @ b; V = crv1(Xd, u, cl); t0 = b[j] / np.sqrt(V[j, j])
    Xr = np.delete(Xd, j, axis=1); br = np.linalg.lstsq(Xr, yd, rcond=None)[0]; ur = yd - Xr @ br; fit_r = Xr @ br
    G, inv = np.unique(cl, return_inverse=True); webb = np.array([-np.sqrt(1.5), -1, -np.sqrt(0.5), np.sqrt(0.5), 1, np.sqrt(1.5)]); ts = np.empty(B)
    for k in range(B):
        v = rng.choice(webb, len(G))[inv]; ys = fit_r + ur * v
        bs = np.linalg.lstsq(Xd, ys, rcond=None)[0]; us = ys - Xd @ bs; Vs = crv1(Xd, us, cl); ts[k] = bs[j] / np.sqrt(Vs[j, j])
    return float(t0), float(np.mean(np.abs(ts) >= abs(t0)))

def did(df, y, regs, fe=("gm", "ym"), boot=None):
    d = df.dropna(subset=[y] + regs).copy()
    Z = demean2(np.column_stack([d[y].values] + [d[r].values for r in regs]), d[fe[0]].values, d[fe[1]].values)
    yd, Xd = Z[:, 0], Z[:, 1:]; b = np.linalg.lstsq(Xd, yd, rcond=None)[0]; u = yd - Xd @ b; V = crv1(Xd, u, d.geo.values)
    out = {r: {"b": round(float(b[i]), 4), "se": round(float(np.sqrt(V[i, i])), 4)} for i, r in enumerate(regs)}
    for r in (boot or []):
        i = regs.index(r); t0, p = wild_cluster(yd, Xd, d.geo.values, i); out[r]["t"] = round(t0, 3); out[r]["p_wild"] = round(p, 4)
    out["_n"] = int(len(d)); out["_clusters"] = int(d.geo.nunique())
    return out

# =============================== 1. continuous-exposure event study ===============================
YEARS = [2015, 2016, 2017, 2018, 2022, 2023, 2024, 2025, 2026]
for yv in YEARS: S[f"ru_{yv}"] = S.ru * (S.year == yv); S[f"cn_{yv}"] = S.cn * (S.year == yv); S[f"us_{yv}"] = S.us * (S.year == yv)
EV = {}
for seg in ("FOR", "TOTAL", "DOM"):
    regs = [f"ru_{y}" for y in YEARS]; ctrl = [f"cn_{y}" for y in YEARS] + [f"us_{y}" for y in YEARS]
    EV[seg] = {"base": did(S, "ln_" + seg, regs), "controls": did(S, "ln_" + seg, regs + ctrl)}
    log("event", seg, {k: v["b"] for k, v in EV[seg]["controls"].items() if k.startswith("ru_")})
R["event"] = EV
# pooled post effect with wild cluster bootstrap
S["post"] = (S.year >= 2022).astype(float); S["ru_post"] = S.ru * S.post; S["cn_post"] = S.cn * S.post; S["us_post"] = S.us * S.post
S["ru_trend_pre"] = S.ru * (S.year - 2019) * (S.year <= 2019)
POOL = {}
for seg in ("FOR", "TOTAL", "DOM"):
    POOL[seg] = {"base": did(S, "ln_" + seg, ["ru_post"], boot=["ru_post"]),
                 "controls": did(S, "ln_" + seg, ["ru_post", "cn_post", "us_post"], boot=["ru_post"])}
    log("pooled", seg, POOL[seg]["controls"]["ru_post"])
R["pooled"] = POOL
# pre-trend test: linear pre-2020 trend in exposure
R["pretrend_FOR"] = did(S[S.year <= 2019].assign(ru_t=lambda d: d.ru * (d.year - 2019)), "ln_FOR", ["ru_t"], boot=["ru_t"])
# robustness
ROB = {}
hi = EXPO.index[EXPO.ru >= 0.10]
S["hi_post"] = S.geo.isin(hi).astype(float) * S.post
ROB["binary_high_exposure"] = did(S, "ln_FOR", ["hi_post", "cn_post", "us_post"], boot=["hi_post"])
ROB["drop_2026"] = did(S[S.year <= 2025], "ln_FOR", ["ru_post", "cn_post", "us_post"])
ROB["drop_CY"] = did(S[S.geo != "CY"], "ln_FOR", ["ru_post", "cn_post", "us_post"])
ROB["drop_FI"] = did(S[S.geo != "FI"], "ln_FOR", ["ru_post", "cn_post", "us_post"])
ROB["eu27_only"] = did(S[~S.geo.isin(["AL", "ME", "MK", "CH", "IS", "NO", "LI", "UK"])], "ln_FOR", ["ru_post", "cn_post", "us_post"])
SA = S_ALL.copy(); SA["post"] = (SA.year >= 2022).astype(float)
for v in ("ru", "cn", "us"): SA[v + "_post"] = SA[v] * SA.post
ROB["include_ME"] = did(SA, "ln_FOR", ["ru_post", "cn_post", "us_post"])
ROB["with_ukraine_exposure"] = did(S.assign(ua_post=S.ua * S.post), "ln_FOR", ["ru_post", "ua_post", "cn_post", "us_post"])
ROB["summer_only"] = did(S[S.month.isin([6, 7, 8, 9])], "ln_FOR", ["ru_post", "cn_post", "us_post"])
ROB["winter_only"] = did(S[S.month.isin([11, 12, 1, 2, 3])], "ln_FOR", ["ru_post", "cn_post", "us_post"])
R["robustness"] = ROB
log("robustness", {k: v[list(v)[0]] for k, v in ROB.items()})

# =============================== 2. synthetic difference-in-differences for the most exposed destinations ===============================
def fw(Yw, zeta, x=None, min_decrease=1e-3, max_iter=10000):
    T0_ = Yw.shape[1] - 1; N = Yw.shape[0]; x = np.full(T0_, 1.0 / T0_) if x is None else x
    Yw = Yw - Yw.mean(axis=0, keepdims=True); Am = Yw[:, :T0_]; bv = Yw[:, T0_]; eta = N * zeta ** 2; Ax = Am @ x; prev = None
    for t in range(max_iter):
        hg = (Ax - bv) @ Am + eta * x; i = int(np.argmin(hg)); dx = -x.copy(); dx[i] += 1.0
        if not dx.any(): break
        derr = Am[:, i] - Ax; step = float(np.clip(-(hg @ dx) / (derr @ derr + eta * (dx @ dx)), 0, 1))
        x = x + step * dx; Ax = Ax + step * derr; err = Ax - bv; val = zeta ** 2 * (x @ x) + (err @ err) / N
        if prev is not None and prev - val < min_decrease ** 2: break
        prev = val
    return x
def sparsify(v):
    v = np.where(v <= v.max() / 4, 0, v); return v / v.sum()
def sdid(Ypre_c, Ypost_c, Ypre_t, Ypost_t):
    N1, T1 = Ypost_t.shape; noise = np.std(np.diff(Ypre_c, axis=1), ddof=1)
    zo = (N1 * T1) ** 0.25 * noise; zl = 1e-6 * noise; md = 1e-5 * noise
    Yl = np.column_stack([Ypre_c, Ypost_c.mean(1)]); lam = fw(Yl, zl, min_decrease=md, max_iter=100); lam = fw(Yl, zl, x=sparsify(lam), min_decrease=md)
    Yo = np.column_stack([Ypre_c.T, Ypre_t.mean(0)]); om = fw(Yo, zo, min_decrease=md, max_iter=100); om = fw(Yo, zo, x=sparsify(om), min_decrease=md)
    path = (Ypost_t.mean(0) - Ypre_t.mean(0) @ lam) - (om @ Ypost_c - om @ (Ypre_c @ lam))
    pre_path = (Ypre_t.mean(0) - Ypre_t.mean(0) @ lam) - (om @ Ypre_c - om @ (Ypre_c @ lam))
    return float(path.mean()), path, pre_path, om, lam
# deseasonalised log foreign nights (country-specific 2015-2019 calendar-month means removed)
F = M[~M.geo.isin(NONCLOSERS | BREAK)].pivot(index="geo", columns="time", values="ln_FOR")
months = [c for c in F.columns if (c <= "2019-12") or ("2022-03" <= c <= "2026-07")]
F = F[months]; base = F[[c for c in months if c <= "2019-12"]]
cal = pd.Series([c[5:] for c in base.columns], index=base.columns)
seas = pd.DataFrame({mm: base[cal.index[cal == mm]].mean(1) for mm in sorted(cal.unique())})
Fd = F - seas[[c[5:] for c in F.columns]].values
nmiss = Fd.isna().sum(1); Fd = Fd[nmiss <= 6].interpolate(axis=1, limit=3, limit_area="inside").bfill(axis=1).ffill(axis=1)   # isolated gaps only
R["sdid_interpolated_cells"] = {g: int(n) for g, n in nmiss[(nmiss > 0) & (nmiss <= 6)].items()}
pre = [c for c in Fd.columns if c <= "2019-12"]; post = [c for c in Fd.columns if c >= "2022-03"]
TREATED = [g for g in EXPO.index[EXPO.ru >= 0.10] if g in Fd.index]
DONORS = [g for g in Fd.index if g in EXPO.index and EXPO.ru[g] < 0.03]
SD = {}
for g in TREATED + ["pooled"]:
    tr = TREATED if g == "pooled" else [g]
    att, path, pre_path, om, lam = sdid(Fd.loc[DONORS, pre].values, Fd.loc[DONORS, post].values, Fd.loc[tr, pre].values, Fd.loc[tr, post].values)
    # in-time placebo: pretend the shock happened in March 2018 using 2015-01 .. 2019-12 only
    p_pre = [c for c in pre if c < "2018-03"]; p_post = [c for c in pre if c >= "2018-03"]
    att_p = sdid(Fd.loc[DONORS, p_pre].values, Fd.loc[DONORS, p_post].values, Fd.loc[tr, p_pre].values, Fd.loc[tr, p_post].values)[0]
    SD[g] = {"att": round(att, 4), "pct": round(100 * (np.exp(att) - 1), 1), "placebo_2018": round(att_p, 4),
             "top_donors": {k: round(float(v), 3) for k, v in pd.Series(om, index=DONORS).sort_values(ascending=False).head(4).items()},
             "path_by_year": {y: round(float(np.mean([v for c, v in zip(post, path) if c.startswith(str(y))])), 4) for y in range(2022, 2027)}}
    if g == "pooled": SD_PATH = (post, path, pre, pre_path)
# in-space placebos: each donor in turn as a pseudo-treated unit (pooled-size comparisons use single units)
plac = {}
for g in DONORS:
    dn = [d for d in DONORS if d != g]
    plac[g] = round(sdid(Fd.loc[dn, pre].values, Fd.loc[dn, post].values, Fd.loc[[g], pre].values, Fd.loc[[g], post].values)[0], 4)
R["sdid"] = SD; R["sdid_space_placebos"] = plac
R["sdid_rank_p"] = {g: round(float((np.sum(np.abs(list(plac.values())) >= abs(SD[g]["att"])) + 1) / (len(plac) + 1)), 3) for g in TREATED}
R["sdid_units"] = {"treated": TREATED, "donors": DONORS}
log("sdid", {g: (SD[g]["att"], SD[g]["placebo_2018"]) for g in SD})

# =============================== 3. mechanism: substitution by other origins (annual) ===============================
An = A.reset_index(); An = An[An.geo.isin(EXPO.index) & ~An.geo.isin(NONCLOSERS | BREAK) & ~An.time.isin([2020, 2021])]
An["nonru"] = An.FORX - An.RU.fillna(0); An["ln_nonru"] = np.log(An.nonru.where(An.nonru > 0)); An["ln_forx"] = np.log(An.FORX.where(An.FORX > 0))
An = An.merge(EXPO, left_on="geo", right_index=True); An["post"] = (An.time >= 2022).astype(float); An["ru_post"] = An.ru * An.post
An["cn_post"] = An.cn * An.post; An["us_post"] = An.us * An.post; An["ym"] = An.time.astype(str); An["gm"] = An.geo
MECH = {"non_russian_foreign": did(An, "ln_nonru", ["ru_post", "cn_post", "us_post"]), "all_foreign": did(An, "ln_forx", ["ru_post", "cn_post", "us_post"])}
# accounting decomposition 2019 -> 2025 for the most exposed destinations
dec = {}
for g in TREATED:
    try:
        r19, r25 = A.loc[(g, 2019)], A.loc[(g, 2025)]
        dec[g] = {"foreign_change_pct": round(100 * (r25.FORX / r19.FORX - 1), 1), "russian_loss_pp": round(100 * (r25.RU - r19.RU) / r19.FORX, 1),
                  "other_origins_pp": round(100 * ((r25.FORX - r25.RU) - (r19.FORX - r19.RU)) / r19.FORX, 1)}
    except KeyError: pass
MECH["decomposition_2019_2025"] = dec
R["mechanism"] = MECH
log("mechanism", MECH["non_russian_foreign"]["ru_post"], dec)

# =============================== 4. heat as a contrasting shock ===============================
Hs = M[~covid & ~M.geo.isin(NONCLOSERS | BREAK)].dropna(subset=["h32"]).copy()
HEAT = {}
for seg in ("FOR", "DOM", "TOTAL"):
    HEAT[seg] = did(Hs, "ln_" + seg, ["h32", "h35"], boot=["h32"])
clim = Hs[Hs.month.isin([7, 8])].groupby("geo").h32.mean(); hotg = clim.index[clim >= clim.median()]
Hs["h32_hot"] = Hs.h32 * Hs.geo.isin(hotg); Hs["h32_cool"] = Hs.h32 * ~Hs.geo.isin(hotg)
Hs["sum"] = Hs.month.isin([6, 7, 8, 9])
for seg in ("FOR", "DOM"):
    HEAT[seg + "_by_climate_summer"] = did(Hs[Hs["sum"]], "ln_" + seg, ["h32_hot", "h32_cool"])
Hs = Hs.sort_values(["geo", "time"]); Hs["h32_l12"] = Hs.groupby("geo").h32.shift(12)
HEAT["FOR_lagged_12"] = did(Hs, "ln_FOR", ["h32", "h32_l12"])
R["heat"] = HEAT
R["heat_scale"] = {"mean_h32_julaug_hot_half": round(float(clim[hotg].mean()), 1), "mean_h32_julaug_cool_half": round(float(clim[~clim.index.isin(hotg)].mean()), 1),
                   "sd_h32_julaug_within": round(float(Hs[Hs.month.isin([7, 8])].groupby("geo").h32.std().mean()), 1)}
log("heat", {k: v.get("h32", v.get("h32_hot")) for k, v in HEAT.items()})

# =============================== 5. diversion: Turkiye and Serbia ===============================
DIV = {}
for g in ("TR", "RS"):
    if g in A.index.get_level_values(0):
        s = A.xs(g, level="geo")
        DIV[g] = {"ru_nights": {int(y): float(s.loc[y, "RU"]) for y in s.index if np.isfinite(s.loc[y, "RU"])},
                  "ru_share_of_foreign": {int(y): round(float(s.loc[y, "RU"] / s.loc[y, "FORX"]), 3) for y in s.index if np.isfinite(s.loc[y, "RU"]) and s.loc[y, "FORX"] > 0}}
Ftr = M[M.geo == "TR"].set_index("time").FOR
DIV["TR_foreign_nights_annual"] = {int(y): float(v) for y, v in Ftr.groupby(Ftr.index.str[:4]).sum().items() if Ftr.groupby(Ftr.index.str[:4]).size()[y] == 12} if len(Ftr) else {}
R["diversion"] = DIV
json.dump(R, open("out/results.json", "w"), indent=1, default=str)
log("results written")

# =============================== figures ===============================
fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.3), gridspec_kw={"width_ratios": [1, 1.25]})
e = EXPO.ru.sort_values(ascending=False).head(14)
axs[0].barh(range(len(e)), e.values * 100, color=["#8c564b" if g in NONCLOSERS else "#1f77b4" for g in e.index]); axs[0].set_yticks(range(len(e))); axs[0].set_yticklabels(e.index)
axs[0].invert_yaxis(); axs[0].set_xlabel("Russian share of foreign nights, 2019 (%)"); axs[0].set_title("(a) Exposure", loc="left")
ax = axs[1]
for g, c in [("CY", "#1f77b4"), ("LV", "#2ca02c"), ("EE", "#17becf"), ("FI", "#9467bd"), ("LT", "#bcbd22"), ("TR", "#8c564b"), ("RS", "#d62728")]:
    if g in ru_idx.index:
        s = ru_idx.loc[g]; b0 = s.get(2019) if np.isfinite(s.get(2019, np.nan)) else s.dropna().loc[:2019].iloc[-1]
        s = s / b0 * 100; ax.plot(s.index, s.values, marker="o", ms=2.5, color=c, label=g, lw=1.4 if g in NONCLOSERS else 1)
ax.axvspan(2019.5, 2021.5, color="grey", alpha=0.15); ax.set_yscale("log"); ax.set_ylabel("Russian nights (pre-2020 base = 100, log)"); ax.set_title("(b) Collapse and diversion", loc="left")
ax.legend(frameon=False, ncol=4, fontsize=6.5); fig.tight_layout(); fig.savefig("fig/fig1_exposure.png"); plt.show()

fig, ax = plt.subplots(figsize=(6.4, 3.4))
for k, (seg, c, off) in enumerate([("FOR", "#d62728", -0.15), ("TOTAL", "#1f77b4", 0), ("DOM", "#2ca02c", 0.15)]):
    ev = EV[seg]["controls"]; b = np.array([ev[f"ru_{y}"]["b"] for y in YEARS] + [0]); se = np.array([ev[f"ru_{y}"]["se"] for y in YEARS] + [0]); yy = np.array(YEARS + [2019])
    o = np.argsort(yy); ax.errorbar(yy[o] + off, b[o] / 10, yerr=1.96 * se[o] / 10, fmt="o-", ms=3, lw=0.9, capsize=0, color=c, label={"FOR": "Foreign", "TOTAL": "Total", "DOM": "Domestic"}[seg])
ax.axhline(0, color="k", lw=0.6); ax.axvspan(2019.5, 2021.5, color="grey", alpha=0.15); ax.axvline(2021.9, color="k", ls=":", lw=0.7)
ax.set_ylabel("Effect of +10 pp Russian exposure (log points)"); ax.set_xlabel("Year (2019 = reference; 2020–2021 excluded)"); ax.legend(frameon=False)
fig.tight_layout(); fig.savefig("fig/fig2_event.png"); plt.show()

post_l, path, pre_l, pre_path = SD_PATH
fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.2), gridspec_kw={"width_ratios": [1.5, 1]})
xp = pd.to_datetime(pre_l); xq = pd.to_datetime(post_l)
axs[0].plot(xp, pre_path, color="#7f7f7f", lw=0.8); axs[0].plot(xq, path, color="#d62728", lw=1.0); axs[0].axhline(0, color="k", lw=0.5)
axs[0].set_ylabel("Treated minus synthetic (log points)"); axs[0].set_title("(a) Pooled SDID gap, most exposed destinations", loc="left")
gs = TREATED; vals = [SD[g]["att"] for g in gs]
axs[1].hist(list(plac.values()), bins=12, color="#c7c7c7", label="Placebo donors")
for g, v in zip(gs, vals): axs[1].axvline(v, color="#d62728", lw=1); axs[1].annotate(g, (v, axs[1].get_ylim()[1] * 0.9), fontsize=6.5, rotation=90, ha="right")
axs[1].set_xlabel("ATT, Mar 2022 – Jul 2026 (log points)"); axs[1].set_ylabel("Number of placebo donors"); axs[1].set_title("(b) Country effects vs placebos", loc="left"); axs[1].legend(frameon=False, fontsize=6.5)
fig.tight_layout(); fig.savefig("fig/fig3_sdid.png"); plt.show()

fig, ax = plt.subplots(figsize=(5.8, 3.4))
g25 = {}
for g in EXPO.index:
    if g in NONCLOSERS | BREAK: continue
    try:
        r19, r25 = A.loc[(g, 2019)], A.loc[(g, 2025)]; g25[g] = (EXPO.ru[g], np.log((r25.FORX - r25.RU) / (r19.FORX - r19.RU)), np.log(r25.FORX / r19.FORX))
    except KeyError: pass
G25 = pd.DataFrame(g25, index=["ru", "nonru", "forx"]).T.dropna()
ax.scatter(G25.ru * 100, G25.nonru * 100, s=14, color="#1f77b4", label="Non-Russian foreign nights"); ax.scatter(G25.ru * 100, G25.forx * 100, s=14, marker="x", color="#d62728", label="All foreign nights")
for g, r in G25.iterrows(): ax.annotate(g, (r.ru * 100, r.nonru * 100), fontsize=5.5, xytext=(2, 2), textcoords="offset points")
ax.axhline(0, color="k", lw=0.5); ax.set_xlabel("Russian share of foreign nights, 2019 (%)"); ax.set_ylabel("Change 2019–2025 (log points × 100)"); ax.legend(frameon=False, fontsize=7)
fig.tight_layout(); fig.savefig("fig/fig4_substitution.png"); plt.show()
R["growth_2019_2025"] = {g: {k: round(float(v), 4) for k, v in r.items()} for g, r in G25.iterrows()}

fig, ax = plt.subplots(figsize=(6.0, 3.0))
rows = [("Russian exposure × post-2022 (+10 pp)", "FOR", POOL["FOR"]["controls"]["ru_post"], 0.1), ("", "TOTAL", POOL["TOTAL"]["controls"]["ru_post"], 0.1), ("", "DOM", POOL["DOM"]["controls"]["ru_post"], 0.1),
        ("Ten more days ≥ 32 °C in a month", "FOR", HEAT["FOR"]["h32"], 10), ("", "TOTAL", HEAT["TOTAL"]["h32"], 10), ("", "DOM", HEAT["DOM"]["h32"], 10)]
cols = {"FOR": "#d62728", "TOTAL": "#1f77b4", "DOM": "#2ca02c"}
for k, (lab, seg, r, sc) in enumerate(rows):
    ax.errorbar(r["b"] * sc * 100, k, xerr=1.96 * r["se"] * sc * 100, fmt="o", color=cols[seg], ms=4, capsize=0)
ax.set_yticks(range(len(rows))); ax.set_yticklabels([f"{l}  {s.title() if s != 'FOR' else 'Foreign'}" if l else ({"TOTAL": "Total", "DOM": "Domestic"}[s]) for l, s, _, _ in rows], fontsize=7)
ax.invert_yaxis(); ax.axvline(0, color="k", lw=0.6); ax.set_xlabel("Effect on nights spent (log points × 100), 95% CI")
fig.tight_layout(); fig.savefig("fig/fig5_shocks.png"); plt.show()
json.dump(R, open("out/results.json", "w"), indent=1, default=str)
log("done")
print(json.dumps({k: R[k] for k in ["sample", "exposure_top", "pooled", "pretrend_FOR", "robustness", "sdid", "sdid_rank_p", "mechanism", "heat", "diversion"]}, indent=1, default=str))

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

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