# -*- coding: utf-8 -*-
"""Did the energy crisis drive Europe's bankruptcy wave? Sector exposure, local projections and the country dimension of business failures

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

# # Did the energy crisis drive Europe's bankruptcy wave? Sector exposure, local projections and the country dimension of business failures
#
# 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: 1–3 dakika.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: Eurostat iflas ve kayıt endeksleri, üretici fiyatları, enerji hesapları ve milli hesaplar). 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/rp15", exist_ok=True); os.chdir("/content/rp15")
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 energy-price and bankruptcy paper into data/ from the Eurostat API:
    quarterly business registration and bankruptcy indices by NACE section (sts_rb_q, 2015Q1-2026Q3), monthly producer price indices
    for energy (sts_inpp_m, MIG_NRG, domestic market), energy use by NACE activity (env_ac_pefasu) and gross value added by
    industry (nama_10_a64) for the energy intensity of sectors, and half-yearly electricity and gas prices for non-household consumers."""
    import pandas as pd
    from common import eurostat, DATA
    if __name__ == "__main__":
        eurostat("sts_rb_q", "unit=I21&s_adj=NSA").to_csv(DATA / "sts_rb_q.csv", index=False)
        eurostat("sts_rb_q", "unit=I21&s_adj=SCA").to_csv(DATA / "sts_rb_q_sca.csv", index=False)
        eurostat("sts_inpp_m", "unit=I21&nace_r2=MIG_NRG&indic_bt=PRC_PRR&s_adj=NSA").to_csv(DATA / "ppi_energy_m.csv", index=False)
        eurostat("sts_inpp_m", "unit=I21&nace_r2=B-E36&indic_bt=PRC_PRR&s_adj=NSA").to_csv(DATA / "ppi_industry_m.csv", index=False)
        secs = "&".join(f"nace_r2={s}" for s in ["B", "C", "D", "E", "F", "G", "H", "I", "J", "K", "L", "M", "N", "P", "Q", "R", "S", "TOTAL"])
        eurostat("env_ac_pefasu", f"stk_flow=ER_USE&prod_nrg=P00&unit=TJ&{secs}").to_csv(DATA / "energy_use_nace.csv", index=False)
        eurostat("nama_10_a64", f"unit=CP_MEUR&na_item=B1G&{secs}").to_csv(DATA / "gva_nace.csv", index=False)
        eurostat("nrg_pc_205", "nrg_cons=MWH2000-19999&tax=X_VAT&currency=EUR&unit=KWH").to_csv(DATA / "elec_price_nonhh.csv", index=False)
        eurostat("nrg_pc_203", "nrg_cons=GJ10000-99999&tax=X_VAT&currency=EUR&unit=GJ_GCV").to_csv(DATA / "gas_price_nonhh.csv", index=False)
        print("done")
print(sorted(os.listdir('data')))

# ## 2. Yerel projeksiyonlar, sağlamlık testleri, ülke ayrıştırması ve şekiller

# Paper 15: energy price shocks and business dynamics in Europe. Quarterly bankruptcy declarations and business registrations by sector and country (Eurostat short-term business statistics, 2015Q1-2026Q2), energy intensity of sectors from physical energy flow accounts and national accounts, and country energy producer prices. Local projections with a shift-share exposure (sector energy intensity x country energy price shock), country-by-quarter and sector-by-quarter fixed effects, asymmetry and state dependence, and the contribution of the 2021-2023 energy price shock to bankruptcies.
import json, sys, warnings
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})
R = {}
AGG = {"EU27_2020", "EA19", "EA20", "EA21"}
SECTORS = {"B-E": ["B", "C", "D", "E"], "F": ["F"], "G": ["G"], "H": ["H"], "I": ["I"], "J": ["J"], "K-N": ["K", "L", "M", "N"], "P-S_X_S94": ["P", "Q", "R", "S"]}
SLAB = {"B-E": "Industry (B–E)", "F": "Construction (F)", "G": "Trade (G)", "H": "Transport (H)", "I": "Accommodation and food (I)",
        "J": "Information and communication (J)", "K-N": "Finance, real estate and business services (K–N)", "P-S_X_S94": "Education, health and other services (P–S)"}
LAST = "2026-Q2"

# ---------------- outcomes: bankruptcy and registration indices (2021 = 100, seasonally and calendar adjusted) ----------------
S = pd.read_csv("data/sts_rb_q_sca.csv")
S = S[S.nace_r2.isin(SECTORS) & ~S.geo.isin(AGG) & (S.time <= LAST)]
S["q"] = pd.PeriodIndex(S.time.str.replace("-", ""), freq="Q")
Y = S.pivot_table(index=["geo", "nace_r2", "q"], columns="indic_bt", values="value").reset_index()
Y.columns.name = None
# keep country-sector series that are strictly positive and complete from 2015Q1 for bankruptcies
ok = Y.groupby(["geo", "nace_r2"]).BKRT.apply(lambda s: (s > 0).all() and s.notna().sum() >= 40)
Y = Y.merge(ok.rename("ok").reset_index(), on=["geo", "nace_r2"]); Y = Y[Y.ok].drop(columns="ok")
Y["lb"] = np.log(Y.BKRT); Y["lr"] = np.log(Y.REG.where(Y.REG > 0))

# ---------------- exposure: energy intensity (TJ per million EUR of gross value added), 2015-2019 average ----------------
E = pd.read_csv("data/energy_use_nace.csv"); G = pd.read_csv("data/gva_nace.csv")
E = E[(E.time >= 2015) & (E.time <= 2019)]; G = G[(G.time >= 2015) & (G.time <= 2019)]
rows = []
for sec, parts in SECTORS.items():
    e = E[E.nace_r2.isin(parts)].groupby(["geo", "time"]).value.sum(min_count=len(parts)); g = G[G.nace_r2.isin(parts)].groupby(["geo", "time"]).value.sum(min_count=len(parts))
    j = pd.concat([e.rename("tj"), g.rename("gva")], axis=1).dropna().groupby("geo").sum()
    j["nace_r2"] = sec; rows.append(j.reset_index())
EI = pd.concat(rows); EI["ei"] = EI.tj / EI.gva
EI = EI[EI.ei > 0]
# EU-wide sector intensity (used for the robustness check)
EU = EI.groupby("nace_r2")[["tj", "gva"]].sum(); EU["ei_eu"] = EU.tj / EU.gva
R["energy_intensity_eu_tj_per_meur"] = {SLAB[k]: round(float(v), 2) for k, v in EU.ei_eu.items()}

# ---------------- shock: country energy producer price index, quarterly mean, 4-quarter log change ----------------
P = pd.read_csv("data/ppi_energy_m.csv"); P = P[~P.geo.isin(AGG)]
P["q"] = pd.PeriodIndex(pd.to_datetime(P.time), freq="Q")
PQ = P.groupby(["geo", "q"]).value.mean().reset_index().sort_values(["geo", "q"])
PQ["lp"] = np.log(PQ.value); PQ["shock"] = PQ.groupby("geo").lp.diff(4)
PI = pd.read_csv("data/ppi_industry_m.csv"); PI = PI[~PI.geo.isin(AGG)]
PI["q"] = pd.PeriodIndex(pd.to_datetime(PI.time), freq="Q")
PIQ = PI.groupby(["geo", "q"]).value.mean().reset_index().sort_values(["geo", "q"])
PIQ["lpi"] = np.log(PIQ.value); PIQ["shock_ind"] = PIQ.groupby("geo").lpi.diff(4)

D = Y.merge(EI[["geo", "nace_r2", "ei"]], on=["geo", "nace_r2"]).merge(EU[["ei_eu"]].reset_index(), on="nace_r2").merge(PQ[["geo", "q", "shock", "lp"]], on=["geo", "q"])
D = D.merge(PIQ[["geo", "q", "shock_ind"]], on=["geo", "q"], how="left")
D = D.sort_values(["geo", "nace_r2", "q"]).reset_index(drop=True)
# standardise exposure: log energy intensity relative to its mean across country-sectors, scaled by its sd
cs = D.drop_duplicates(["geo", "nace_r2"])
D["x"] = (np.log(D.ei) - np.log(cs.ei).mean()) / np.log(cs.ei).std()
D["x_eu"] = (np.log(D.ei_eu) - np.log(cs.ei_eu).mean()) / np.log(cs.ei_eu).std()
D["z"] = D.x * D.shock; D["z_eu"] = D.x_eu * D.shock
D["zpos"] = D.x * D.shock.clip(lower=0); D["zneg"] = D.x * D.shock.clip(upper=0)
D["cq"] = D.geo + "_" + D.q.astype(str); D["sq"] = D.nace_r2 + "_" + D.q.astype(str); D["csid"] = D.geo + "_" + D.nace_r2
g_ = D.groupby("csid")
D["dlb_l1"] = g_.lb.diff().groupby(D.csid).shift(1); D["dlb_l2"] = g_.lb.diff().groupby(D.csid).shift(2)
D["z_l1"] = g_.z.shift(1); D["z_l2"] = g_.z.shift(2)
R["sample"] = {"countries": int(D.geo.nunique()), "sectors": int(D.nace_r2.nunique()), "country_sectors": int(D.csid.nunique()), "obs": int(len(D)),
               "first": str(D.q.min()), "last": str(D.q.max()), "countries_list": sorted(D.geo.unique().tolist())}
R["shock_sd"] = round(float(D.drop_duplicates(["geo", "q"]).shock.std()), 4)
R["shock_2022_mean"] = round(float(D[(D.q >= pd.Period("2022Q1")) & (D.q <= pd.Period("2022Q4"))].drop_duplicates(["geo", "q"]).shock.mean()), 4)
R["ei_by_sector_median"] = {SLAB[k]: round(float(v), 2) for k, v in cs.groupby("nace_r2").ei.median().items()}

# ---------------- local projections ----------------
HZ = list(range(0, 9))
def lp(y="lb", shock="z", extra="", data=None, fe="csid + cq + sq", hz=HZ, ctrl="dlb_l1 + dlb_l2 + z_l1 + z_l2"):
    data = D if data is None else data
    out = {}
    for h in hz:
        d = data.copy(); d["yh"] = d.groupby("csid")[y].shift(-h) - d.groupby("csid")[y].shift(1)
        d = d.dropna(subset=["yh", shock] + [c.strip() for c in ctrl.split("+") if c.strip()])
        fml = f"yh ~ {shock}{extra} + {ctrl} | {fe}"
        m = pf.feols(fml, data=d, vcov={"CRV1": "geo"}); t = m.tidy()
        out[h] = {k: {"b": round(float(t.loc[k, "Estimate"]), 4), "se": round(float(t.loc[k, "Std. Error"]), 4)} for k in [shock] + [e.strip() for e in extra.split("+") if e.strip()]}
        out[h]["nobs"] = int(m._N)
    return out
R["lp_bankruptcy"] = lp()
R["lp_registration"] = lp(y="lr", ctrl="z_l1 + z_l2")
R["lp_bankruptcy_eu_intensity"] = lp(shock="z_eu", ctrl="dlb_l1 + dlb_l2")
R["lp_asymmetry"] = lp(shock="zpos", extra=" + zneg", ctrl="dlb_l1 + dlb_l2")
# state dependence: monetary tightening (2022Q3-2024Q2, ECB deposit rate above 2%) versus other quarters
D["tight"] = ((D.q >= pd.Period("2022Q3")) & (D.q <= pd.Period("2024Q2"))).astype(int)
D["z_tight"] = D.z * D.tight; D["z_loose"] = D.z * (1 - D.tight)
R["lp_state"] = lp(shock="z_tight", extra=" + z_loose", ctrl="dlb_l1 + dlb_l2")
# placebo: industrial producer prices excluding energy are not available; use the industrial PPI as an alternative shock
D["z_ind"] = D.x * D.shock_ind
R["lp_industrial_ppi"] = lp(shock="z_ind", ctrl="dlb_l1 + dlb_l2", data=D.dropna(subset=["shock_ind"]))
R["lp_registration_eu_intensity"] = lp(y="lr", shock="z_eu", ctrl="z_l1 + z_l2")
R["lp_registration_excl_covid"] = lp(y="lr", ctrl="z_l1 + z_l2", data=D[(D.q < pd.Period("2020Q1")) | (D.q > pd.Period("2021Q4"))])
R["lp_registration_asymmetry"] = lp(y="lr", shock="zpos", extra=" + zneg", ctrl="z_l1 + z_l2")
# alternative shock: electricity price for non-household consumers (2,000-19,999 MWh, excl. VAT), half-yearly, 2-semester log change
EP = pd.read_csv("data/elec_price_nonhh.csv"); EP = EP[~EP.geo.isin(AGG) & (EP.value > 0)].copy()
EP["sem"] = EP.time.str[-1].astype(int); EP["yr"] = EP.time.str[:4].astype(int); EP = EP.sort_values(["geo", "yr", "sem"])
EP["le"] = np.log(EP.value); EP["dle"] = EP.groupby("geo").le.diff(2)
eq = []
for _, r_ in EP.iterrows():
    for qq in ((1, 2) if r_.sem == 1 else (3, 4)): eq.append((r_.geo, pd.Period(f"{r_.yr}Q{qq}"), r_.dle))
EQ = pd.DataFrame(eq, columns=["geo", "q", "shock_el"]).drop_duplicates(["geo", "q"])
D = D.merge(EQ, on=["geo", "q"], how="left"); D["z_el"] = D.x * D.shock_el
R["lp_bankruptcy_electricity"] = lp(shock="z_el", ctrl="dlb_l1 + dlb_l2", data=D.dropna(subset=["shock_el"]))
R["lp_registration_electricity"] = lp(y="lr", shock="z_el", ctrl="dlb_l1 + dlb_l2", data=D.dropna(subset=["shock_el", "lr"]))
# excluding the pandemic (2020Q1-2021Q4)
R["lp_excl_covid"] = lp(data=D[(D.q < pd.Period("2020Q1")) | (D.q > pd.Period("2021Q4"))])
# pre-period falsification: leads of the exposure (shock at t+4 on outcome change t-1..t)
d = D.copy(); d["z_lead4"] = d.groupby("csid").z.shift(-4); d["dy"] = d.groupby("csid").lb.diff()
m = pf.feols("dy ~ z_lead4 + dlb_l1 + dlb_l2 | csid + cq + sq", data=d.dropna(subset=["z_lead4", "dy", "dlb_l1", "dlb_l2"]), vcov={"CRV1": "geo"})
R["lead_placebo"] = {"b": round(float(m.coef()["z_lead4"]), 4), "se": round(float(m.se()["z_lead4"]), 4)}

# ---------------- catch-up: bankruptcies missing during the pandemic and the excess since 2022 ----------------
def qrange(a, b): return (D.q >= pd.Period(a)) & (D.q <= pd.Period(b))
D["t"] = (D.q - pd.Period("2015Q1")).apply(lambda x: x.n)
cf = {}
for k, g in D.groupby("csid"):
    pre = g[(g.q >= pd.Period("2017Q1")) & (g.q <= pd.Period("2019Q4"))]
    tr = g[(g.q >= pd.Period("2015Q1")) & (g.q <= pd.Period("2019Q4"))]
    b1, b0 = np.polyfit(tr.t, tr.lb, 1)
    cf[k] = (pre.lb.mean(), b0, b1)
D["cf_flat"] = D.csid.map(lambda k: cf[k][0]); D["cf_trend"] = D.csid.map(lambda k: cf[k][1]) + D.csid.map(lambda k: cf[k][2]) * D.t
for v in ["flat", "trend"]:
    D[f"dev_{v}"] = D.lb - D[f"cf_{v}"]
    D[f"lvl_gap_{v}"] = D.BKRT - np.exp(D[f"cf_{v}"])
WIN = {"2022": ("2022Q1", "2022Q4"), "2023": ("2023Q1", "2023Q4"), "2024": ("2024Q1", "2024Q4"), "2025": ("2025Q1", "2025Q4"), "2026H1": ("2026Q1", "2026Q2"), "2022-2026H1": ("2022Q1", "2026Q2")}
CU = {}
for v in ["flat", "trend"]:
    cs_ = pd.DataFrame({"short": -D[qrange("2020Q2", "2021Q4")].groupby("csid")[f"dev_{v}"].mean(),
                        "cum_short": -D[qrange("2020Q2", "2021Q4")].groupby("csid")[f"lvl_gap_{v}"].sum()})
    for w, (a, b) in WIN.items():
        cs_[f"ex_{w}"] = D[qrange(a, b)].groupby("csid")[f"dev_{v}"].mean()
        cs_[f"cum_ex_{w}"] = D[qrange(a, b)].groupby("csid")[f"lvl_gap_{v}"].sum()
    cs_["ratio"] = cs_["cum_ex_2022-2026H1"] / cs_.cum_short
    cs_["geo"] = cs_.index.str.split("_").str[0]; cs_["nace_r2"] = cs_.index.str.split("_", n=1).str[1]
    CU[v] = cs_
C0 = CU["flat"]
R["catchup"] = {v: {"median_shortfall_logpts": round(float(100 * CU[v].short.median()), 1),
                    "median_excess_logpts": {w: round(float(100 * CU[v][f"ex_{w}"].median()), 1) for w in WIN},
                    "share_csid_shortfall_positive": round(float((CU[v].short > 0).mean()), 3),
                    "aggregate_ratio": round(float(CU[v]["cum_ex_2022-2026H1"].sum() / CU[v].cum_short.sum()), 2),
                    "median_ratio_if_shortfall_positive": round(float(CU[v][CU[v].cum_short > 0].ratio.median()), 2),
                    "ratio_by_sector": {SLAB[s]: round(float(g["cum_ex_2022-2026H1"].sum() / g.cum_short.sum()), 2) for s, g in CU[v].groupby("nace_r2")},
                    "shortfall_by_sector_logpts": {SLAB[s]: round(float(100 * g.short.median()), 1) for s, g in CU[v].groupby("nace_r2")},
                    "excess_by_sector_logpts": {SLAB[s]: round(float(100 * g["ex_2022-2026H1"].median()), 1) for s, g in CU[v].groupby("nace_r2")}} for v in CU}
# cross-sectional regressions: excess bankruptcies on the pandemic shortfall and on energy-price exposure
lp2022 = PQ[(PQ.q >= pd.Period("2022Q1")) & (PQ.q <= pd.Period("2022Q4"))].groupby("geo").lp.mean()
lp2019 = PQ[(PQ.q >= pd.Period("2019Q1")) & (PQ.q <= pd.Period("2019Q4"))].groupby("geo").lp.mean()
dP = (lp2022 - lp2019).rename("dP")
R["energy_price_increase_2019_2022_logpts"] = {k: round(float(100 * v), 1) for k, v in dP.sort_values().items() if k in set(D.geo)}
XS = {}
for v in ["flat", "trend"]:
    c = CU[v].join(D.drop_duplicates("csid").set_index("csid")[["x"]]).join(dP, on="geo")
    c["energy"] = c.x * c.dP
    c["short_s"] = c.short / c.short.std(); c["energy_s"] = c.energy / c.energy.std()
    XS[v] = {}
    for w in ["2022-2026H1", "2023", "2024", "2025", "2026H1"]:
        c["y"] = c[f"ex_{w}"]
        m = pf.feols("y ~ energy | geo + nace_r2", data=c.dropna(subset=["y", "energy"]), vcov={"CRV1": "geo"}); t = m.tidy()
        XS[v][w] = {k: {"b": round(float(t.loc[k, "Estimate"]), 4), "se": round(float(t.loc[k, "Std. Error"]), 4), "p": round(float(t.loc[k, "Pr(>|t|)"]), 4)} for k in ["energy"]}
        XS[v][w]["nobs"] = int(m._N); XS[v][w]["r2_within"] = round(float(m._r2_within), 3)
    XS[v]["energy_sd"] = round(float(c.energy.std()), 4)
    if v == "flat": CS_FLAT = c.copy()
R["cross_section"] = XS
# variance decomposition of excess bankruptcies (flat counterfactual): country, sector and residual shares
c = CS_FLAT.dropna(subset=["ex_2022-2026H1"]); yv = c["ex_2022-2026H1"]
def r2(groups):
    Xd = pd.get_dummies(c[groups], drop_first=False, dtype=float).values; bb, *_ = np.linalg.lstsq(Xd, yv.values, rcond=None)
    return round(float(1 - ((yv.values - Xd @ bb) ** 2).sum() / ((yv.values - yv.mean()) ** 2).sum()), 3)
R["excess_variance_explained"] = {"country": r2(["geo"]), "sector": r2(["nace_r2"]), "country_and_sector": r2(["geo", "nace_r2"])}
# all-sector bankruptcies by country: 2025 relative to 2017-2019 (log points), all countries with data
allb = pd.read_csv("data/sts_rb_q_sca.csv"); allb = allb[(allb.indic_bt == "BKRT") & (allb.nace_r2 == "B-S_X_O_S94")].copy(); allb["y"] = allb.time.str[:4]
pre_ = allb[allb.y.isin(["2017", "2018", "2019"])].groupby("geo").value.mean(); y25 = allb[allb.y == "2025"].groupby("geo").value.mean()
R["country_2025_vs_2017_2019_logpts"] = {k: round(float(v), 0) for k, v in (100 * np.log(y25 / pre_)).dropna().sort_values().items()}
# deviation paths by sector (median across countries)
R["dev_path_flat"] = {SLAB[s]: {str(q): round(float(100 * v), 1) for q, v in g.groupby("q").dev_flat.median().items()} for s, g in D.groupby("nace_r2")}

# ---------------- descriptives ----------------
eu = S[S.geo == "EU27_2020"] if (S.geo == "EU27_2020").any() else None
raw = pd.read_csv("data/sts_rb_q_sca.csv"); raw = raw[(raw.geo == "EU27_2020") & (raw.indic_bt == "BKRT") & raw.nace_r2.isin(SECTORS) & (raw.time <= LAST)]
R["eu_bankruptcy_index"] = {SLAB[s]: {t: round(float(v), 1) for t, v in g.set_index("time").value.items() if t in ("2019-Q4", "2021-Q4", "2023-Q4", "2024-Q4", "2025-Q4", "2026-Q2")} for s, g in raw.groupby("nace_r2")}
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
def irf(ax, res, key, lab, color, off=0.0):
    b = np.array([res[h][key]["b"] for h in HZ]); se = np.array([res[h][key]["se"] for h in HZ])
    ax.plot(np.array(HZ) + off, 100 * b, "o-", color=color, ms=3, lw=1, label=lab)
    ax.fill_between(np.array(HZ) + off, 100 * (b - 1.96 * se), 100 * (b + 1.96 * se), color=color, alpha=0.15, lw=0)
fig, axs = plt.subplots(1, 3, figsize=(11, 3.4))
irf(axs[0], R["lp_bankruptcy"], "z", "Bankruptcies", "#d62728"); irf(axs[0], R["lp_registration"], "z", "Registrations", "#1f77b4")
axs[0].set_title("(a) Exposure × energy price shock"); axs[0].legend(frameon=False)
irf(axs[1], R["lp_asymmetry"], "zpos", "Price increases", "#d62728"); irf(axs[1], R["lp_asymmetry"], "zneg", "Price decreases", "#2ca02c")
axs[1].set_title("(b) Asymmetry (bankruptcies)"); axs[1].legend(frameon=False)
irf(axs[2], R["lp_state"], "z_tight", "Monetary tightening (2022Q3–2024Q2)", "#9467bd"); irf(axs[2], R["lp_state"], "z_loose", "Other quarters", "#7f7f7f")
axs[2].set_title("(c) State dependence (bankruptcies)"); axs[2].legend(frameon=False, fontsize=7)
for a in axs: a.axhline(0, color="grey", lw=0.6); a.set_xlabel("Quarters after the shock")
axs[0].set_ylabel("Response (log points × 100) to a 1-sd higher\nexposure and a 100% energy price increase")
fig.tight_layout(); fig.savefig("fig/fig5_irf.png"); plt.show()
# Fig 1: EU bankruptcy index by sector and energy PPI
raw2 = pd.read_csv("data/sts_rb_q_sca.csv"); raw2 = raw2[(raw2.geo == "EU27_2020") & (raw2.indic_bt == "BKRT") & raw2.nace_r2.isin(SECTORS) & (raw2.time <= LAST)]
pe = pd.read_csv("data/ppi_energy_m.csv"); pe = pe[pe.geo == "EU27_2020"]
fig, axs = plt.subplots(1, 2, figsize=(11, 3.6))
cols = plt.cm.tab10(np.linspace(0, 1, 10))
for (s, g), c_ in zip(raw2.groupby("nace_r2"), cols):
    axs[0].plot(range(len(g)), g.value.values, color=c_, lw=1.1 if s in ("B-E", "I", "H") else 0.8, label=SLAB[s])
tl = raw2.drop_duplicates("time").time.tolist(); axs[0].set_xticks(range(0, len(tl), 8)); axs[0].set_xticklabels([tl[i][:4] for i in range(0, len(tl), 8)])
axs[0].set_title("(a) EU bankruptcy declarations by sector (2021 = 100)"); axs[0].legend(frameon=False, fontsize=6)
if len(pe):
    pe = pe[pe.time >= "2015-01"]; axs[1].plot(pd.to_datetime(pe.time), pe.value, color="k", lw=1)
    axs[1].set_title("(b) EU producer price index, energy (domestic market, 2021 = 100)")
fig.tight_layout(); fig.savefig("fig/fig1_descriptive.png"); plt.show()
# Fig 3: deviation of bankruptcies from the 2017-2019 level by sector (median country)
fig, ax = plt.subplots(figsize=(9, 4.6))
for (s, g), c_ in zip(D.groupby("nace_r2"), plt.cm.tab10(np.linspace(0, 1, 10))):
    p_ = g.groupby("q").dev_flat.median(); ax.plot(range(len(p_)), 100 * p_.values, color=c_, lw=1.4 if s in ("B-E", "I", "H") else 0.9, label=SLAB[s])
ql = [str(q) for q in p_.index]; ax.set_xticks(range(0, len(ql), 4)); ax.set_xticklabels([q[:4] for q in ql[::4]])
ax.axhline(0, color="grey", lw=0.6); ax.axvspan(ql.index("2020Q2") - 0.5, ql.index("2021Q4") + 0.5, color="#dddddd", alpha=0.6, lw=0)
ax.set_ylabel("Bankruptcies relative to 2017–2019\n(log points × 100, median country)"); ax.legend(frameon=False, fontsize=7, ncol=4, loc="upper center", bbox_to_anchor=(0.5, -0.1))
fig.tight_layout(); fig.savefig("fig/fig3_deviation.png"); plt.show()

print(json.dumps({k: R[k] for k in R if not k.startswith("lp_")}, indent=1, default=str)[:6000])
for k in [k for k in R if k.startswith("lp_")]:
    print(k, {h: R[k][h] for h in (0, 4, 8)})

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

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