# -*- coding: utf-8 -*-
"""Do carbon taxes bend national emission trajectories? Pre-adoption trends and evidence from ten staggered-adoption estimators, 1980–2025

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

# # Do carbon taxes bend national emission trajectories? Pre-adoption trends and evidence from ten staggered-adoption estimators, 1980–2025
#
# 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 (on tahmin edicinin 199 tekrarlı ülke önyüklemesi).
#
# *Hızlı deneme:* ilk kod hücresinden sonra `import sys; sys.argv.append('--quick')` çalıştırılırsa önyükleme tekrar sayısı azaltılır (sonuçlar makaledekinden farklı olur).
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: EDGAR 2026, Dünya Bankası, Dünya Karbon Fiyatlama Veritabanı (GOMD API üzerinden)). 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 openpyxl
# [Colab Magic] !pip install -q --upgrade gomd

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

import json
if not os.path.exists('entities.json'):
    json.dump([["TUR", "Turkiye", "country"], ["AFG", "Afghanistan", "country"], ["ALB", "Albania", "country"], ["DZA", "Algeria", "country"], ["DEU", "Germany", "country"], ["AND", "Andorra", "country"], ["AGO", "Angola", "country"], ["ATG", "Antigua and Barbuda", "country"], ["ARG", "Argentina", "country"], ["ARM", "Armenia", "country"], ["AUS", "Australia", "country"], ["AUT", "Austria", "country"], ["AZE", "Azerbaijan", "country"], ["BHS", "Bahamas, The", "country"], ["BHR", "Bahrain", "country"], ["BGD", "Bangladesh", "country"], ["BRB", "Barbados", "country"], ["BLR", "Belarus", "country"], ["BEL", "Belgium", "country"], ["BLZ", "Belize", "country"], ["BEN", "Benin", "country"], ["BTN", "Bhutan", "country"], ["BOL", "Bolivia", "country"], ["BIH", "Bosnia and Herzegovina", "country"], ["BWA", "Botswana", "country"], ["BRA", "Brazil", "country"], ["BRN", "Brunei Darussalam", "country"], ["BGR", "Bulgaria", "country"], ["BFA", "Burkina Faso", "country"], ["BDI", "Burundi", "country"], ["CPV", "Cabo Verde", "country"], ["KHM", "Cambodia", "country"], ["CMR", "Cameroon", "country"], ["CAN", "Canada", "country"], ["CAF", "Central African Republic", "country"], ["TCD", "Chad", "country"], ["CHL", "Chile", "country"], ["TWN", "Chinese Taipei", "country"], ["CHN", "China", "country"], ["COL", "Colombia", "country"], ["COM", "Comoros", "country"], ["COG", "Congo, Rep.", "country"], ["CRI", "Costa Rica", "country"], ["CIV", "Cote d'Ivoire", "country"], ["HRV", "Croatia", "country"], ["CUB", "Cuba", "country"], ["CYP", "Greek Cypriot Administration of Southern Cyprus (GCA)", "country"], ["CZE", "Czechia", "country"], ["COD", "Democratic Republic of Congo", "country"], ["DNK", "Denmark", "country"], ["DJI", "Djibouti", "country"], ["DMA", "Dominica", "country"], ["DOM", "Dominican Republic", "country"], ["ECU", "Ecuador", "country"], ["EGY", "Egypt", "country"], ["SLV", "El Salvador", "country"], ["GNQ", "Equatorial Guinea", "country"], ["ERI", "Eritrea", "country"], ["EST", "Estonia", "country"], ["SWZ", "Eswatini", "country"], ["ETH", "Ethiopia", "country"], ["FJI", "Fiji", "country"], ["FIN", "Finland", "country"], ["FRA", "France", "country"], ["GAB", "Gabon", "country"], ["GMB", "Gambia, The", "country"], ["GEO", "Georgia", "country"], ["GHA", "Ghana", "country"], ["GRD", "Grenada", "country"], ["GTM", "Guatemala", "country"], ["GIN", "Guinea", "country"], ["GNB", "Guinea-Bissau", "country"], ["GUY", "Guyana", "country"], ["HTI", "Haiti", "country"], ["HND", "Honduras", "country"], ["HKG", "Hong Kong SAR, China", "country"], ["HUN", "Hungary", "country"], ["ISL", "Iceland", "country"], ["IND", "India", "country"], ["IDN", "Indonesia", "country"], ["IRN", "Iran, Islamic Rep.", "country"], ["IRQ", "Iraq", "country"], ["IRL", "Ireland", "country"], ["ESP", "Spain", "country"], ["ISR", "Israel", "country"], ["ITA", "Italy", "country"], ["JAM", "Jamaica", "country"], ["JPN", "Japan", "country"], ["JOR", "Jordan", "country"], ["KAZ", "Kazakhstan", "country"], ["KEN", "Kenya", "country"], ["KIR", "Kiribati", "country"], ["KOR", "Korea Rep", "country"], ["PRK", "Korea, Dem. People's Rep.", "country"], ["KWT", "Kuwait", "country"], ["KGZ", "Kyrgyzstan", "country"], ["LAO", "Laos", "country"], ["LVA", "Latvia", "country"], ["LBN", "Lebanon", "country"], ["LSO", "Lesotho", "country"], ["LBR", "Liberia", "country"], ["LBY", "Libya", "country"], ["LIE", "Liechtenstein", "country"], ["LTU", "Lithuania", "country"], ["LUX", "Luxembourg", "country"], ["MAC", "Macao SAR, China", "country"], ["MDG", "Madagascar", "country"], ["MWI", "Malawi", "country"], ["MYS", "Malaysia", "country"], ["MDV", "Maldives", "country"], ["MLI", "Mali", "country"], ["MLT", "Malta", "country"], ["MHL", "Marshall Islands", "country"], ["MRT", "Mauritania", "country"], ["MUS", "Mauritius", "country"], ["MEX", "Mexico", "country"], ["FSM", "Micronesia, Fed. Sts.", "country"], ["MDA", "Moldova", "country"], ["MCO", "Monaco", "country"], ["MNG", "Mongolia", "country"], ["MNE", "Montenegro", "country"], ["MAR", "Morocco", "country"], ["MOZ", "Mozambique", "country"], ["MMR", "Myanmar", "country"], ["NAM", "Namibia", "country"], ["NRU", "Nauru", "country"], ["NPL", "Nepal", "country"], ["NLD", "Netherlands", "country"], ["NZL", "New Zealand", "country"], ["NIC", "Nicaragua", "country"], ["NER", "Niger", "country"], ["NGA", "Nigeria", "country"], ["MKD", "North Macedonia", "country"], ["NOR", "Norway", "country"], ["OMN", "Oman", "country"], ["PAK", "Pakistan", "country"], ["PLW", "Palau", "country"], ["PAN", "Panama", "country"], ["PNG", "Papua New Guinea", "country"], ["PRY", "Paraguay", "country"], ["PER", "Peru", "country"], ["PHL", "Philippines", "country"], ["POL", "Poland", "country"], ["PRT", "Portugal", "country"], ["PRI", "Puerto Rico", "country"], ["QAT", "Qatar", "country"], ["ROU", "Romania", "country"], ["RUS", "Russian Federation", "country"], ["RWA", "Rwanda", "country"], ["WSM", "Samoa", "country"], ["SMR", "San Marino", "country"], ["STP", "Sao Tome and Principe", "country"], ["SAU", "Saudi Arabia", "country"], ["SEN", "Senegal", "country"], ["SRB", "Serbia", "country"], ["SYC", "Seychelles", "country"], ["SLE", "Sierra Leone", "country"], ["SGP", "Singapore", "country"], ["SVK", "Slovakia", "country"], ["SVN", "Slovenia", "country"], ["SLB", "Solomon Islands", "country"], ["SOM", "Somalia", "country"], ["ZAF", "South Africa", "country"], ["SSD", "South Sudan", "country"], ["LKA", "Sri Lanka", "country"], ["KNA", "St. Kitts and Nevis", "country"], ["LCA", "St. Lucia", "country"], ["VCT", "St. Vincent and the Grenadines", "country"], ["SDN", "Sudan", "country"], ["SUR", "Suriname", "country"], ["SWE", "Sweden", "country"], ["CHE", "Switzerland", "country"], ["SYR", "Syrian Arab Republic", "country"], ["TJK", "Tajikistan", "country"], ["TZA", "Tanzania", "country"], ["THA", "Thailand", "country"], ["TLS", "Timor-Leste", "country"], ["TGO", "Togo", "country"], ["TON", "Tonga", "country"], ["TTO", "Trinidad and Tobago", "country"], ["TUN", "Tunisia", "country"], ["TKM", "Turkmenistan", "country"], ["TUV", "Tuvalu", "country"], ["UGA", "Uganda", "country"], ["GBR", "UK", "country"], ["UKR", "Ukraine", "country"], ["ARE", "United Arab Emirates", "country"], ["URY", "Uruguay", "country"], ["USA", "USA", "country"], ["UZB", "Uzbekistan", "country"], ["VUT", "Vanuatu", "country"], ["VAT", "Vatican", "country"], ["VEN", "Venezuela, RB", "country"], ["VNM", "Vietnam", "country"], ["PSE", "West Bank and Gaza", "country"], ["ESH", "Western Sahara", "country"], ["YEM", "Yemen", "country"], ["GRC", "Greece", "country"], ["ZMB", "Zambia", "country"], ["ZWE", "Zimbabwe", "country"]], open('entities.json', 'w', encoding='utf-8'), ensure_ascii=False)

# Downloads the inputs from the original sources (only when VERI_KAYNAGI == 'api').
if VERI_KAYNAGI == 'api':
    """Re-downloads the inputs of the carbon tax paper into data/: EDGAR 2026 release (CO2 and CH4, 1970-2025), World Bank WDI
    (population, GDP per capita PPP) and the emissions-weighted carbon tax rate of the World Carbon Pricing Database (via the GOMD API)."""
    import json
    from common import edgar, worldbank, gomd_panel
    edgar("IEA_EDGAR_CO2_1970_2025.zip"); edgar("EDGAR_CH4_1970_2025.zip")
    for ind in ("SP.POP.TOTL", "NY.GDP.PCAP.PP.KD"): worldbank(ind)
    gomd_panel("carbon_tax_price_weighted", json.load(open("entities.json", encoding="utf-8")), "data/ctax.csv")
print(sorted(os.listdir('data')))

# ## 2. Panel, on tahmin edici, plasebo testleri, heterojenlik ve şekiller

# Paper 6: carbon tax adoption and national CO2 emissions, 1980-2025. Ten staggered-adoption estimators are compared under one falsification test (adoption moved three years earlier, pre-adoption data only). Data: EDGAR 2026 release (CO2 and CH4, 1970-2025), World Carbon Pricing Database emissions-weighted carbon tax price (Dolphin & Xiahou, 2022; via GOMD), World Bank WDI (population, GDP per capita PPP).
import json, sys, time, warnings
import numpy as np, pandas as pd
import matplotlib
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
B_MAIN = 5 if QUICK else 199; B_HET = 5 if QUICK else 99
rng = np.random.default_rng(2026); R = {}; T0 = time.time()
Y0, Y1 = 1980, 2025; EMIN, EHI, SHIFT = -10, 10, 3
def log(*a): print(f"[{time.time() - T0:7.0f}s]", *a, flush=True)

# =============================== data ===============================
def edgar(path, sheet="TOTALS BY COUNTRY"):
    return pd.read_excel(path, sheet_name=sheet, header=9)
co2_tot = edgar("data/IEA_EDGAR_CO2_1970_2025.xlsx")
co2_sec = edgar("data/IEA_EDGAR_CO2_1970_2025.xlsx", "IPCC 2006")
ch4_sec = edgar("data/EDGAR_CH4_1970_2025.xlsx", "IPCC 2006")
YC = [f"Y_{y}" for y in range(1970, 2026)]
def long(df, name):
    d = df.groupby("Country_code_A3")[YC].sum(min_count=1).stack().reset_index()
    d.columns = ["iso", "yc", name]; d["year"] = d.yc.str[2:].astype(int); d[name] = d[name] / 1000.0   # Gg -> Mt
    return d.drop(columns="yc")
code = co2_sec.ipcc_code_2006_for_standard_report.astype(str)
P = long(co2_tot, "co2")
SECTORS = {"co2_power": code.str.startswith("1.A.1"), "co2_industry": code.str.startswith("1.A.2") | code.str.startswith("2."),
           "co2_transport": code.str.startswith("1.A.3"), "co2_buildings": code.str.startswith("1.A.4")}
for nm, mask in SECTORS.items():
    P = P.merge(long(co2_sec[mask], nm), on=["iso", "year"], how="left")
ccode = ch4_sec.ipcc_code_2006_for_standard_report.astype(str)
P = P.merge(long(ch4_sec[ccode.str.startswith("3.")], "ch4_agri"), on=["iso", "year"], how="left")
names = co2_tot.set_index("Country_code_A3").Name.to_dict()
wb = {k: pd.read_csv(f"data/wb_{k}.csv") for k in ("SP.POP.TOTL", "NY.GDP.PCAP.PP.KD")}
P = P.merge(wb["SP.POP.TOTL"].rename(columns={"iso3": "iso", "value": "pop"})[["iso", "year", "pop"]], on=["iso", "year"], how="left")
P = P.merge(wb["NY.GDP.PCAP.PP.KD"].rename(columns={"iso3": "iso", "value": "gdppc"})[["iso", "year", "gdppc"]], on=["iso", "year"], how="left")
tx = pd.read_csv("data/ctax.csv", dtype={"period": str}); tx = tx[tx.etype == "country"]; tx["year"] = tx.period.astype(int)
TX = tx.pivot(index="entity", columns="year", values="value")
first_any = TX.gt(0).idxmax(axis=1)[TX.gt(0).any(axis=1)]
price2025 = TX[2025]
P = P[(P.year >= Y0) & (P.year <= Y1)]
pop15 = P[P.year == 2015].set_index("iso")["pop"]
keep = [i for i in P.iso.unique() if pop15.get(i, 0) >= 1e6 and i in TX.index]
P = P[P.iso.isin(keep)].copy()
full = P.groupby("iso").co2.apply(lambda s: s.notna().sum() == (Y1 - Y0 + 1) and (s > 0).all())
P = P[P.iso.isin(full[full].index)].copy()
for v in ["co2", "co2_power", "co2_industry", "co2_transport", "co2_buildings", "ch4_agri"]:
    P["ln_" + v] = np.log(P[v].where(P[v] > 0))
P["g_any"] = P.iso.map(first_any).fillna(0).astype(int)
P["name"] = P.iso.map(names)
P = P.sort_values(["iso", "year"]).reset_index(drop=True)
P.to_csv("out/panel.csv", index=False)
ADOPT = {i: int(g) for i, g in P.groupby("iso").g_any.first().items() if g > 0}
EUR = {"AUT", "BEL", "BGR", "HRV", "CYP", "CZE", "DNK", "EST", "FIN", "FRA", "DEU", "GRC", "HUN", "IRL", "ITA", "LVA", "LTU", "LUX", "MLT", "NLD", "POL", "PRT", "ROU", "SVK", "SVN",
       "ESP", "SWE", "GBR", "NOR", "ISL", "CHE", "ALB", "BIH", "MNE", "MKD", "SRB", "UKR", "MDA", "BLR", "TUR", "XKX"}
TRANSITION = {"ALB", "BGR", "BIH", "BLR", "CZE", "EST", "HRV", "HUN", "LTU", "LVA", "MDA", "MKD", "MNE", "POL", "ROU", "RUS", "SRB", "SVK", "SVN", "UKR",
              "ARM", "AZE", "GEO", "KAZ", "KGZ", "TJK", "TKM", "UZB", "MNG"}
gdp95 = P[P.year == 1995].set_index("iso").gdppc
R["N_countries"] = int(P.iso.nunique()); R["N_never"] = int((P.groupby("iso").g_any.first() == 0).sum()); R["years"] = [Y0, Y1]
R["adopters"] = ADOPT
R["price2025_adopters"] = {i: round(float(price2025.get(i, np.nan)), 2) for i in ADOPT}
R["price_mean_0_10"] = {i: round(float(TX.loc[i, g:min(g + EHI, 2025)].mean()), 2) for i, g in ADOPT.items()}
R["desc"] = {"co2_2025_total_Mt": round(float(P[P.year == 2025].co2.sum()), 1),
             "share_adopters_2025": round(float(P[(P.year == 2025) & (P.g_any > 0)].co2.sum() / P[P.year == 2025].co2.sum()), 4),
             "adopter_growth_1990_2025": round(float(np.log(P[(P.year == 2025) & (P.g_any > 0)].co2.sum() / P[(P.year == 1990) & (P.g_any > 0)].co2.sum())), 4),
             "never_growth_1990_2025": round(float(np.log(P[(P.year == 2025) & (P.g_any == 0)].co2.sum() / P[(P.year == 1990) & (P.g_any == 0)].co2.sum())), 4)}
pre_growth = {}
for i, g in ADOPT.items():
    s = P[(P.iso == i)].set_index("year").ln_co2
    if g - 11 >= Y0: pre_growth[i] = float((s[g - 1] - s[g - 11]) / 10)
never_growth = P[P.g_any == 0].groupby("iso").apply(lambda d: d.set_index("year").ln_co2.diff().mean())
R["pre_growth_adopters_mean"] = round(float(np.mean(list(pre_growth.values()))), 4)
R["growth_never_mean_1980_2025"] = round(float(never_growth.mean()), 4)
log("panel", R["N_countries"], "countries,", len(ADOPT), "adopters")

# =============================== estimators ===============================
# All estimators take (df, y) with columns iso, year, y, g_any and return an event-time Series (index e).
def wide(df, y):
    W = df.pivot(index="iso", columns="year", values=y); G = df.groupby("iso").g_any.first().reindex(W.index).fillna(0).astype(int)
    return W, G

def cs_attgt(df, y, control="never"):
    """Callaway & Sant'Anna (2021) ATT(g,t), unconditional, universal base period g-1."""
    W, G = wide(df, y); yrs = W.columns.values; out = []
    for g in sorted(G[G > 0].unique()):
        if g - 1 < yrs.min(): continue
        tr = W.loc[G == g]
        for t in yrs:
            e = t - g
            if e < EMIN or e > EHI or t == g - 1: continue
            cn = W.loc[G == 0] if control == "never" else W.loc[(G == 0) | (G > max(t, g - 1))]
            dt = (tr[t] - tr[g - 1]).dropna(); dc = (cn[t] - cn[g - 1]).dropna()
            if len(dt) == 0 or len(dc) < 5: continue
            out.append((g, t, e, len(dt), float(dt.mean() - dc.mean())))
    return pd.DataFrame(out, columns=["g", "t", "e", "n_g", "att"])

def agg(a):
    return a.groupby("e").apply(lambda d: np.average(d.att, weights=d.n_g)).rename("att")

def est_cs_never(df, y): return agg(cs_attgt(df, y, "never"))
def est_cs_notyet(df, y): return agg(cs_attgt(df, y, "notyet"))
def est_cs_comparable(df, y):
    """Callaway-Sant'Anna with the comparison pool restricted to countries with GDP per capita (PPP) of at least 15,000 in 1995."""
    rich = set(gdp95[gdp95 >= 15000].index)
    base = df.iso.str.split("_").str[0]
    return agg(cs_attgt(df[base.isin(rich) | (df.g_any > 0)], y, "never"))

def est_cs_trend(df, y, K=10):
    """Cohort-specific linear trends: the treated-minus-never-treated gap of each cohort is net of a linear trend fitted on its K pre-adoption years."""
    W, G = wide(df, y); yrs = W.columns.values; out = []
    co = W.loc[G == 0]
    for g in sorted(G[G > 0].unique()):
        tr = W.loc[G == g]; Dg = tr.mean() - co.mean(); fy = [t for t in yrs if g - K <= t <= g - 1 and np.isfinite(Dg.get(t, np.nan))]
        if len(fy) < K: continue
        b = np.polyfit(fy, Dg[fy].values, 1)
        for t in yrs:
            e = t - g
            if EMIN <= e <= EHI and np.isfinite(Dg[t]): out.append((g, t, e, len(tr), float(Dg[t] - np.polyval(b, t))))
    return agg(pd.DataFrame(out, columns=["g", "t", "e", "n_g", "att"]))

def est_twfe(df, y):
    """Dynamic two-way fixed effects event study (binned endpoints), estimated by within transformation."""
    d = df.dropna(subset=[y]).copy(); g = d.g_any
    rel = np.where(g > 0, (d.year - g).clip(EMIN, EHI), -1).astype(int)
    ks = [k for k in range(EMIN, EHI + 1) if k != -1]
    X = np.column_stack([(rel == k).astype(float) for k in ks]); yv = d[y].values
    def dm(a):
        a = pd.DataFrame(a); a = a - a.groupby(d.iso.values).transform("mean"); a = a - a.groupby(d.year.values).transform("mean")
        return a.values
    Xd, yd = dm(X), dm(yv[:, None])[:, 0]
    for _ in range(20):          # alternating projections for the unbalanced case
        Xd, yd = dm(Xd), dm(yd[:, None])[:, 0]
    b = np.linalg.lstsq(Xd, yd, rcond=None)[0]
    return pd.Series(dict(zip(ks, b))).sort_index()

def fe_fit(Y, M, iters=300):
    Ym = np.where(M, Y, np.nan); mu = np.nanmean(Ym); a = np.zeros(Y.shape[0]); b = np.zeros(Y.shape[1])
    for _ in range(iters):
        a2 = np.nan_to_num(np.nanmean(Ym - mu - b[None, :], axis=1)); b2 = np.nan_to_num(np.nanmean(Ym - mu - a2[:, None], axis=0))
        d = np.max(np.abs(a2 - a)) + np.max(np.abs(b2 - b)); a, b = a2, b2
        if d < 1e-10: break
    return mu + a[:, None] + b[None, :]

def fe_full(Z):
    mu = Z.mean(); return mu + (Z.mean(1) - mu)[:, None] + (Z.mean(0) - mu)[None, :]

def cf_fit(Y, M, r=0, lam=None, iters=1000, tol=1e-8):
    """Counterfactual Y(0) from untreated cells (M=True): two-way FE (r=0), interactive FE with r factors (Xu, 2017; Liu, Wang & Xu, 2024)
    or matrix completion with nuclear-norm penalty lam (Athey et al., 2021), by iterative imputation and (soft-)thresholded SVD."""
    M = M & np.isfinite(Y); fit = fe_fit(Y, M); Z = np.where(M, Y, fit)
    if r == 0 and lam is None: return fit
    for it in range(iters):
        fe = fe_full(Z); U, s, Vt = np.linalg.svd(Z - fe, full_matrices=False)
        s2 = np.maximum(s - lam, 0) if lam is not None else np.where(np.arange(len(s)) < r, s, 0)
        new = fe + (U * s2) @ Vt; Znew = np.where(M, Y, new)
        ch = np.sqrt(np.mean((Znew - Z) ** 2)); Z = Znew
        if ch < tol: break
    return new

def cf_event(df, y, r=0, lam=None):
    W, G = wide(df, y); Y = W.values.astype(float); yrs = W.columns.values; Gv = G.values
    D = (Gv[:, None] > 0) & (yrs[None, :] >= Gv[:, None]); fit = cf_fit(Y, ~D, r, lam)
    E = yrs[None, :] - Gv[:, None]; T = Gv > 0; res = {}
    for e in range(EMIN, EHI + 1):
        m = (E == e) & T[:, None] & np.isfinite(Y)
        if m.any(): res[e] = float(np.mean((Y - fit)[m]))
    return pd.Series(res), pd.Series({W.index[i]: float(np.nanmean(np.where((E[i] >= 0) & (E[i] <= EHI), Y[i] - fit[i], np.nan))) for i in np.where(T)[0]})

def block_cv(df, y, cands, mode="r", nobs=3, k=10, seed=7):
    """Cross-validation in the spirit of Liu, Wang & Xu (2024): hold out blocks of consecutive pre-adoption observations of treated units."""
    W, G = wide(df, y); Y = W.values.astype(float); yrs = W.columns.values; Gv = G.values
    M = ~((Gv[:, None] > 0) & (yrs[None, :] >= Gv[:, None])); rr = np.random.default_rng(seed); err = {c: [] for c in cands}
    for f in range(k):
        M2 = M.copy(); held = []
        for i in np.where(Gv > 0)[0]:
            pre = np.where(M[i])[0]
            if len(pre) < nobs + 5: continue
            s = rr.integers(5, len(pre) - nobs + 1); cells = pre[s:s + nobs]; M2[i, cells] = False; held += [(i, j) for j in cells]
        h = np.array(held)
        for c in cands:
            fit = cf_fit(Y, M2, r=c if mode == "r" else 0, lam=c if mode == "lam" else None)
            err[c].append(float(np.mean((Y[h[:, 0], h[:, 1]] - fit[h[:, 0], h[:, 1]]) ** 2)))
    return {c: float(np.mean(v)) for c, v in err.items()}

def fw(Yw, zeta, x=None, min_decrease=1e-3, max_iter=10000):
    """Frank-Wolfe solver for simplex-constrained ridge weights with intercept, as in the synthdid reference implementation."""
    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); A = Yw[:, :T0]; b = Yw[:, T0]; eta = N * zeta ** 2; Ax = A @ x; prev = None
    for t in range(max_iter):
        hg = (Ax - b) @ A + eta * x; i = int(np.argmin(hg)); dx = -x.copy(); dx[i] += 1.0
        if not dx.any(): break
        derr = A[:, 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 - b; 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_unit_time(co_pre, co_post, tr_pre, tr_post):
    N1, T1 = tr_post.shape; noise = np.std(np.diff(co_pre, axis=1), ddof=1)
    z_om = (N1 * T1) ** 0.25 * noise; z_la = 1e-6 * noise; md = 1e-5 * noise
    Yl = np.column_stack([co_pre, co_post.mean(1)]); lam = fw(Yl, z_la, min_decrease=md, max_iter=100); lam = fw(Yl, z_la, x=sparsify(lam), min_decrease=md)
    Yo = np.column_stack([co_pre.T, tr_pre.mean(0)]); om = fw(Yo, z_om, min_decrease=md, max_iter=100); om = fw(Yo, z_om, x=sparsify(om), min_decrease=md)
    return om, lam

def sdid_block(W, trmask, comask, g, ehi):
    yrs = W.columns.values; pre = [t for t in yrs if t < g]; post = [t for t in yrs if g <= t <= g + ehi]
    cols = [t for t in yrs if t <= g + ehi]
    tr = W.loc[trmask, cols].dropna(); co = W.loc[comask, cols].dropna()
    if len(pre) < 5 or len(post) == 0 or len(tr) == 0: return None
    om, lam = sdid_unit_time(co[pre].values, co[post].values, tr[pre].values, tr[post].values)
    trp = tr[pre].values.mean(0) @ lam; cop = om @ (co[pre].values @ lam)
    ev = {int(t - g): float((tr[t].mean() - trp) - (om @ co[t].values - cop)) for t in cols if t >= g + EMIN}
    return {"ev": ev, "n": len(tr), "T1": len(post), "tau": float(np.mean([ev[e] for e in range(0, ehi + 1) if e in ev])), "omega": pd.Series(om, index=co.index)}

def est_sdid(df, y, by_unit=False, return_cohorts=False):
    """Synthetic difference-in-differences (Arkhangelsky et al., 2021), estimated cohort by cohort against never-treated countries and
    aggregated with weights proportional to treated units x post periods (Clarke et al., 2023)."""
    W, G = wide(df, y); res = {}
    groups = [(i, W.index == i, int(G[i])) for i in G[G > 0].index] if by_unit else [(g, (G == g).values, int(g)) for g in sorted(G[G > 0].unique())]
    for key, trmask, g in groups:
        ehi = min(EHI, int(df[df.g_any == g].year.max()) - g) if (df.g_any == g).any() else EHI
        r = sdid_block(W, trmask, (G == 0).values, g, ehi)
        if r: res[key] = r
    if return_cohorts: return res
    evs = {}
    for r in res.values():
        for e, v in r["ev"].items(): evs.setdefault(e, []).append((v, r["n"]))
    return pd.Series({e: np.average([a for a, _ in l], weights=[n for _, n in l]) for e, l in evs.items()}).sort_index()

def est_lpdid(df, y, L=4):
    """Local-projections difference-in-differences (Dube, Girardi, Jorda & Taylor, 2025) with clean (not-yet-treated) controls,
    year effects and L lags of the outcome growth to absorb pre-existing trajectories."""
    W, G = wide(df, y); Y = W.values.astype(float); yrs = W.columns.values; Gv = G.values; N, T = Y.shape; res = {}
    dY = np.full_like(Y, np.nan); dY[:, 1:] = Y[:, 1:] - Y[:, :-1]
    for h in range(0, EHI + 1):
        rows = []
        for j in range(L + 1, T - h):
            t = yrs[j]
            treat = Gv == t; clean = (Gv == 0) | (Gv > t + h)
            use = (treat | clean) & np.isfinite(Y[:, j + h]) & np.isfinite(Y[:, j - 1])
            lags = np.column_stack([dY[:, j - l] for l in range(1, L + 1)]); use &= np.isfinite(lags).all(1)
            for i in np.where(use)[0]: rows.append((j, float(treat[i]), Y[i, j + h] - Y[i, j - 1], *lags[i]))
        if not rows: continue
        A = np.array(rows); tt = A[:, 0]
        if A[:, 1].sum() == 0: continue
        Z = A[:, 1:]
        Zd = Z - pd.DataFrame(Z).groupby(tt).transform("mean").values
        Xd = np.column_stack([Zd[:, 0]] + [Zd[:, 2 + l] for l in range(L)]); yd = Zd[:, 1]
        res[h] = float(np.linalg.lstsq(Xd, yd, rcond=None)[0][0])
    return pd.Series(res)

R_IFE = {}; LAM_MC = None
def est_ife(df, y): return cf_event(df, y, r=R_IFE[y])[0]
def est_mc(df, y): return cf_event(df, y, lam=LAM_MC)[0]
def est_bjs(df, y): return cf_event(df, y, r=0)[0]

ESTIMATORS = {
    "twfe": ("Two-way fixed effects (dynamic)", est_twfe, "naive"),
    "cs_never": ("Callaway–Sant'Anna, never treated", est_cs_never, "naive"),
    "cs_notyet": ("Callaway–Sant'Anna, not yet treated", est_cs_notyet, "naive"),
    "bjs": ("Imputation (Borusyak et al.)", est_bjs, "naive"),
    "cs_comparable": ("Callaway–Sant'Anna, comparable-income controls", est_cs_comparable, "naive"),
    "lpdid": ("LP-DiD with lagged outcome growth", est_lpdid, "trend"),
    "cs_trend": ("Callaway–Sant'Anna, cohort-specific trends", est_cs_trend, "trend"),
    "sdid": ("Synthetic difference-in-differences", est_sdid, "trend"),
    "ife": ("Interactive fixed effects counterfactual", est_ife, "trend"),
    "mc": ("Matrix completion counterfactual", est_mc, "trend"),
}

# =============================== inference and falsification ===============================
def att(ev, lo=0, hi=EHI):
    v = ev.loc[(ev.index >= lo) & (ev.index <= hi)]; return float(v.mean()) if len(v) else np.nan

def shift_placebo(df, s=SHIFT):
    """Placebo design: each adopter's adoption is moved s years earlier and its observations from the true adoption year onward are dropped."""
    d = df.copy(); tr = d.g_any > 0
    d = d[~(tr & (d.year >= d.g_any))].copy(); d.loc[d.g_any > 0, "g_any"] = d.loc[d.g_any > 0, "g_any"] - s
    return d

def boot(df, fn, B):
    isos = df.iso.unique(); groups = {i: g for i, g in df.groupby("iso")}; out = []
    for b in range(B):
        s = rng.choice(isos, len(isos), replace=True)
        bd = pd.concat([groups[i].assign(iso=f"{i}_{k}") for k, i in enumerate(s)], ignore_index=True)
        try: out.append(fn(bd))
        except Exception as ex: pass
    return out

def summarize(point, draws):
    a = np.array([d for d in draws if np.isfinite(d)])
    return {"est": round(point, 4), "se": round(float(a.std(ddof=1)), 4), "lo": round(float(np.quantile(a, 0.025)), 4), "hi": round(float(np.quantile(a, 0.975)), 4),
            "pct": round(100 * (np.exp(point) - 1), 2), "B": int(len(a))}

def run_estimator(key, df, y, B, keep_event=False):
    lab, fn, fam = ESTIMATORS[key]
    ev = fn(df, y); dfp = shift_placebo(df); evp = fn(dfp, y)
    def stat(d):
        e = fn(d, y); p = fn(shift_placebo(d), y)
        return e, att(e), att(e, 0, 5), att(p, 0, SHIFT - 1)
    dr = boot(df[["iso", "year", y, "g_any"]].dropna(), stat, B)
    out = {"label": lab, "family": fam, "att_0_10": summarize(att(ev), [d[1] for d in dr]), "att_0_5": summarize(att(ev, 0, 5), [d[2] for d in dr]),
           "placebo": summarize(att(evp, 0, SHIFT - 1), [d[3] for d in dr])}
    out["placebo_pass"] = bool(out["placebo"]["lo"] <= 0 <= out["placebo"]["hi"])
    if keep_event:
        E = pd.DataFrame([d[0] for d in dr]); out["event"] = pd.DataFrame({"est": ev, "lo": E.quantile(0.025).reindex(ev.index), "hi": E.quantile(0.975).reindex(ev.index)})
    return out

# =============================== tuning of the factor models ===============================
CV_ALL = {}
for yv in ["ln_co2", "ln_co2_power", "ln_co2_industry", "ln_co2_transport", "ln_co2_buildings", "ln_ch4_agri"]:
    dd = P.dropna(subset=[yv]); dd = dd[dd.groupby("iso")[yv].transform("size") == dd.groupby("iso").year.transform("size")]
    CV_ALL[yv] = block_cv(dd, yv, [0, 1, 2, 3, 4, 5]); R_IFE[yv] = int(min(CV_ALL[yv], key=CV_ALL[yv].get))
cv_r = CV_ALL["ln_co2"]
W0, G0 = wide(P, "ln_co2"); Yv = W0.values; D0 = (G0.values[:, None] > 0) & (W0.columns.values[None, :] >= G0.values[:, None])
smax = np.linalg.svd(np.nan_to_num(Yv - fe_fit(Yv, ~D0)), compute_uv=False)[0]
lams = [round(smax * f, 4) for f in (0.3, 0.2, 0.1, 0.05, 0.03, 0.02, 0.01)]
cv_l = block_cv(P, "ln_co2", lams, mode="lam"); LAM_MC = float(min(cv_l, key=cv_l.get))
R["cv"] = {"ife_mspe": {str(k): round(v, 5) for k, v in cv_r.items()}, "r_selected": R_IFE["ln_co2"], "r_by_outcome": R_IFE,
           "mc_mspe": {str(k): round(v, 5) for k, v in cv_l.items()}, "lambda_selected": LAM_MC, "lambda_over_smax": round(LAM_MC / smax, 3)}
log("CV: r =", R_IFE, " lambda =", LAM_MC)

# =============================== main comparison ===============================
MAIN = {}
for k in ESTIMATORS:
    MAIN[k] = run_estimator(k, P, "ln_co2", B_MAIN, keep_event=True)
    MAIN[k]["event"].to_csv(f"out/event_{k}.csv")
    log(k, MAIN[k]["att_0_10"], "placebo", MAIN[k]["placebo"])
R["main"] = {k: {kk: vv for kk, vv in v.items() if kk != "event"} for k, v in MAIN.items()}
credible = [k for k, v in MAIN.items() if v["placebo_pass"]]
R["credible"] = credible

# smoothness-restriction sensitivity (Rambachan & Roth, 2023) for the cohort-specific-trend estimator: worst-case bias of the linear
# extrapolation when second differences of the differential trend are bounded by M is M (e+1)(e+2)/2 at horizon e.
mt = MAIN["cs_trend"]["att_0_10"]; bfac = np.mean([(e + 1) * (e + 2) / 2 for e in range(0, EHI + 1)])
R["smoothness"] = {str(M): {"lo": round(mt["est"] - M * bfac - 1.96 * mt["se"], 4), "hi": round(mt["est"] + M * bfac + 1.96 * mt["se"], 4)} for M in (0, 0.0005, 0.001, 0.002, 0.004)}
R["smoothness_bias_factor"] = round(float(bfac), 2)

# =============================== heterogeneity, sectors and robustness (trend-robust estimators) ===============================
HET_EST = ["cs_trend", "ife"]
hi_iso = [i for i in ADOPT if (price2025.get(i, 0) or 0) >= 20]; lo_iso = [i for i in ADOPT if (price2025.get(i, 0) or 0) < 20]
R["high_price"] = hi_iso; R["low_price"] = lo_iso
SPECS = {
    "high_price": (P[~P.iso.isin(lo_iso)], "ln_co2"),
    "low_price": (P[~P.iso.isin(hi_iso)], "ln_co2"),
    "excl_transition": (P[~P.iso.isin(TRANSITION)], "ln_co2"),
    "europe": (P[P.iso.isin(EUR)], "ln_co2"),
    "since_1990": (P[P.year >= 1990], "ln_co2"),
    "power": (P, "ln_co2_power"), "industry": (P, "ln_co2_industry"), "transport": (P, "ln_co2_transport"), "buildings": (P, "ln_co2_buildings"),
    "placebo_ch4_agri": (P, "ln_ch4_agri"),
}
HET = {}
for s, (d, y) in SPECS.items():
    d = d.dropna(subset=[y]); d = d[d.groupby("iso")[y].transform("size") == d.groupby("iso").year.transform("size")]
    HET[s] = {}
    for k in HET_EST:
        r = run_estimator(k, d, y, B_HET); HET[s][k] = {kk: vv for kk, vv in r.items() if kk != "event"}
    log("het", s, {k: (HET[s][k]["att_0_10"]["est"], HET[s][k]["placebo"]["est"]) for k in HET_EST})
R["het"] = HET

# =============================== country-level effects and price ===============================
coh = est_sdid(P, "ln_co2", by_unit=True, return_cohorts=True)
_, ife_unit = cf_event(P, "ln_co2", r=R_IFE["ln_co2"])
CU = pd.DataFrame({"sdid": {i: c["tau"] for i, c in coh.items()}, "ife": ife_unit}).dropna()
CU["price"] = [R["price_mean_0_10"][i] for i in CU.index]; CU["g"] = [ADOPT[i] for i in CU.index]
CU.to_csv("out/country_effects.csv")
def wls_slope(x, yv):
    X = np.column_stack([np.ones(len(x)), x]); b = np.linalg.lstsq(X, yv, rcond=None)[0]; return b
dose = {}
for k in ("sdid", "ife"):
    b = wls_slope(np.log1p(CU.price.values), CU[k].values)
    bd = []
    for _ in range(999):
        idx = rng.integers(0, len(CU), len(CU)); bd.append(wls_slope(np.log1p(CU.price.values[idx]), CU[k].values[idx])[1])
    dose[k] = {"intercept": round(float(b[0]), 4), "slope_log1p_price": round(float(b[1]), 4), "lo": round(float(np.quantile(bd, 0.025)), 4), "hi": round(float(np.quantile(bd, 0.975)), 4),
               "spearman": round(float(CU[[k, "price"]].corr(method="spearman").iloc[0, 1]), 3)}
R["dose"] = dose
R["country_effects"] = {i: {"sdid": round(float(r.sdid), 4), "ife": round(float(r.ife), 4), "price": float(r.price)} for i, r in CU.iterrows()}
sw = coh.get("SWE")
R["sdid_weights_example"] = {k: round(float(v), 3) for k, v in coh["SWE"]["omega"].sort_values(ascending=False).head(6).items()} if sw else None
json.dump(R, open("out/results.json", "w"), indent=1, default=str)
log("results written")

# =============================== figures ===============================
fig, ax = plt.subplots(figsize=(7.2, 3.9))
ad = pd.Series(ADOPT).sort_values(kind="stable")
for k, (i, g) in enumerate(ad.items()):
    s = TX.loc[i, g:2025]; ax.scatter(s.index, [k] * len(s), c=np.log1p(s.clip(lower=0).values), cmap="viridis", vmin=0, vmax=np.log1p(160), s=9, marker="s")
ax.set_yticks(range(len(ad))); ax.set_yticklabels([names.get(i, i) for i in ad.index], fontsize=6.5); ax.set_xlabel("Year"); ax.set_xlim(1988, 2026)
sm = plt.cm.ScalarMappable(cmap="viridis", norm=matplotlib.colors.Normalize(0, np.log1p(160))); cb = fig.colorbar(sm, ax=ax, pad=0.01)
cb.set_ticks(np.log1p([0, 5, 20, 50, 150])); cb.set_ticklabels(["0", "5", "20", "50", "150"]); cb.set_label("Emissions-weighted tax rate (2021 US$/tCO$_2$e)")
fig.tight_layout(); fig.savefig("fig/fig1_adoption.png"); plt.show()

fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.4), sharey=True)
pal = {"twfe": "#7f7f7f", "cs_never": "#1f77b4", "cs_notyet": "#17becf", "bjs": "#9467bd", "cs_comparable": "#8c564b",
       "lpdid": "#2ca02c", "cs_trend": "#bcbd22", "sdid": "#d62728", "ife": "#ff7f0e", "mc": "#e377c2"}
for ax, fam, ttl in [(axs[0], "naive", "(a) Parallel-trends estimators"), (axs[1], "trend", "(b) Trend-robust estimators")]:
    ks = [k for k in ESTIMATORS if ESTIMATORS[k][2] == fam]
    for j, k in enumerate(ks):
        d = MAIN[k]["event"]; off = (j - (len(ks) - 1) / 2) * 0.14
        ax.errorbar(d.index + off, d.est * 100, yerr=[((d.est - d.lo) * 100).clip(lower=0).fillna(0), ((d.hi - d.est) * 100).clip(lower=0).fillna(0)],
                    fmt="o", ms=2.2, lw=0.6, color=pal[k], label=ESTIMATORS[k][0])
    ax.axhline(0, color="k", lw=0.6); ax.axvline(-0.5, color="k", lw=0.6, ls=":"); ax.set_title(ttl, loc="left"); ax.set_xlabel("Years relative to adoption")
    ax.legend(frameon=False, fontsize=6, loc="lower left")
axs[0].set_ylabel("Effect on CO$_2$ (log points × 100)"); axs[0].set_ylim(-95, 45)
fig.tight_layout(); fig.savefig("fig/fig2_event.png"); plt.show()

fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.4), sharey=True, gridspec_kw={"width_ratios": [1.3, 1]})
ks = list(ESTIMATORS)
for j, k in enumerate(ks):
    m = MAIN[k]; c = "#d62728" if m["placebo_pass"] else "#7f7f7f"
    for ax, key in [(axs[0], "att_0_10"), (axs[1], "placebo")]:
        r = m[key]; ax.plot([r["lo"] * 100, r["hi"] * 100], [j, j], color=c, lw=1.2); ax.plot(r["est"] * 100, j, "o", color=c, ms=4)
for ax, ttl in [(axs[0], "(a) Effect, years 0–10"), (axs[1], "(b) Placebo, 3 years early")]:
    ax.axvline(0, color="k", lw=0.6); ax.set_title(ttl, loc="left"); ax.set_xlabel("Log points × 100")
axs[0].set_yticks(range(len(ks))); axs[0].set_yticklabels([ESTIMATORS[k][0] for k in ks], fontsize=7); axs[0].invert_yaxis()
fig.tight_layout(); fig.savefig("fig/fig3_estimators.png"); plt.show()

fig, ax = plt.subplots(figsize=(7.2, 4.2))
labels = {"high_price": "Tax ≥ US$20 in 2025", "low_price": "Tax < US$20 in 2025", "excl_transition": "Excluding transition economies", "europe": "European countries only",
          "since_1990": "Sample 1990–2025", "power": "Power sector CO$_2$", "industry": "Industry CO$_2$", "transport": "Transport CO$_2$", "buildings": "Buildings CO$_2$",
          "placebo_ch4_agri": "Placebo outcome: agricultural CH$_4$"}
rows = [("All adopters, total CO$_2$", {k: R["main"][k] for k in HET_EST})] + [(labels[s], HET[s]) for s in SPECS]
mk = {"cs_trend": ("o", "#bcbd22", "Cohort-specific trends"), "ife": ("s", "#ff7f0e", "Interactive fixed effects")}
for j, (lab, d) in enumerate(rows):
    for q, k in enumerate(HET_EST):
        r = d[k]["att_0_10"]; yy = j + (q - 0.5) * 0.26
        if not d[k].get("placebo_pass", True): ax.plot(r["est"] * 100, yy, "x", color="k", ms=5, mew=0.8)
        ax.plot([r["lo"] * 100, r["hi"] * 100], [yy, yy], color=mk[k][1], lw=0.9); ax.plot(r["est"] * 100, yy, mk[k][0], color=mk[k][1], ms=3.5, label=mk[k][2] if j == 0 else None)
ax.axvline(0, color="k", lw=0.6); ax.set_yticks(range(len(rows))); ax.set_yticklabels([r[0] for r in rows], fontsize=7.5); ax.invert_yaxis()
ax.set_xlabel("Average effect over years 0–10 (log points × 100), 95% bootstrap CI"); ax.legend(frameon=False, loc="lower right")
fig.tight_layout(); fig.savefig("fig/fig4_heterogeneity.png"); plt.show()

fig, axs = plt.subplots(1, 2, figsize=(7.5, 3.2), sharey=True)
for ax, k, ttl in [(axs[0], "sdid", "(a) SDID, country by country"), (axs[1], "ife", "(b) Interactive fixed effects")]:
    x = np.log1p(CU.price.values); ax.scatter(x, CU[k] * 100, s=14, color="#1f77b4")
    for i, r in CU.iterrows(): ax.annotate(i, (np.log1p(r.price), r[k] * 100), fontsize=5.5, xytext=(2, 2), textcoords="offset points")
    xs = np.linspace(0, x.max(), 50); ax.plot(xs, (dose[k]["intercept"] + dose[k]["slope_log1p_price"] * xs) * 100, color="#d62728", lw=1)
    ax.axhline(0, color="k", lw=0.5); ax.set_title(ttl, loc="left")
    ax.set_xticks(np.log1p([0, 5, 20, 50, 100])); ax.set_xticklabels(["0", "5", "20", "50", "100"]); ax.set_xlabel("Mean tax rate in years 0–10 (2021 US$/t, log scale)")
axs[0].set_ylabel("Country effect, years 0–10 (log points × 100)")
fig.tight_layout(); fig.savefig("fig/fig5_dose.png"); plt.show()
log("done")
print(json.dumps({k: R[k] for k in ["N_countries", "N_never", "cv", "credible", "smoothness", "dose", "desc", "pre_growth_adopters_mean", "growth_never_mean_1980_2025"]}, 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/rp06_outputs", "zip", ".", ".")
files.download("/content/rp06_outputs.zip")