# -*- coding: utf-8 -*-
"""Front-loading and selective re-sourcing: Early trade effects of the EU Carbon Border Adjustment Mechanism

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

# # Front-loading and selective re-sourcing: Early trade effects of the EU Carbon Border Adjustment Mechanism
#
# 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: 30–45 dakika.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: Eurostat Comext aylık dış ticaret verisi ve Ember aylık elektrik verisi). 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

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp12", exist_ok=True); os.chdir("/content/rp12")
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 CBAM paper into data/: Eurostat Comext monthly extra-EU imports (DS-045409, reporter EU, flow import)
    by HS4 product and partner, quantity (100 kg) and value (EUR), January 2015 - August 2026, for CBAM goods and comparison goods;
    plus partner carbon intensity (Ember electricity emissions intensity, OWID CO2 per GDP)."""
    import io
    import pandas as pd
    from common import get, DATA
    CBAM = {"Cement": ["2523"], "Electricity": ["2716"], "Fertilisers": ["2808", "2814", "3102", "3105"],
            "Iron and steel": ["7201", "7203", "7205"] + [str(x) for x in range(7206, 7230)] + [str(x) for x in range(7301, 7312)] + ["7318", "7326"],
            "Aluminium": ["7601"] + [str(x) for x in range(7603, 7615)] + ["7616"]}
    CONTROL = {"Ferrous scrap and other steel articles": ["7204"] + [str(x) for x in range(7312, 7318)] + [str(x) for x in range(7319, 7326)],
               "Aluminium scrap and household articles": ["7602", "7615"],
               "Other base metals": [str(x) for x in range(7401, 7420)] + [str(x) for x in range(7501, 7509)] + ["7801", "7802", "7804", "7806"] +
                                    [str(x) for x in range(7901, 7908)] + ["8001", "8002", "8003", "8007"] + [str(x) for x in range(8101, 8114)],
               "Other fertilisers": ["3101", "3103", "3104"],
               "Other inorganic chemicals": [str(x) for x in range(2801, 2854) if str(x) not in ("2808", "2814", "2834")],
               "Lime, plaster, glass and ceramics": ["2520", "2522"] + [str(x) for x in range(7001, 7020)] + [str(x) for x in range(6901, 6915)],
               "Paper": [str(x) for x in range(4801, 4824)]}
    URL = "https://ec.europa.eu/eurostat/api/comext/dissemination/sdmx/2.1/data/DS-045409/M.EU..{p}.1.QUANTITY_IN_100KG+VALUE_IN_EUROS?startPeriod=2015-01&endPeriod=2026-08&format=SDMX-CSV"

    if __name__ == "__main__":
        rows = []
        for grp, d, cb in [(g, CBAM, 1) for g in CBAM] + [(g, CONTROL, 0) for g in CONTROL]:
            for p in d[grp]:
                r = get(URL.format(p=p))
                t = pd.read_csv(io.StringIO(r.text), dtype={"product": str})
                if len(t) == 0: continue
                t = t[["partner", "product", "indicators", "TIME_PERIOD", "OBS_VALUE"]]
                t["group"] = grp; t["cbam"] = cb; rows.append(t); print(grp, p, len(t), flush=True)
        d = pd.concat(rows)
        d = d.pivot_table(index=["group", "cbam", "product", "partner", "TIME_PERIOD"], columns="indicators", values="OBS_VALUE").reset_index()
        d.columns.name = None
        d.rename(columns={"TIME_PERIOD": "month", "QUANTITY_IN_100KG": "qty_100kg", "VALUE_IN_EUROS": "value_eur"}).to_csv(DATA / "comext_imports.csv.gz", index=False, compression="gzip")
        print("done", d.shape)
print(sorted(os.listdir('data')))

# ## 2. CBAM: üçlü fark (PPML), karbon yoğunluğu, olay çalışmaları, birim değerler, Türkiye, sentetik fark-farklar ve plasebo grupları, şekiller

# Paper 12: early trade effects of the EU Carbon Border Adjustment Mechanism (CBAM). Monthly extra-EU imports by HS4 product and partner (Eurostat Comext, January 2019 - August 2026): Poisson pseudo-maximum-likelihood triple differences with high-dimensional fixed effects, continuous treatment by partner carbon intensity, event studies, unit values, synthetic difference-in-differences on product aggregates and partner-specific results for Türkiye.
import json, sys, gc, warnings, pickle, os
import numpy as np, pandas as pd
import pyfixest as pf
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})
RESUME = "--resume" in sys.argv and os.path.exists("out/models.pkl")
R = {}

EU27 = set("AT BE BG CY CZ DE DK EE ES FI FR GR HR HU IE IT LT LU LV MT NL PL PT RO SE SI SK".split())
EXEMPT = {"NO", "IS", "LI", "CH"}                         # EEA EFTA states and Switzerland: outside the scope of CBAM
DROP = EU27 | {"EXT_EU", "INT_EU", "RU", "BY", "UA", "XC", "XL", "QP", "QV", "QW", "QY", "QZ", "AQ", "BV", "HM", "GS", "TF", "UM", "IO", "CX", "CC"}
d = pd.read_csv("data/comext_imports.csv.gz", dtype={"product": str})
d = d[(d.group != "Electricity") & (~d.partner.isin(DROP)) & (d.month >= "2019-01") & (d.month <= "2026-08")]
d["qty"] = d.qty_100kg.fillna(0) / 10.0                    # tonnes
d["value"] = d.value_eur.fillna(0)
d = d.groupby(["group", "cbam", "product", "partner", "month"], as_index=False)[["qty", "value"]].sum()
months = sorted(d.month.unique())
# balanced panel of product-partner pairs with imports in at least 6 months of 2019-2022
pre = d[d.month <= "2022-12"]
act = pre[pre.qty > 0].groupby(["product", "partner"]).size()
pairs = act[act >= 6].index
idx = pd.MultiIndex.from_tuples([(p, j, m) for (p, j) in pairs for m in months], names=["product", "partner", "month"])
info = d.drop_duplicates("product").set_index("product")[["group", "cbam"]]
P = d.set_index(["product", "partner", "month"])[["qty", "value"]].reindex(idx, fill_value=0).reset_index()
P = P.join(info, on="product")
R["panel"] = {"pairs": int(len(pairs)), "products": int(P["product"].nunique()), "cbam_products": int(P[P.cbam == 1]["product"].nunique()),
              "partners": int(P.partner.nunique()), "months": len(months), "obs": int(len(P)), "first": months[0], "last": months[-1]}
del d, pre; gc.collect()

# partner carbon intensity: electricity emissions intensity 2022 (Ember, gCO2/kWh)
E = pd.read_csv("data/ember_monthly.csv.gz")
E = E[(E["Area type"] == "Country or economy") & (E["Electricity source"] == "Total generation") & (E.Date.str[:4] == "2022")]
ci = E.groupby("ISO 3 code").agg(em=("Emissions (MtCO2e)", "sum"), gen=("Generation (TWh)", "sum"))
ci = (ci.em / ci.gen * 1000).rename("ci")               # gCO2/kWh
iso = pd.read_csv("data/iso_codes.csv")                  # alpha-2 to alpha-3
ci = iso.set_index("alpha3").join(ci, how="inner").set_index("alpha2").ci
P["ci"] = P.partner.map(ci)
R["ci_coverage"] = {"partners_with_ci": int(P.drop_duplicates("partner").ci.notna().sum()),
                    "share_of_cbam_imports_with_ci": round(float(P[(P.cbam == 1) & P.ci.notna()].qty.sum() / P[P.cbam == 1].qty.sum()), 3)}

# timing
P["liable"] = (~P.partner.isin(EXEMPT)).astype(int)
P["T1"] = ((P.month >= "2023-10") & (P.month <= "2025-09")).astype(int)      # transitional period (reporting only)
P["T2"] = ((P.month >= "2025-10") & (P.month <= "2025-12")).astype(int)      # last quarter before the definitive period
P["T3"] = (P.month >= "2026-01").astype(int)                                # definitive period
for k in ["T1", "T2", "T3"]: P[f"D_{k}"] = P.cbam * P.liable * P[k]
P["pair"] = (P["product"] + "_" + P.partner).astype("category").cat.codes
P["pm"] = (P["product"] + "_" + P.month).astype("category").cat.codes
P["jm"] = (P.partner + "_" + P.month).astype("category").cat.codes
P["q"] = pd.PeriodIndex(P.month, freq="Q").astype(str)
VC = {"CRV1": "partner"}
FE = "pair + pm + jm"

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

L = P[(P.liable == 1) & P.ci.notna()].copy()
L["ci_s"] = (L.ci - L.drop_duplicates("partner").ci.mean()) / L.drop_duplicates("partner").ci.std()
for k in ["T1", "T2", "T3"]: L[f"C_{k}"] = L.cbam * L.ci_s * L[k]
quarters = sorted(P.q.unique()); ref = "2023Q2"
if RESUME:
    R.update(pickle.load(open('out/models.pkl', 'rb'))); A = R['triple_difference']; B = R['carbon_intensity']; ev = R['event_study']
else:
    # ---------------- A. triple differences (CBAM goods x liable partners x period) ----------------
    A = {}
    m = pf.fepois(f"qty ~ D_T1 + D_T2 + D_T3 | {FE}", data=P, vcov=VC); A["all"] = tidy(m, ["D_T1", "D_T2", "D_T3"]); A["all"]["nobs"] = int(m._N)
    for sec in ["Iron and steel", "Aluminium", "Fertilisers", "Cement"]:
        sub = P[(P.cbam == 0) | (P.group == sec)]
        m = pf.fepois(f"qty ~ D_T1 + D_T2 + D_T3 | {FE}", data=sub, vcov=VC); A[sec] = tidy(m, ["D_T1", "D_T2", "D_T3"])
        print("A", sec, A[sec], flush=True)
    R["triple_difference"] = A

    # ---------------- B. continuous treatment: partner carbon intensity among liable partners ----------------
    B = {}
    m = pf.fepois(f"qty ~ C_T1 + C_T2 + C_T3 | {FE}", data=L, vcov=VC); B["all"] = tidy(m, ["C_T1", "C_T2", "C_T3"]); B["all"]["nobs"] = int(m._N)
    for sec in ["Iron and steel", "Aluminium", "Fertilisers", "Cement"]:
        sub = L[(L.cbam == 0) | (L.group == sec)]
        m = pf.fepois(f"qty ~ C_T1 + C_T2 + C_T3 | {FE}", data=sub, vcov=VC); B[sec] = tidy(m, ["C_T1", "C_T2", "C_T3"])
        print("B", sec, B[sec], flush=True)
    R["carbon_intensity"] = B
    R["ci_sd_gco2_kwh"] = round(float(L.drop_duplicates("partner").ci.std()), 1)
    R["ci_mean_gco2_kwh"] = round(float(L.drop_duplicates("partner").ci.mean()), 1)

    # ---------------- C. event studies by quarter (reference 2023Q2) ----------------
    ev = {}
    for name, data, var in [("liable", P, "liable"), ("carbon_intensity", L, "ci_s")]:
        cols = []
        for qq in quarters:
            if qq == ref: continue
            c = f"E_{qq}"; data[c] = (data.cbam * data[var] * (data.q == qq)).astype("float32"); cols.append(c)
        m = pf.fepois(f"qty ~ {' + '.join(cols)} | {FE}", data=data, vcov=VC)
        t = m.tidy(); ev[name] = {qq: {"b": round(float(t.loc[f"E_{qq}", "Estimate"]), 4), "se": round(float(t.loc[f"E_{qq}", "Std. Error"]), 4)} for qq in quarters if qq != ref}
        data.drop(columns=cols, inplace=True); gc.collect()
        print("event", name, "done", flush=True)
    R["event_study"] = ev

    # ---------------- D. unit values (positive flows) ----------------
    U = L[(L.qty > 0) & (L.value > 0)].copy(); U["luv"] = np.log(U.value / U.qty)
    m = pf.feols(f"luv ~ C_T1 + C_T2 + C_T3 | {FE}", data=U, vcov=VC); R["unit_values_ci"] = tidy(m, ["C_T1", "C_T2", "C_T3"])
    U2 = P[(P.qty > 0) & (P.value > 0)].copy(); U2["luv"] = np.log(U2.value / U2.qty)
    m = pf.feols(f"luv ~ D_T1 + D_T2 + D_T3 | {FE}", data=U2, vcov=VC); R["unit_values_liable"] = tidy(m, ["D_T1", "D_T2", "D_T3"])
    del U, U2; gc.collect()

    # ---------------- E. Türkiye ----------------
    P["TR"] = (P.partner == "TR").astype(int)
    for k in ["T1", "T2", "T3"]: P[f"TR_{k}"] = P.cbam * P.TR * P[k]
    m = pf.fepois(f"qty ~ D_T1 + D_T2 + D_T3 + TR_T1 + TR_T2 + TR_T3 | {FE}", data=P, vcov=VC)
    R["turkiye"] = tidy(m, ["TR_T1", "TR_T2", "TR_T3"])
    tr = P[(P.cbam == 1)].assign(yr=P.month.str[:4]).groupby(["yr", "group"]).apply(lambda g: g[g.partner == "TR"].qty.sum() / g.qty.sum()).unstack()
    R["turkiye_share_of_eu_cbam_imports"] = {g: {y: round(float(v), 4) for y, v in tr[g].items()} for g in tr.columns}
    R["turkiye_ci_gco2_kwh"] = round(float(ci.get("TR", np.nan)), 1)


    pickle.dump(dict(R), open('out/models.pkl', 'wb'))

# ---------------- H. definitive period split at 1 July 2026 (Regulation (EU) 2026/1384 replaces the steel safeguard) ----------------
if os.path.exists("out/models_h.pkl"):
    R["t3_split"] = pickle.load(open("out/models_h.pkl", "rb"))
else:
    P["T3a"] = ((P.month >= "2026-01") & (P.month <= "2026-06")).astype(int); P["T3b"] = (P.month >= "2026-07").astype(int)
    for k in ["T3a", "T3b"]: P[f"D_{k}"] = P.cbam * P.liable * P[k]
    H = {}
    for sec in ["all", "Iron and steel", "Aluminium", "Fertilisers", "Cement"]:
        sub = P if sec == "all" else P[(P.cbam == 0) | (P.group == sec)]
        m = pf.fepois(f"qty ~ D_T1 + D_T2 + D_T3a + D_T3b | {FE}", data=sub, vcov=VC); H[sec] = tidy(m, ["D_T1", "D_T2", "D_T3a", "D_T3b"])
        print("H", sec, H[sec], flush=True); del m; gc.collect()
    P.drop(columns=["T3a", "T3b", "D_T3a", "D_T3b"], inplace=True)
    R["t3_split"] = H; pickle.dump(H, open("out/models_h.pkl", "wb"))

# ---------------- H2. net effect from October 2025 to August 2026 (front-loading plus definitive period) ----------------
if os.path.exists("out/models_h2.pkl"):
    R["net_oct2025_aug2026"] = pickle.load(open("out/models_h2.pkl", "rb"))
else:
    P["T23"] = (P.month >= "2025-10").astype(int); P["D_T23"] = P.cbam * P.liable * P.T23
    H2 = {}
    for sec in ["all", "Iron and steel", "Aluminium", "Fertilisers", "Cement"]:
        sub = P if sec == "all" else P[(P.cbam == 0) | (P.group == sec)]
        m = pf.fepois(f"qty ~ D_T1 + D_T23 | {FE}", data=sub, vcov=VC); H2[sec] = tidy(m, ["D_T1", "D_T23"])
        print("H2", sec, H2[sec], flush=True); del m; gc.collect()
    P.drop(columns=["T23", "D_T23"], inplace=True)
    R["net_oct2025_aug2026"] = H2; pickle.dump(H2, open("out/models_h2.pkl", "wb"))

# ---------------- F. descriptive: volumes and sourcing by carbon intensity ----------------
P["yr"] = P.month.str[:4]
vol = P[P.cbam == 1].groupby(["group", "yr"]).qty.sum().unstack() / 1e6
R["cbam_imports_mt"] = {g: {y: round(float(v), 2) for y, v in vol.loc[g].items()} for g in vol.index}
med = L.drop_duplicates("partner").ci.median()
hi = L[L.cbam == 1].assign(hi=lambda x: x.ci > med, yr=lambda x: x.month.str[:4]).groupby(["group", "yr"]).apply(lambda g: g[g.hi].qty.sum() / g.qty.sum()).unstack()
R["share_high_ci_partners"] = {g: {y: round(float(v), 4) for y, v in hi.loc[g].items()} for g in hi.index}
R["ci_median_gco2_kwh"] = round(float(med), 1)
R["exempt_share"] = {g: {y: round(float(v), 4) for y, v in s.items()} for g, s in
                     P[P.cbam == 1].groupby(["group", "yr"]).apply(lambda g: g[g.liable == 0].qty.sum() / g.qty.sum()).unstack().iterrows()}

# ---------------- G. synthetic difference-in-differences on product aggregates (imports from liable partners) ----------------
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))
    return float(path.mean()), path, om
agg = P[P.liable == 1].groupby(["product", "month"]).qty.sum().unstack()
agg = agg[(agg[[m_ for m_ in months if m_ <= "2022-12"]] > 0).all(1)]
Y = np.log(agg.replace(0, np.nan)).T.interpolate(limit=2).T.dropna()
# deseasonalise with 2019-2022 calendar-month means
base = Y[[m_ for m_ in months if m_ <= "2022-12"]]
seas = pd.DataFrame({mm: base[[c for c in base.columns if c[5:] == mm]].mean(1) for mm in sorted({c[5:] for c in base.columns})})
Yd = Y - seas[[c[5:] for c in Y.columns]].values
pre_c = [m_ for m_ in Y.columns if m_ <= "2023-09"]; post_c = [m_ for m_ in Y.columns if m_ >= "2023-10"]
treated = [p for p in Yd.index if info.loc[p, "cbam"] == 1]; donors = [p for p in Yd.index if info.loc[p, "cbam"] == 0]
SD = {}
for sec in ["all", "Iron and steel", "Aluminium", "Fertilisers", "Cement"]:
    tr_ = treated if sec == "all" else [p for p in treated if info.loc[p, "group"] == sec]
    if not tr_: continue
    att, path, om = sdid(Yd.loc[donors, pre_c].values, Yd.loc[donors, post_c].values, Yd.loc[tr_, pre_c].values, Yd.loc[tr_, post_c].values)
    SD[sec] = {"n_treated": len(tr_), "att_post": round(att, 4),
               "att_transitional": round(float(np.mean([v for c, v in zip(post_c, path) if c <= "2025-09"])), 4),
               "att_q4_2025": round(float(np.mean([v for c, v in zip(post_c, path) if "2025-10" <= c <= "2025-12"])), 4),
               "att_2026": round(float(np.mean([v for c, v in zip(post_c, path) if c >= "2026-01"])), 4),
               "path": [round(float(v), 4) for v in path]}
    if sec == "all": SD_ALL = (post_c, path)
# placebo: control products as pseudo-treated, one at a time (2026 effect)
rng = np.random.default_rng(12); PLB = {}
for sec in SD:
    k_ = SD[sec]["n_treated"]; paths = []; avg = {"transitional": [], "q4_2025": [], "2026": []}
    for r_ in range(200):
        pt = list(rng.choice(donors, size=k_, replace=False)); dn = [x for x in donors if x not in pt]
        _, pth, _ = sdid(Yd.loc[dn, pre_c].values, Yd.loc[dn, post_c].values, Yd.loc[pt, pre_c].values, Yd.loc[pt, post_c].values)
        paths.append(pth)
        avg["transitional"].append(np.mean([v for c, v in zip(post_c, pth) if c <= "2025-09"]))
        avg["q4_2025"].append(np.mean([v for c, v in zip(post_c, pth) if "2025-10" <= c <= "2025-12"]))
        avg["2026"].append(np.mean([v for c, v in zip(post_c, pth) if c >= "2026-01"]))
        avg.setdefault("oct2025_aug2026", []).append(np.mean([v for c, v in zip(post_c, pth) if c >= "2025-10"]))
    paths = np.array(paths)
    PLB[sec] = {"q05_path": np.quantile(paths, 0.05, axis=0), "q95_path": np.quantile(paths, 0.95, axis=0)}
    SD[sec]["att_oct2025_aug2026"] = round(float(np.mean([v for c, v in zip(post_c, SD[sec]["path"]) if c >= "2025-10"])), 4)
    for kk, key in [("transitional", "att_transitional"), ("q4_2025", "att_q4_2025"), ("2026", "att_2026"), ("oct2025_aug2026", "att_oct2025_aug2026")]:
        a = np.array(avg[kk]); SD[sec][f"placebo_p_{kk}"] = round(float((np.abs(a) >= abs(SD[sec][key])).mean()), 3)
        SD[sec][f"placebo_sd_{kk}"] = round(float(a.std()), 4)
    print("placebo", sec, {k: v for k, v in SD[sec].items() if k != "path"}, flush=True)
R["sdid"] = SD; R["sdid_donors"] = len(donors); R["sdid_placebo_draws"] = 200
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
fig, axs = plt.subplots(1, 2, figsize=(10, 3.6))
for a_, (name, ttl) in zip(axs, [("liable", "(a) CBAM goods × liable partners"), ("carbon_intensity", "(b) CBAM goods × partner carbon intensity (per sd)")]):
    e = ev[name]; qs = list(e); b = np.array([e[q]["b"] for q in qs]); se = np.array([e[q]["se"] for q in qs])
    allq = quarters; xb = [allq.index(q) for q in qs]
    a_.errorbar(xb, 100 * b, yerr=196 * se, fmt="o", ms=3, color="k", capsize=0, lw=0.8)
    a_.axhline(0, color="grey", lw=0.6); a_.axvline(allq.index(ref), color="grey", ls=":", lw=0.8)
    for x0, lab in [("2023Q4", "Transitional"), ("2026Q1", "Definitive")]:
        a_.axvline(allq.index(x0) - 0.5, color="#d62728", lw=0.8, ls="--"); a_.text(allq.index(x0) - 0.4, a_.get_ylim()[1] * 0.9 if a_.get_ylim()[1] > 0 else 1, lab, fontsize=7, color="#d62728")
    a_.set_xticks(range(0, len(allq), 4)); a_.set_xticklabels([allq[i][:4] for i in range(0, len(allq), 4)]); a_.set_title(ttl); a_.set_ylabel("Effect on imports (log points × 100)")
fig.tight_layout(); fig.savefig("fig/fig2_event.png"); plt.show()
# Fig 1: EU imports of CBAM goods (index 2019 = 100) and share from exempt partners / high-intensity partners
mq = P[P.cbam == 1].groupby(["group", "q"]).qty.sum().unstack()
fig, axs = plt.subplots(1, 2, figsize=(10, 3.6))
for g, c_ in zip(mq.index, ["#7f7f7f", "#1f77b4", "#2ca02c", "#d62728"]):
    s = mq.loc[g]; s = s[s.index <= "2026Q2"]; axs[0].plot(range(len(s)), 100 * s / s[[q for q in s.index if q.startswith("2019")]].mean(), color=c_, label=g)
axs[0].set_xticks(range(0, len(s), 4)); axs[0].set_xticklabels([s.index[i][:4] for i in range(0, len(s), 4)]); axs[0].axhline(100, color="grey", lw=0.6)
axs[0].set_ylabel("Imports, 2019 = 100"); axs[0].set_title("(a) EU imports of CBAM goods from non-EU partners"); axs[0].legend(frameon=False, fontsize=7)
for g, c_ in zip(hi.index, ["#7f7f7f", "#1f77b4", "#2ca02c", "#d62728"]):
    axs[1].plot(hi.columns, 100 * hi.loc[g], "o-", ms=3, color=c_, label=g)
axs[1].set_ylabel("Share of imports from partners with above-median\nelectricity carbon intensity (%)"); axs[1].set_title("(b) Sourcing from high-intensity partners")
fig.tight_layout(); fig.savefig("fig/fig1_volumes.png"); plt.show()
# Fig 3: coefficients by sector (triple difference and carbon intensity)
fig, axs = plt.subplots(1, 2, figsize=(10, 3.4))
secs = ["Iron and steel", "Aluminium", "Fertilisers", "Cement"]
for a_, (res, ttl, key) in zip(axs, [(A, "(a) Liable vs exempt partners", "D"), (B, "(b) Per sd of partner carbon intensity", "C")]):
    for j, (k, lab, c_) in enumerate([("T1", "Transitional (Oct 2023 – Sep 2025)", "#9ecae1"), ("T2", "Oct – Dec 2025", "#4292c6"), ("T3", "Definitive (Jan – Aug 2026)", "#08306b")]):
        b = [res[s][f"{key}_{k}"]["b"] for s in secs]; se = [res[s][f"{key}_{k}"]["se"] for s in secs]
        a_.bar(np.arange(4) + (j - 1) * 0.26, 100 * np.array(b), 0.26, yerr=196 * np.array(se), color=c_, capsize=2, label=lab)
    a_.axhline(0, color="grey", lw=0.6); a_.set_xticks(range(4)); a_.set_xticklabels(secs, fontsize=8); a_.set_title(ttl); a_.set_ylabel("Log points × 100")
axs[0].legend(frameon=False, fontsize=7)
fig.tight_layout(); fig.savefig("fig/fig3_sectors.png"); plt.show()
# Fig 4: SDID path, all CBAM products
post_c, path = SD_ALL
fig, ax = plt.subplots(figsize=(9, 3.4))
ax.plot(range(len(path)), 100 * np.array(path), color="k", lw=1)
ax.fill_between(range(len(path)), 100 * PLB["all"]["q05_path"], 100 * PLB["all"]["q95_path"], color="#c6dbef", alpha=0.7, lw=0, label="5–95% range of 200 placebo groups of 59 control products")
ax.axhline(0, color="grey", lw=0.6); ax.set_xticks(range(0, len(post_c), 6)); ax.set_xticklabels([post_c[i] for i in range(0, len(post_c), 6)], fontsize=7)
ax.axvline(post_c.index("2026-01") - 0.5, color="#d62728", ls="--", lw=0.8); ax.set_ylabel("Effect on imports (log points × 100)"); ax.legend(frameon=False, fontsize=7)
ax.set_title("Synthetic difference-in-differences: CBAM products, imports from liable partners", fontsize=9)
fig.tight_layout(); fig.savefig("fig/fig4_sdid.png"); plt.show()
print(json.dumps({k: R[k] for k in R if k not in ("event_study",)}, indent=1, default=str)[:9000])

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

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