# -*- coding: utf-8 -*-
"""Hydropower droughts and fossil backfill: A Bayesian hierarchical and instrumental-variable analysis of 57 power systems, 2015–2026

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

# # Hydropower droughts and fossil backfill: A Bayesian hierarchical and instrumental-variable analysis of 57 power systems, 2015–2026
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür (replikasyon paketindeki betiklerle aynı kod, hücrelere bölünmüş hâlde).
#
# **Çalıştırma:** Çalışma zamanı → Tümünü çalıştır (*Runtime → Run all*). Tahmini süre: 5–10 dakika.
#
# **Veri:** Varsayılan olarak veriler özgün kaynakların açık API'lerinden canlı indirilir (`VERI_KAYNAGI = "api"`: Ember aylık elektrik verisi, WRI santral veritabanı ve NOAA CPC PREC/L yağış ızgarası (77 MB)). 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 pymc arviz xarray netCDF4 statsmodels

VERI_KAYNAGI = "api"   # "api": download from the original sources; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp14", exist_ok=True); os.chdir("/content/rp14")
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 hydropower-drought paper into data/: Ember monthly electricity data (generation by source, demand,
    net imports, emissions; to October 2026), the WRI Global Power Plant Database (hydropower plant locations and capacities) and
    NOAA CPC PREC/L monthly gauge-based land precipitation on a 1-degree grid (1948 - October 2026)."""
    from common import get, DATA
    if __name__ == "__main__":
        (DATA / "ember_monthly.csv").write_bytes(get("https://files.ember-energy.org/public-downloads/generation/outputs/release_generation_monthly_global.csv").content)
        (DATA / "gppd.csv").write_bytes(get("https://raw.githubusercontent.com/wri/global-power-plant-database/master/output_database/global_power_plant_database.csv").content)
        (DATA / "precl_1x1.nc").write_bytes(get("https://downloads.psl.noaa.gov/Datasets/precl/1.0deg/precip.mon.mean.1x1.nc").content)
        import pandas as pd
        pd.read_csv(DATA / "ember_monthly.csv").to_csv(DATA / "ember_monthly.csv.gz", index=False, compression="gzip"); (DATA / "ember_monthly.csv").unlink()
        print("done")
print(sorted(os.listdir('data')))

# ## 2. Standartlaştırılmış yağış indeksleri, Bayesçi hiyerarşik model, araç değişkenli ikame marjları, kuraklığa atfedilen fosil üretim ve CO2, tahmin, şekiller

# Paper 14: hydropower droughts and fossil backfill. Monthly electricity data for 60 countries with a hydropower share of at least 5% (Ember, to October 2026), hydropower-capacity-weighted precipitation (NOAA CPC PREC/L gauge analysis, WRI plant locations) and standardised precipitation indices; Bayesian hierarchical response of hydropower to precipitation, instrumental-variable estimates of the substitution margins, drought-attributable fossil generation and CO2, and the out-of-sample predictability of hydropower anomalies.
import json, sys, warnings
import numpy as np, pandas as pd
import xarray as xr
from scipy import stats
import matplotlib
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore")
plt.rcParams.update({"font.family": "DejaVu Serif", "font.size": 9, "axes.titlesize": 10, "legend.fontsize": 8,
                     "savefig.dpi": 300, "axes.spines.top": False, "axes.spines.right": False})
QUICK = "--quick" in sys.argv
R = {}

# ---------------- electricity data ----------------
E = pd.read_csv("data/ember_monthly.csv.gz")
E = E[E["Area type"] == "Country or economy"]
src = {"Hydro": "hydro", "Fossil": "fossil", "Coal": "coal", "Gas": "gas", "Wind and solar": "vre", "Total generation": "gen", "Demand": "demand", "Net imports": "netimp",
       "Nuclear": "nuclear", "Bioenergy": "bio", "Other renewables": "othren"}
X = E[E["Electricity source"].isin(src)].pivot_table(index=["ISO 3 code", "Area", "Date"], columns="Electricity source", values="Generation (TWh)").rename(columns=src).reset_index()
em = E[E["Electricity source"] == "Fossil"].set_index(["ISO 3 code", "Date"])["Emissions (MtCO2e)"].rename("fossil_co2")
X = X.join(em, on=["ISO 3 code", "Date"]).rename(columns={"ISO 3 code": "iso", "Area": "country", "Date": "date"})
X["date"] = pd.to_datetime(X.date)
X = X[(X.date >= "2015-01-01")]
share = X.groupby("iso").apply(lambda g: g.hydro.sum() / g.gen.sum())
keep = share[(share >= 0.05)].index.difference(["LUX"])                         # Luxembourg's hydro is pumped storage
n = X[X.iso.isin(keep)].groupby("iso").hydro.count()
keep = n[n >= 60].index
X = X[X.iso.isin(keep)].sort_values(["iso", "date"]).reset_index(drop=True)
for c in ["fossil", "netimp", "vre", "nuclear", "bio", "othren", "coal", "gas"]: X[c] = X[c].fillna(0)
R["countries"] = int(X.iso.nunique()); R["obs"] = int(len(X)); R["last_month"] = str(X.date.max().date())
R["hydro_share"] = {k: round(float(v), 3) for k, v in share[keep].sort_values(ascending=False).items()}

# ---------------- hydropower-weighted precipitation ----------------
g = pd.read_csv("data/gppd.csv", low_memory=False)
h = g[(g.primary_fuel == "Hydro") & g.country.isin(keep)].copy()
h["la"] = np.floor(h.latitude) + 0.5; h["lo"] = (np.floor(h.longitude) + 0.5) % 360
cells = h.groupby(["country", "la", "lo"]).capacity_mw.sum().reset_index().sort_values(["country", "capacity_mw"], ascending=[True, False])
cells["cum"] = cells.groupby("country").capacity_mw.cumsum() / cells.groupby("country").capacity_mw.transform("sum")
cells = cells[(cells.cum - cells.capacity_mw / cells.groupby("country").capacity_mw.transform("sum")) < 0.9]
cells = cells.groupby("country").head(40)
import os
if os.path.exists("data/precl_1x1.nc"):
    ds = xr.open_dataset("data/precl_1x1.nc")
    pr = ds.precip.sel(time=slice("1981-01-01", None))
    series = {}
    for c, gc_ in cells.groupby("country"):
        acc = None; wsum = 0.0
        for _, r in gc_.iterrows():
            s = pr.sel(lat=r.la, lon=r.lo, method="nearest").to_series()
            if s.isna().mean() > 0.2:                                   # coastal cell without land data: nearest valid neighbour
                box = pr.sel(lat=slice(r.la + 1.6, r.la - 1.6), lon=slice(r.lo - 1.6, r.lo + 1.6)).mean(["lat", "lon"]).to_series()
                s = box
            if s.isna().mean() > 0.2: continue
            acc = s * r.capacity_mw if acc is None else acc + s * r.capacity_mw; wsum += r.capacity_mw
        if acc is not None: series[c] = acc / wsum
    PR = pd.DataFrame(series)
    PR.index = pd.to_datetime(PR.index)
    PR.to_csv("out/precip_hydro_weighted.csv")
    PR = pd.read_csv("out/precip_hydro_weighted.csv", index_col=0, parse_dates=True)   # identical values in both paths
else:                                                       # replication package: capacity-weighted series derived from PREC/L
    PR = pd.read_csv("data/precip_hydro_weighted.csv", index_col=0, parse_dates=True)
R["precip_cells"] = {c: int(len(gc_)) for c, gc_ in cells.groupby("country")}
clim = PR[(PR.index.year >= 1991) & (PR.index.year <= 2020)]
mx = clim.groupby(clim.index.month).max(); mn = clim.groupby(clim.index.month).min()
rec = PR[PR.index >= "2025-01-01"]
outside = ((rec.values > mx.loc[rec.index.month].values * 1.5) | (rec.values < 0)).mean()
R["precip_screen"] = {"share_recent_months_above_1_5x_record": round(float(outside), 4), "last_month": str(PR.dropna(how="all").index.max().date())}
# standardised precipitation indices (empirical CDF against 1991-2020 by calendar month, Gringorten plotting positions)
def spi(P_, k):
    A = P_.rolling(k).mean(); out = pd.DataFrame(index=A.index, columns=A.columns, dtype=float)
    for m in range(1, 13):
        sel = A.index.month == m; base = A[(A.index.month == m) & (A.index.year >= 1991) & (A.index.year <= 2020)]
        for c in A.columns:
            b = np.sort(base[c].dropna().values); nb = len(b)
            if nb < 20: continue
            r_ = np.searchsorted(b, A.loc[sel, c].values, side="right")
            p = (r_ + 0.5 - 0.44) / (nb + 0.12)
            out.loc[sel, c] = stats.norm.ppf(np.clip(p, 0.5 / (nb + 1), 1 - 0.5 / (nb + 1)))
    return out
SPI = {k: spi(PR, k) for k in (1, 3, 6, 12)}
for k in SPI:
    s = SPI[k].stack().rename(f"spi{k}"); s.index.names = ["date", "iso"]; s = s.reset_index()
    X = X.merge(s, on=["date", "iso"], how="left")
X = X.dropna(subset=["spi3", "spi12"]).reset_index(drop=True)
R["countries_with_precip"] = int(X.iso.nunique())
X.to_csv("out/panel.csv", index=False)

# ---------------- residualise: country x calendar-month means and country-specific linear trends ----------------
X["moy"] = X.date.dt.month; X["t"] = (X.date.dt.year - 2015) * 12 + X.date.dt.month
X["lhydro"] = np.log(X.hydro.clip(lower=1e-4))
def resid(col):
    out = pd.Series(index=X.index, dtype=float)
    for c, g_ in X.groupby("iso"):
        D = pd.get_dummies(g_.moy, prefix="m", dtype=float); D["t"] = g_.t.values / 100
        b, *_ = np.linalg.lstsq(D.values, g_[col].values, rcond=None)
        out.loc[g_.index] = g_[col].values - D.values @ b
    return out
for col in ["lhydro", "spi1", "spi3", "spi6", "spi12", "hydro", "fossil", "netimp", "demand", "vre", "nuclear", "bio", "othren", "coal", "gas"]:
    X[f"r_{col}"] = resid(col)

# ---------------- A. pooled and Bayesian hierarchical response of hydropower to precipitation ----------------
import statsmodels.api as sm
def cluster_ols(y, Xm, groups):
    m = sm.OLS(y, Xm).fit(cov_type="cluster", cov_kwds={"groups": groups}); return m
pool = cluster_ols(X.r_lhydro, X[["r_spi3", "r_spi12"]], X.iso.astype("category").cat.codes)
R["pooled"] = {k: {"b": round(float(pool.params[k]), 4), "se": round(float(pool.bse[k]), 4)} for k in pool.params.index}
pool_multi = cluster_ols(X.r_lhydro, X[["r_spi1", "r_spi3", "r_spi6", "r_spi12"]], X.iso.astype("category").cat.codes)
R["pooled_multiscale"] = {k: {"b": round(float(pool_multi.params[k]), 4), "se": round(float(pool_multi.bse[k]), 4)} for k in pool_multi.params.index}
import pymc as pm
cidx, clist = pd.factorize(X.iso)
with pm.Model() as model:
    mu = pm.Normal("mu", 0, 0.5, shape=2); tau = pm.HalfNormal("tau", 0.3, shape=2)
    z = pm.Normal("z", 0, 1, shape=(len(clist), 2)); beta = pm.Deterministic("beta", mu + z * tau)
    sig = pm.HalfNormal("sigma", 0.5, shape=len(clist)); nu = pm.Exponential("nu_m1", 1 / 10) + 1
    mean = beta[cidx, 0] * X.r_spi3.values + beta[cidx, 1] * X.r_spi12.values
    pm.StudentT("y", nu=nu, mu=mean, sigma=sig[cidx], observed=X.r_lhydro.values)
    idata = pm.sample(500 if QUICK else 1500, tune=500 if QUICK else 1500, chains=4, cores=2, target_accept=0.9, random_seed=14, progressbar=False)
post = idata.posterior
B = post["beta"].stack(s=("chain", "draw")).values                     # (country, 2, samples)
import arviz as az
R["hier"] = {"mu": [round(float(v), 4) for v in post["mu"].mean(("chain", "draw")).values],
             "mu_sd": [round(float(v), 4) for v in post["mu"].std(("chain", "draw")).values],
             "tau": [round(float(v), 4) for v in post["tau"].mean(("chain", "draw")).values],
             "rhat_max": round(float(az.summary(idata, var_names=["mu", "tau", "beta"]).r_hat.max()), 3),
             "ess_bulk_min": int(az.summary(idata, var_names=["mu", "tau", "beta"]).ess_bulk.min())}
eff = pd.DataFrame({"iso": clist, "b3": B[:, 0, :].mean(1), "b12": B[:, 1, :].mean(1),
                    "b3_lo": np.quantile(B[:, 0, :], 0.05, axis=1), "b3_hi": np.quantile(B[:, 0, :], 0.95, axis=1),
                    "b12_lo": np.quantile(B[:, 1, :], 0.05, axis=1), "b12_hi": np.quantile(B[:, 1, :], 0.95, axis=1)}).set_index("iso")
# elasticity of a one-sd drought on both scales (SPI3 = SPI12 = -1), in % of hydropower generation
dr = -(B[:, 0, :] + B[:, 1, :])
eff["drought_effect_pct"] = 100 * (np.exp(-(B[:, 0, :] + B[:, 1, :])) - 1).mean(1)
eff["drought_lo"] = 100 * (np.exp(np.quantile(-(B[:, 0, :] + B[:, 1, :]), 0.05, axis=1)) - 1)
eff["drought_hi"] = 100 * (np.exp(np.quantile(-(B[:, 0, :] + B[:, 1, :]), 0.95, axis=1)) - 1)
eff["country"] = X.drop_duplicates("iso").set_index("iso").country
eff["hydro_share"] = share
eff.to_csv("out/country_effects.csv")
R["drought_effect_pct_quantiles"] = {k: round(float(v), 1) for k, v in eff.drought_effect_pct.quantile([0.1, 0.5, 0.9]).items()}
R["drought_effect_selected"] = {c: round(float(eff.loc[c, "drought_effect_pct"]), 1) for c in ["TUR", "NOR", "BRA", "CHN", "ITA", "ESP", "CHE", "CAN", "COL", "VNM"] if c in eff.index}

# ---------------- B. substitution margins: instrumental variables ----------------
# All flows are residualised (country x calendar month, country trend) and divided by the country's mean monthly generation, so that
# coefficients are pooled across countries of different size. The instrument is the precipitation-predicted hydropower deviation
# (country posterior slopes x residualised SPI-3/SPI-12 x the country's mean hydropower in that calendar month), a function of
# precipitation only; a second specification uses SPI-3 and SPI-12 interacted with the hydropower share.
X = X.join(eff[["b3", "b12"]], on="iso")
X["dlog"] = X.b3 * X.r_spi3 + X.b12 * X.r_spi12                               # precipitation-driven log deviation of hydropower
gbar = X.groupby("iso").gen.transform("mean"); mh = X.groupby(["iso", "moy"]).hydro.transform("mean")
X["zpred"] = X.dlog * mh / gbar
X["hs"] = X.iso.map(share)
for col in ["hydro", "fossil", "netimp", "demand", "vre", "nuclear", "bio", "othren", "coal", "gas"]:
    X[f"n_{col}"] = X[f"r_{col}"] / gbar
def iv(y, x, Z, groups):
    """2SLS of y on x (one endogenous regressor), instruments Z; cluster-robust by country; first-stage cluster-robust F."""
    Zm = np.column_stack(Z); pi, *_ = np.linalg.lstsq(Zm, x, rcond=None); xh = Zm @ pi
    b = (xh @ y) / (xh @ x); u = y - b * x
    meat = sum(((xh[groups == g_] * u[groups == g_]).sum()) ** 2 for g_ in np.unique(groups))
    se = np.sqrt(meat / (xh @ xh) ** 2)
    fs = sm.OLS(x, Zm).fit(cov_type="cluster", cov_kwds={"groups": groups})
    F = float(fs.wald_test(np.eye(Zm.shape[1]), scalar=True).statistic) / Zm.shape[1]
    return round(float(b), 4), round(float(se), 4), round(F, 1)
grp = X.iso.astype("category").cat.codes.values
LABS = [("fossil", "Fossil generation"), ("coal", "Coal"), ("gas", "Gas"), ("netimp", "Net imports"), ("demand", "Demand"),
        ("vre", "Wind and solar"), ("nuclear", "Nuclear"), ("bio", "Bioenergy"), ("othren", "Other renewables")]
margins, margins2 = {}, {}
for col, lab in LABS:
    b, se, F = iv(X[f"n_{col}"].values, X.n_hydro.values, [X.zpred.values], grp); margins[lab] = {"b": b, "se": se, "first_stage_F": F}
    b, se, F = iv(X[f"n_{col}"].values, X.n_hydro.values, [X.r_spi3.values * X.hs.values, X.r_spi12.values * X.hs.values], grp); margins2[lab] = {"b": b, "se": se, "first_stage_F": F}
R["margins_iv"] = margins; R["margins_iv_spi_x_share"] = margins2
ols = cluster_ols(X.n_fossil, X[["n_hydro"]], grp); R["margins_ols_fossil"] = {"b": round(float(ols.params.n_hydro), 4), "se": round(float(ols.bse.n_hydro), 4)}
# by region
EUROPE = {"NOR", "AUT", "CHE", "MNE", "LVA", "SWE", "HRV", "BIH", "SVN", "SRB", "ROU", "MKD", "FIN", "PRT", "ITA", "SVK", "ESP", "FRA", "LTU", "BGR", "GRC", "ISL", "MDA", "TUR", "GEO", "ARM", "AZE"}
AMER = {"COL", "ECU", "CRI", "BRA", "CAN", "PER", "URY", "ARG", "BOL", "CHL", "SLV", "MEX", "USA", "DOM"}
X["region"] = np.where(X.iso.isin(EUROPE), "Europe", np.where(X.iso.isin(AMER), "Americas", "Asia, Africa and Oceania"))
reg = {}
for rg, g_ in X.groupby("region"):
    gg = g_.iso.astype("category").cat.codes.values
    reg[rg] = {lab: dict(zip(["b", "se", "F"], iv(g_[f"n_{col}"].values, g_.n_hydro.values, [g_.zpred.values], gg)))
               for col, lab in [("fossil", "Fossil generation"), ("coal", "Coal"), ("gas", "Gas"), ("netimp", "Net imports"), ("demand", "Demand")]}
    reg[rg]["countries"] = int(g_.iso.nunique())
R["margins_by_region"] = reg

# ---------------- C. drought-attributable hydropower shortfall, fossil generation and CO2 ----------------
bF = X.region.map({rg: reg[rg]["Fossil generation"]["b"] for rg in reg})      # region-specific dFossil / dHydro (negative)
R["attribution_fossil_margin_by_region"] = {rg: reg[rg]["Fossil generation"]["b"] for rg in reg}
X["hydro_cf"] = X.hydro * np.exp(-X.dlog)                                      # hydropower with normal precipitation
X["dhydro"] = X.hydro - X.hydro_cf                                             # TWh, negative in dry months
X["dfossil"] = bF * X.dhydro                                                   # TWh
X["ef"] = (X.fossil_co2 / X.fossil.replace(0, np.nan)).fillna(X.groupby("iso").fossil_co2.transform("sum") / X.groupby("iso").fossil.transform("sum"))
X["dco2"] = X.dfossil * X.ef                                                   # MtCO2
X["year"] = X.date.dt.year
bal = X[X.date >= "2019-01-01"]
yr = bal.groupby("year")[["dhydro", "dfossil", "dco2"]].sum()
R["attributable_by_year"] = {int(y): {"hydro_twh": round(float(v.dhydro), 1), "fossil_twh": round(float(v.dfossil), 1), "co2_mt": round(float(v.dco2), 1)} for y, v in yr.iterrows()}
R["attributable_countries_in_balanced_window"] = int(bal.iso.nunique())
cy = X.groupby(["iso", "year"])[["dhydro", "dco2"]].sum()
top = cy.dco2.sort_values(ascending=False).head(12)
R["attributable_top_country_years"] = [[i, int(y), round(float(v), 1), round(float(cy.loc[(i, y), "dhydro"]), 1)] for (i, y), v in top.items()]
R["turkiye"] = {int(y): {"hydro_twh": round(float(cy.loc[("TUR", y), "dhydro"]), 2), "co2_mt": round(float(cy.loc[("TUR", y), "dco2"]), 2)} for y in sorted(X[X.iso == "TUR"].year.unique())} if "TUR" in X.iso.values else None
X.to_csv("out/panel.csv", index=False)

# ---------------- D. predictability: out-of-sample forecasts of hydropower anomalies from lagged precipitation ----------------
fc = {}
for lead in [1, 2, 3, 6]:
    rows = []
    for c, g_ in X.groupby("iso"):
        g_ = g_.sort_values("date").copy()
        g_["x1"] = g_.r_spi3.shift(lead); g_["x2"] = g_.r_spi12.shift(lead); g_["ylag"] = g_.r_lhydro.shift(lead)
        rows.append(g_[["iso", "date", "r_lhydro", "x1", "x2", "ylag"]].dropna())
    Fd = pd.concat(rows)
    errs = {"climatology": [], "persistence": [], "precipitation": [], "persistence + precipitation": []}
    for t in pd.date_range("2021-01-01", X.date.max(), freq="MS"):
        tr = Fd[Fd.date < t - pd.DateOffset(months=lead - 1)]; te = Fd[Fd.date == t]
        if len(te) == 0 or len(tr) < 500: continue
        for name, cols in [("persistence", ["ylag"]), ("precipitation", ["x1", "x2"]), ("persistence + precipitation", ["ylag", "x1", "x2"])]:
            b, *_ = np.linalg.lstsq(tr[cols].values, tr.r_lhydro.values, rcond=None)
            errs[name] += list((te.r_lhydro.values - te[cols].values @ b) ** 2)
        errs["climatology"] += list(te.r_lhydro.values ** 2)
    base_ = np.mean(errs["climatology"])
    fc[lead] = {k: round(float(1 - np.mean(v) / base_), 3) for k, v in errs.items() if k != "climatology"}
R["oos_r2_vs_climatology"] = fc
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
# Fig 1: country drought effects (posterior mean and 90% interval)
ef = eff.sort_values("drought_effect_pct")
fig, ax = plt.subplots(figsize=(6.6, 8.2))
y_ = np.arange(len(ef))
ax.errorbar(ef.drought_effect_pct, y_, xerr=[ef.drought_effect_pct - ef.drought_lo, ef.drought_hi - ef.drought_effect_pct], fmt="o", ms=3, color="k", lw=0.8)
ax.axvline(0, color="grey", lw=0.6); ax.axvline(100 * (np.exp(-sum(R["hier"]["mu"])) - 1), color="#d62728", ls="--", lw=0.8)
ax.set_yticks(y_); ax.set_yticklabels(ef.country, fontsize=6.5); ax.set_xlabel("Change in hydropower generation for SPI-3 = SPI-12 = −1 (%)")
fig.tight_layout(); fig.savefig("fig/fig1_country_effects.png"); plt.show()
# Fig 2: examples
ex = [c for c in ["TUR", "BRA", "CHN", "ITA", "NOR", "COL"] if c in X.iso.values]
fig, axs = plt.subplots(3, 2, figsize=(9, 6.6), sharex=True)
for a_, c in zip(axs.flat, ex):
    g_ = X[X.iso == c]
    a_.plot(g_.date, 100 * g_.r_lhydro, color="#1f77b4", lw=1, label="Hydropower anomaly (%)")
    a2 = a_.twinx(); a2.plot(g_.date, g_.spi12, color="#d62728", lw=0.9, label="SPI-12"); a2.axhline(0, color="grey", lw=0.5); a2.set_ylim(-3.5, 3.5)
    a_.set_title(g_.country.iloc[0]); a_.axhline(0, color="grey", lw=0.5)
axs[0, 0].legend(frameon=False, loc="lower left", fontsize=7)
fig.autofmt_xdate(); fig.tight_layout(); fig.savefig("fig/fig2_examples.png"); plt.show()
# Fig 3: substitution margins (IV)
fig, ax = plt.subplots(figsize=(6.4, 3.4))
labs = ["Fossil generation", "Net imports", "Demand", "Wind and solar", "Nuclear", "Bioenergy", "Other renewables"]
b = np.array([-margins[l]["b"] for l in labs[:2]] + [margins["Demand"]["b"]] + [-margins[l]["b"] for l in labs[3:]]); se = np.array([margins[l]["se"] for l in labs])
ax.barh(range(len(labs)), 100 * b, xerr=196 * se, color=["#7f7f7f", "#1f77b4", "#2ca02c", "#ff7f0e", "#9467bd", "#8c564b", "#e377c2"], capsize=2)
ax.set_yticks(range(len(labs))); ax.set_yticklabels(["More fossil generation", "More net imports", "Lower demand", "More wind and solar", "More nuclear", "More bioenergy", "More other renewables"])
ax.invert_yaxis(); ax.axvline(0, color="grey", lw=0.6); ax.set_xlabel("Share of a 1 TWh hydropower shortfall (%)")
fig.tight_layout(); fig.savefig("fig/fig3_margins.png"); plt.show()
# Fig 4: attributable CO2 by year and region
ry = bal.groupby(["year", "region"]).dco2.sum().unstack().fillna(0)
fig, ax = plt.subplots(figsize=(6.6, 3.4))
pos = np.zeros(len(ry)); neg = np.zeros(len(ry))
for rg, c_ in zip(ry.columns, ["#1f77b4", "#2ca02c", "#d62728"]):
    v = ry[rg].values; ax.bar(ry.index, v, bottom=np.where(v >= 0, pos, neg), color=c_, label=rg)
    pos = pos + np.where(v >= 0, v, 0); neg = neg + np.where(v < 0, v, 0)
ax.plot(ry.index, ry.sum(1), "ko-", ms=3, lw=0.8, label="Net")
ax.axhline(0, color="grey", lw=0.6); ax.set_ylabel("Precipitation-driven change in\npower-sector CO$_2$ (Mt)")
lo_, hi_ = ax.get_ylim(); ax.set_ylim(lo_, hi_ + 0.35 * (hi_ - lo_)); ax.legend(frameon=False, fontsize=7, ncol=4, loc="upper center")
fig.tight_layout(); fig.savefig("fig/fig4_attributable.png"); plt.show()
print(json.dumps({k: R[k] for k in R if k not in ("hydro_share",)}, indent=1, default=str)[:8000])

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

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