# -*- coding: utf-8 -*-
"""Local, global and hierarchical forecasting of monthly maritime CO₂ emissions for 140 economies

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

# # Local, global and hierarchical forecasting of monthly maritime CO₂ emissions for 140 economies
#
# Bu not defteri makaledeki **tüm sayıları, tabloları ve şekilleri** üreten kodun Colab sürümüdür (makalenin 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: 25–40 dakika (çapraz doğrulama ve LightGBM).
#
# **Veri:** Varsayılan olarak veriler GOMD açık veri API'sinden canlı çekilir (`VERI_KAYNAGI = "api"`). Makaledeki sayıları birebir yeniden üretmek için replikasyon paketi ZIP'ini yükleyebilirsiniz (`VERI_KAYNAGI = "zip"`): API'deki seriler kaynak kurumlar tarafından güncellenirse sonuçlar küçük farklar gösterebilir. Veri çekim tarihi makalede 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 and the GOMD Python client. If Colab asks to restart the session, restart and run all cells again.
# [Colab Magic] !pip install -q statsforecast==2.1.1 mlforecast==1.1.0 hierarchicalforecast==1.5.3 lightgbm pycountry pycountry-convert
# [Colab Magic] !pip install -q --upgrade gomd

VERI_KAYNAGI = "api"   # "api": live from GOMD; "zip": upload the replication package ZIP (exact reproduction)
import os
os.makedirs("/content/rp01", exist_ok=True); os.chdir("/content/rp01")
for d in ("data", "out", "fig"): os.makedirs(d, exist_ok=True)
if VERI_KAYNAGI == "zip":
    import zipfile, glob, shutil
    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/"): z.extract(n, ".")
    print(sorted(os.listdir("data")))
print("Calisma klasoru:", os.getcwd())

# ## 1. Veri

# Data from the public GOMD open data API (entity lists as used in the paper). Requests run in a few parallel threads with retries.
import json, time
from concurrent.futures import ThreadPoolExecutor
import pandas as pd
from gomd import GOMDClient
client = GOMDClient("https://gomd.pirireis.edu.tr")
JOBS = [["OECD_MARITIME_CO2_RESIDENCE_M", "oecd_maritime_co2_monthly"]]
ENTITIES = json.loads("{\"OECD_MARITIME_CO2_RESIDENCE_M\": [[\"TUR\", \"Sovereign Coastal State TUR\", \"country\"], [\"ASM\", \"American Samoa\", \"country\"], [\"AIA\", \"Anguilla\", \"country\"], [\"ATG\", \"Antigua and Barbuda\", \"country\"], [\"ABW\", \"Aruba\", \"country\"], [\"AZE\", \"Azerbaijan\", \"country\"], [\"BHS\", \"Bahamas, The\", \"country\"], [\"GGY\", \"Bailiwick of Guernsey\", \"country\"], [\"BRB\", \"Barbados\", \"country\"], [\"BLZ\", \"Belize\", \"country\"], [\"BMU\", \"Bermuda\", \"country\"], [\"VGB\", \"British Virgin Islands\", \"country\"], [\"BRN\", \"Brunei Darussalam\", \"country\"], [\"CYM\", \"Cayman Islands\", \"country\"], [\"COK\", \"Cook Islands\", \"country\"], [\"CUW\", \"Curacao\", \"country\"], [\"GNQ\", \"Equatorial Guinea\", \"country\"], [\"FLK\", \"Falkland Islands (Malvinas)\", \"country\"], [\"FRO\", \"Faroe Islands\", \"country\"], [\"PYF\", \"French Polynesia\", \"country\"], [\"GIB\", \"Gibraltar\", \"country\"], [\"GRL\", \"Greenland\", \"country\"], [\"GLP\", \"Guadeloupe\", \"country\"], [\"GTM\", \"Guatemala\", \"country\"], [\"GUY\", \"Guyana\", \"country\"], [\"HND\", \"Honduras\", \"country\"], [\"HKG\", \"Hong Kong SAR, China\", \"country\"], [\"IMN\", \"Isle of Man\", \"country\"], [\"JEY\", \"Jersey\", \"country\"], [\"KAZ\", \"Kazakhstan\", \"country\"], [\"KIR\", \"Kiribati\", \"country\"], [\"PRK\", \"Korea, Dem. People's Rep.\", \"country\"], [\"LUX\", \"Luxembourg\", \"country\"], [\"MAC\", \"Macao SAR, China\", \"country\"], [\"MHL\", \"Marshall Islands\", \"country\"], [\"FSM\", \"Micronesia, Fed. Sts.\", \"country\"], [\"MCO\", \"Monaco\", \"country\"], [\"MMR\", \"Myanmar\", \"country\"], [\"NRU\", \"Naoero\", \"country\"], [\"NCL\", \"New Caledonia\", \"country\"], [\"OECD\", \"Organisation for Economic Co-operation and Development (OECD)\", \"aggregate\"], [\"W_O\", \"Other countries\", \"aggregate\"], [\"PLW\", \"Palau\", \"country\"], [\"PRY\", \"Paraguay\", \"country\"], [\"PRI\", \"Puerto Rico (US)\", \"country\"], [\"WSM\", \"Samoa\", \"country\"], [\"SRB\", \"Serbia\", \"country\"], [\"SLB\", \"Solomon Islands\", \"country\"], [\"AGO\", \"Sovereign Coastal State AGO\", \"country\"], [\"ALB\", \"Sovereign Coastal State ALB\", \"country\"], [\"ARE\", \"Sovereign Coastal State ARE\", \"country\"], [\"ARG\", \"Sovereign Coastal State ARG\", \"country\"], [\"AUS\", \"Sovereign Coastal State AUS\", \"country\"], [\"BEL\", \"Sovereign Coastal State BEL\", \"country\"], [\"BGD\", \"Sovereign Coastal State BGD\", \"country\"], [\"BGR\", \"Sovereign Coastal State BGR\", \"country\"], [\"BHR\", \"Sovereign Coastal State BHR\", \"country\"], [\"BRA\", \"Sovereign Coastal State BRA\", \"country\"], [\"CAN\", \"Sovereign Coastal State CAN\", \"country\"], [\"CHL\", \"Sovereign Coastal State CHL\", \"country\"], [\"CHN\", \"Sovereign Coastal State CHN\", \"country\"], [\"CIV\", \"Sovereign Coastal State CIV\", \"country\"], [\"CMR\", \"Sovereign Coastal State CMR\", \"country\"], [\"COG\", \"Sovereign Coastal State COG\", \"country\"], [\"COL\", \"Sovereign Coastal State COL\", \"country\"], [\"COM\", \"Sovereign Coastal State COM\", \"country\"], [\"CPV\", \"Sovereign Coastal State CPV\", \"country\"], [\"CUB\", \"Sovereign Coastal State CUB\", \"country\"], [\"CYP\", \"Greek Cypriot Administration of Southern Cyprus (GCA)\", \"country\"], [\"DEU\", \"Sovereign Coastal State DEU\", \"country\"], [\"DNK\", \"Sovereign Coastal State DNK\", \"country\"], [\"DOM\", \"Sovereign Coastal State DOM\", \"country\"], [\"DZA\", \"Sovereign Coastal State DZA\", \"country\"], [\"ECU\", \"Sovereign Coastal State ECU\", \"country\"], [\"EGY\", \"Sovereign Coastal State EGY\", \"country\"], [\"ESP\", \"Sovereign Coastal State ESP\", \"country\"], [\"EST\", \"Sovereign Coastal State EST\", \"country\"], [\"FIN\", \"Sovereign Coastal State FIN\", \"country\"], [\"FJI\", \"Sovereign Coastal State FJI\", \"country\"], [\"FRA\", \"Sovereign Coastal State FRA\", \"country\"], [\"GAB\", \"Sovereign Coastal State GAB\", \"country\"], [\"GBR\", \"Sovereign Coastal State GBR\", \"country\"], [\"GEO\", \"Sovereign Coastal State GEO\", \"country\"], [\"GHA\", \"Sovereign Coastal State GHA\", \"country\"], [\"GIN\", \"Sovereign Coastal State GIN\", \"country\"], [\"GRC\", \"Sovereign Coastal State GRC\", \"country\"], [\"HRV\", \"Sovereign Coastal State HRV\", \"country\"], [\"IDN\", \"Sovereign Coastal State IDN\", \"country\"], [\"IND\", \"Sovereign Coastal State IND\", \"country\"], [\"IRL\", \"Sovereign Coastal State IRL\", \"country\"], [\"IRN\", \"Sovereign Coastal State IRN\", \"country\"], [\"IRQ\", \"Sovereign Coastal State IRQ\", \"country\"], [\"ISL\", \"Sovereign Coastal State ISL\", \"country\"], [\"ITA\", \"Sovereign Coastal State ITA\", \"country\"], [\"JAM\", \"Sovereign Coastal State JAM\", \"country\"], [\"JOR\", \"Sovereign Coastal State JOR\", \"country\"], [\"JPN\", \"Sovereign Coastal State JPN\", \"country\"], [\"KEN\", \"Sovereign Coastal State KEN\", \"country\"], [\"KHM\", \"Sovereign Coastal State KHM\", \"country\"], [\"KOR\", \"Sovereign Coastal State KOR\", \"country\"], [\"KWT\", \"Sovereign Coastal State KWT\", \"country\"], [\"LBN\", \"Sovereign Coastal State LBN\", \"country\"], [\"LBR\", \"Sovereign Coastal State LBR\", \"country\"], [\"LKA\", \"Sovereign Coastal State LKA\", \"country\"], [\"LTU\", \"Sovereign Coastal State LTU\", \"country\"], [\"LVA\", \"Sovereign Coastal State LVA\", \"country\"], [\"MAR\", \"Sovereign Coastal State MAR\", \"country\"], [\"MDG\", \"Sovereign Coastal State MDG\", \"country\"], [\"MDV\", \"Sovereign Coastal State MDV\", \"country\"], [\"MEX\", \"Sovereign Coastal State MEX\", \"country\"], [\"MLT\", \"Sovereign Coastal State MLT\", \"country\"], [\"MNE\", \"Sovereign Coastal State MNE\", \"country\"], [\"MOZ\", \"Sovereign Coastal State MOZ\", \"country\"], [\"MRT\", \"Sovereign Coastal State MRT\", \"country\"], [\"MUS\", \"Sovereign Coastal State MUS\", \"country\"], [\"MYS\", \"Sovereign Coastal State MYS\", \"country\"], [\"NAM\", \"Sovereign Coastal State NAM\", \"country\"], [\"NGA\", \"Sovereign Coastal State NGA\", \"country\"], [\"NLD\", \"Sovereign Coastal State NLD\", \"country\"], [\"NOR\", \"Sovereign Coastal State NOR\", \"country\"], [\"NZL\", \"Sovereign Coastal State NZL\", \"country\"], [\"OMN\", \"Sovereign Coastal State OMN\", \"country\"], [\"PAK\", \"Sovereign Coastal State PAK\", \"country\"], [\"PAN\", \"Sovereign Coastal State PAN\", \"country\"], [\"PER\", \"Sovereign Coastal State PER\", \"country\"], [\"PHL\", \"Sovereign Coastal State PHL\", \"country\"], [\"PNG\", \"Sovereign Coastal State PNG\", \"country\"], [\"POL\", \"Sovereign Coastal State POL\", \"country\"], [\"PRT\", \"Sovereign Coastal State PRT\", \"country\"], [\"QAT\", \"Sovereign Coastal State QAT\", \"country\"], [\"ROU\", \"Sovereign Coastal State ROU\", \"country\"], [\"RUS\", \"Sovereign Coastal State RUS\", \"country\"], [\"SAU\", \"Sovereign Coastal State SAU\", \"country\"], [\"SEN\", \"Sovereign Coastal State SEN\", \"country\"], [\"SGP\", \"Sovereign Coastal State SGP\", \"country\"], [\"SLE\", \"Sovereign Coastal State SLE\", \"country\"], [\"SVN\", \"Sovereign Coastal State SVN\", \"country\"], [\"SWE\", \"Sovereign Coastal State SWE\", \"country\"], [\"SYC\", \"Sovereign Coastal State SYC\", \"country\"], [\"SYR\", \"Sovereign Coastal State SYR\", \"country\"], [\"THA\", \"Sovereign Coastal State THA\", \"country\"], [\"TTO\", \"Sovereign Coastal State TTO\", \"country\"], [\"TUN\", \"Sovereign Coastal State TUN\", \"country\"], [\"TWN\", \"Sovereign Coastal State TWN\", \"country\"], [\"TZA\", \"Sovereign Coastal State TZA\", \"country\"], [\"UKR\", \"Sovereign Coastal State UKR\", \"country\"], [\"URY\", \"Sovereign Coastal State URY\", \"country\"], [\"USA\", \"Sovereign Coastal State USA\", \"country\"], [\"VEN\", \"Sovereign Coastal State VEN\", \"country\"], [\"VNM\", \"Sovereign Coastal State VNM\", \"country\"], [\"WLD\", \"Sovereign Coastal State WLD\", \"country\"], [\"YEM\", \"Sovereign Coastal State YEM\", \"country\"], [\"ZAF\", \"Sovereign Coastal State ZAF\", \"country\"], [\"KNA\", \"St. Kitts and Nevis\", \"country\"], [\"VCT\", \"St. Vincent and the Grenadines\", \"country\"], [\"SUR\", \"Suriname\", \"country\"], [\"CHE\", \"Switzerland\", \"country\"], [\"TLS\", \"Timor-Leste\", \"country\"], [\"TON\", \"Tonga\", \"country\"], [\"TKM\", \"Turkmenistan\", \"country\"], [\"TCA\", \"Turks and Caicos Islands\", \"country\"], [\"VUT\", \"Vanuatu\", \"country\"], [\"VIR\", \"Virgin Islands (U.S.)\", \"country\"]]}")
CAP = {"OECD_MARITIME_CO2_RESIDENCE_M": "2026-06"}   # last period used in the paper (later months, if published, are excluded)

def fetch(indicator, entity, retries=6):
    for k in range(retries):
        try:
            rows = client.get_series(indicator=indicator, entity=entity, limit=50000, as_df=False)
            return [(str(r["period"]), float(r["value"])) for r in rows if str(r.get("entity_code") or "") == entity and r.get("value") is not None]
        except Exception as e:
            if k == retries - 1: raise
            time.sleep(2 ** k)

def fetch_panel(indicator, out_csv):
    ents = ENTITIES[indicator]
    with ThreadPoolExecutor(max_workers=4) as ex:
        res = list(ex.map(lambda e: fetch(indicator, e[0]), ents))
    recs = [(c, n, t, p, v) for (c, n, t), rows in zip(ents, res) for p, v in rows]
    df = pd.DataFrame(recs, columns=["entity", "name", "etype", "period", "value"]).drop_duplicates(["entity", "period"])
    if indicator in CAP: df = df[df.period <= CAP[indicator]]
    df.to_csv(out_csv, index=False)
    return df

if VERI_KAYNAGI == "api":
    for ind, f in JOBS:
        t = time.time(); df = fetch_panel(ind, f"data/{f}.csv")
        print(f"{f:28s} {len(df):7d} rows  {df.entity.nunique():4d} entities  {df.period.min()}-{df.period.max()}  ({time.time() - t:.0f} s)")

# ## 2. Hiyerarşinin kurulması (140 ekonomi + kalan)

# Build the hierarchy World > Region > Economy from the OECD residence-principle maritime CO2 panel.
import pandas as pd, numpy as np, pycountry, pycountry_convert as pcc
df = pd.read_csv("data/oecd_maritime_co2_monthly.csv", dtype={"period": str})
piv = df.pivot(index="period", columns="entity", values="value").sort_index()
econ = [c for c in piv.columns if c not in ("WLD", "W_O", "OECD")]
full = [c for c in econ if piv[c].notna().all()]
MAN = {"TWN": "AS", "HKG": "AS", "MAC": "AS", "XKX": "EU", "GGY": "EU", "JEY": "EU", "IMN": "EU", "FLK": "SA", "ANT": "NA", "TLS": "AS", "SXM": "NA", "CUW": "NA", "BES": "NA", "VAT": "EU"}
REG = {"AF": "Africa", "AS": "Asia", "EU": "Europe", "NA": "Americas", "SA": "Americas", "OC": "Oceania"}
def region(iso3):
    if iso3 in MAN: return REG[MAN[iso3]]
    try:
        a2 = pycountry.countries.get(alpha_3=iso3).alpha_2
        return REG[pcc.country_alpha2_to_continent_code(a2)]
    except Exception:
        return None
reg = {c: region(c) for c in full}
miss = [c for c, r in reg.items() if r is None]
print("unmapped:", miss)
rows = []
for c in full:
    for p, v in piv[c].items():
        rows.append(("World", reg[c], c, p, v))
rest = piv["WLD"] - piv[full].sum(axis=1)
for p, v in rest.items():
    rows.append(("World", "Rest", "REST", p, v))
bot = pd.DataFrame(rows, columns=["World", "Region", "Economy", "period", "y"])
bot["ds"] = pd.to_datetime(bot.period + "-01")
bot.to_csv("data/bottom.csv", index=False)
print("economies with full 90 months:", len(full), "| rest share of world:", round((rest / piv.WLD).mean(), 4), "min rest", rest.min())
print(bot.groupby("Region").Economy.nunique())
names = df[["entity", "name"]].drop_duplicates().set_index("entity").name
names.to_csv("data/names.csv")

# ## 3. Yerel ve küresel modellerle kayan başlangıçlı tahminler

# Stage 1: rolling-origin forecasts from local statistical and global machine-learning models. 13 forecast origins (2024-06 ... 2025-06), horizons 1-12, every node of the World > Region > Economy hierarchy.
import os, warnings, time
import numpy as np, pandas as pd
warnings.filterwarnings("ignore")
os.environ["NIXTLA_ID_AS_COL"] = "1"
from hierarchicalforecast.utils import aggregate
from statsforecast import StatsForecast
from statsforecast.models import Naive, SeasonalNaive, AutoETS, AutoARIMA, AutoTheta
from mlforecast import MLForecast
from mlforecast.target_transforms import LocalStandardScaler, GlobalSklearnTransformer
from mlforecast.lag_transforms import RollingMean
from mlforecast.utils import PredictionIntervals
from sklearn.preprocessing import FunctionTransformer
from sklearn.linear_model import Ridge
import lightgbm as lgb

H, W, LEVEL = 12, 13, [80, 95]
bot = pd.read_csv("data/bottom.csv", parse_dates=["ds"])
spec = [["World"], ["World", "Region"], ["World", "Region", "Economy"]]
Y, S, tags = aggregate(bot[["World", "Region", "Economy", "ds", "y"]], spec)
Y = Y.sort_values(["unique_id", "ds"]).reset_index(drop=True)
Y.to_csv("out/Y.csv", index=False)
S_df = S if isinstance(S, pd.DataFrame) else pd.DataFrame(S)
S_df.to_csv("out/S.csv", index=False)
pd.Series({k: list(v) for k, v in tags.items()}).to_json("out/tags.json")
print("nodes:", Y.unique_id.nunique(), {k: len(v) for k, v in tags.items()})

if not os.path.exists("out/cv_local.csv"):
    t = time.time()
    sf = StatsForecast(models=[Naive(), SeasonalNaive(season_length=12), AutoETS(season_length=12),
                               AutoARIMA(season_length=12), AutoTheta(season_length=12)], freq="MS", n_jobs=-1)
    cv = sf.cross_validation(df=Y, h=H, n_windows=W, step_size=1, level=LEVEL, fitted=True)
    fit = sf.cross_validation_fitted_values()
    cv.to_csv("out/cv_local.csv", index=False); fit.to_csv("out/cv_local_fitted.csv", index=False)
    print("local done", round(time.time() - t, 1), "s", cv.shape)

if not os.path.exists("out/cv_global.csv"):
    t = time.time()
    log_tf = GlobalSklearnTransformer(FunctionTransformer(func=np.log, inverse_func=np.exp))
    mlf = MLForecast(
        models={"LightGBM": lgb.LGBMRegressor(n_estimators=400, learning_rate=0.03, num_leaves=15, min_child_samples=20,
                                              subsample=0.8, subsample_freq=1, colsample_bytree=0.8, random_state=1, verbosity=-1),
                "GlobalAR": Ridge(alpha=1.0)},
        freq="MS", lags=[1, 2, 3, 6, 12],
        lag_transforms={1: [RollingMean(window_size=3), RollingMean(window_size=12)]},
        date_features=["month"], target_transforms=[log_tf, LocalStandardScaler()], num_threads=2)
    cvg = mlf.cross_validation(df=Y, n_windows=W, h=H, step_size=1, refit=True,
                               prediction_intervals=PredictionIntervals(n_windows=2, h=H), level=LEVEL)
    cvg.to_csv("out/cv_global.csv", index=False)
    print("global done", round(time.time() - t, 1), "s", cvg.shape)

# ## 4. Mevsimsel farklı LightGBM

# Stage 1b: global LightGBM on seasonally differenced log series (robustness variant).
import os, warnings, numpy as np, pandas as pd
warnings.filterwarnings("ignore"); os.environ["NIXTLA_ID_AS_COL"] = "1"
from mlforecast import MLForecast
from mlforecast.target_transforms import LocalStandardScaler, GlobalSklearnTransformer, Differences
from mlforecast.lag_transforms import RollingMean
from mlforecast.utils import PredictionIntervals
from sklearn.preprocessing import FunctionTransformer
import lightgbm as lgb
Y = pd.read_csv("out/Y.csv", parse_dates=["ds"])
log_tf = GlobalSklearnTransformer(FunctionTransformer(func=np.log, inverse_func=np.exp))
mlf = MLForecast(models={"LightGBM_SD": lgb.LGBMRegressor(n_estimators=400, learning_rate=0.03, num_leaves=15, min_child_samples=20,
                    subsample=0.8, subsample_freq=1, colsample_bytree=0.8, random_state=1, verbosity=-1)},
                 freq="MS", lags=[1, 2, 3, 6, 12], lag_transforms={1: [RollingMean(window_size=3), RollingMean(window_size=12)]},
                 date_features=["month"], target_transforms=[log_tf, Differences([12]), LocalStandardScaler()], num_threads=2)
cv = mlf.cross_validation(df=Y, n_windows=13, h=12, step_size=1, refit=True, prediction_intervals=PredictionIntervals(n_windows=2, h=12), level=[80, 95])
cv.drop(columns=["y"]).to_csv("out/cv_global_b.csv", index=False); print(cv.shape)

# ## 5. Doğruluk, testler, uzlaştırma ve şekiller

# Stage 2: accuracy, calibration, significance tests, reconciliation and series-feature analysis; figures and results.json.
import os, json, warnings
import numpy as np, pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
from statsmodels.tsa.seasonal import STL
warnings.filterwarnings("ignore")
os.environ["NIXTLA_ID_AS_COL"] = "1"
from hierarchicalforecast.core import HierarchicalReconciliation
from hierarchicalforecast.methods import BottomUp, TopDown, MinTrace

plt.rcParams.update({"font.family": "DejaVu Serif", "font.size": 9, "axes.titlesize": 10, "axes.labelsize": 9,
                     "legend.fontsize": 8, "figure.dpi": 100, "savefig.dpi": 300, "axes.spines.top": False, "axes.spines.right": False})
R = {}
Y = pd.read_csv("out/Y.csv", parse_dates=["ds"])
S_df = pd.read_csv("out/S.csv")
tags = {k: np.array(v) for k, v in json.load(open("out/tags.json")).items()}
lvl = {u: "World" for u in tags["World"]}
lvl.update({u: "Region" for u in tags["World/Region"]})
lvl.update({u: "Economy" for u in tags["World/Region/Economy"]})
loc = pd.read_csv("out/cv_local.csv", parse_dates=["ds", "cutoff"])
glo = pd.read_csv("out/cv_global.csv", parse_dates=["ds", "cutoff"])
cv = loc.merge(glo.drop(columns=["y"]), on=["unique_id", "ds", "cutoff"], how="inner").merge(pd.read_csv("out/cv_global_b.csv", parse_dates=["ds", "cutoff"]), on=["unique_id", "ds", "cutoff"], how="inner")
NAMES = {"Naive": "Naive", "SeasonalNaive": "sNaive", "AutoETS": "ETS", "AutoARIMA": "ARIMA", "AutoTheta": "Theta",
         "Comb": "Comb (ETS+ARIMA+Theta)", "LightGBM": "LightGBM (global)", "LightGBM_SD": "LightGBM-SD (global)", "GlobalAR": "Ridge-AR (global)", "Hybrid": "Hybrid (Comb+LightGBM)"}
for suf in ["", "-lo-80", "-hi-80", "-lo-95", "-hi-95"]:
    cv["Comb" + suf] = cv[["AutoETS" + suf, "AutoARIMA" + suf, "AutoTheta" + suf]].mean(axis=1)
    cv["Hybrid" + suf] = cv[["Comb" + suf, "LightGBM" + suf]].mean(axis=1)
MODELS = list(NAMES)
cv["h"] = ((cv.ds.dt.year - cv.cutoff.dt.year) * 12 + (cv.ds.dt.month - cv.cutoff.dt.month)).astype(int)
cv["level"] = cv.unique_id.map(lvl)

# ---- scales from training data (seasonal naive in-sample MAE / MSE) ----
sc = []
for c in sorted(cv.cutoff.unique()):
    tr = Y[Y.ds <= c].copy()
    tr["lag12"] = tr.groupby("unique_id").y.shift(12)
    d = (tr.y - tr.lag12)
    g = pd.DataFrame({"unique_id": tr.unique_id, "ae": d.abs(), "se": d ** 2}).dropna().groupby("unique_id").mean()
    g["cutoff"] = c
    sc.append(g.reset_index())
sc = pd.concat(sc)
cv = cv.merge(sc, on=["unique_id", "cutoff"])

def interval_score(y, lo, hi, a):
    return (hi - lo) + (2 / a) * np.clip(lo - y, 0, None) + (2 / a) * np.clip(y - hi, 0, None)

rows = []
for m in MODELS:
    e = cv[["unique_id", "cutoff", "h", "level", "ae", "se", "y"]].copy()
    e["model"] = m
    e["ase"] = (cv[m] - cv.y).abs() / cv.ae
    e["sse"] = (cv[m] - cv.y) ** 2 / cv.se
    e["sis95"] = interval_score(cv.y, cv[m + "-lo-95"], cv[m + "-hi-95"], 0.05) / cv.ae
    e["cov80"] = ((cv.y >= cv[m + "-lo-80"]) & (cv.y <= cv[m + "-hi-80"])).astype(float)
    e["cov95"] = ((cv.y >= cv[m + "-lo-95"]) & (cv.y <= cv[m + "-hi-95"])).astype(float)
    e["ape"] = (cv[m] - cv.y).abs() / cv.y * 100
    rows.append(e)
E = pd.concat(rows)
summ = E.groupby(["level", "model"]).agg(MASE=("ase", "mean"), RMSSE=("sse", lambda s: np.sqrt(s.mean())), MSIS=("sis95", "mean"),
                                          cov80=("cov80", "mean"), cov95=("cov95", "mean"), MAPE=("ape", "mean")).reset_index()
R["n_nodes"] = {k: int(len(v)) for k, v in tags.items()}
R["n_origins"] = int(cv.cutoff.nunique()); R["origins"] = [str(cv.cutoff.min().date()), str(cv.cutoff.max().date())]
R["n_forecasts_per_model"] = int(len(cv))
R["summary"] = {f"{r.level}|{r.model}": {k: round(float(getattr(r, k)), 4) for k in ["MASE", "RMSSE", "MSIS", "cov80", "cov95", "MAPE"]} for r in summ.itertuples()}
byh = E[E.level == "Economy"].groupby(["model", "h"]).agg(MASE=("ase", "mean"), MSIS=("sis95", "mean")).reset_index()
R["mase_by_h_economy"] = {m: [round(float(v), 4) for v in byh[byh.model == m].sort_values("h").MASE] for m in MODELS}

# ---- ranks, Friedman and Nemenyi (economy level; series as blocks) ----
ser = E[E.level == "Economy"].groupby(["unique_id", "model"]).ase.mean().unstack()[MODELS]
ranks = ser.rank(axis=1)
k, N = len(MODELS), len(ser)
q = {2: 1.960, 3: 2.343, 4: 2.569, 5: 2.728, 6: 2.850, 7: 2.949, 8: 3.031, 9: 3.102, 10: 3.164}[k]  # Demsar (2006), alpha = 0.05
CD = q * np.sqrt(k * (k + 1) / (6 * N))
fr = stats.friedmanchisquare(*[ser[m] for m in MODELS])
R["ranks"] = {m: round(float(ranks[m].mean()), 3) for m in MODELS}
R["nemenyi_CD"] = round(float(CD), 3); R["friedman"] = {"chi2": round(float(fr.statistic), 2), "p": float(fr.pvalue), "N": N, "k": k}
best = min(R["ranks"], key=R["ranks"].get)
R["best_rank_model"] = best
R["not_sig_diff_from_best"] = [m for m in MODELS if R["ranks"][m] - R["ranks"][best] <= CD]
# pairwise Wilcoxon: LightGBM vs Comb, LightGBM vs ETS (economy level, series means)
for a, b in [("LightGBM", "Comb"), ("LightGBM", "AutoETS"), ("Comb", "AutoETS"), ("GlobalAR", "Comb"), ("LightGBM_SD", "LightGBM"), ("LightGBM_SD", "Comb"), ("Hybrid", "Comb"), ("Hybrid", "AutoETS")]:
    w = stats.wilcoxon(ser[a], ser[b])
    R[f"wilcoxon_{a}_vs_{b}"] = {"median_ratio": round(float((ser[a] / ser[b]).median()), 4), "share_a_better": round(float((ser[a] < ser[b]).mean()), 4), "p": float(w.pvalue)}

# ---- reconciliation ----
fit = pd.read_csv("out/cv_local_fitted.csv", parse_dates=["ds", "cutoff"])
rec_rows = []
for c in sorted(cv.cutoff.unique()):
    base = cv[cv.cutoff == c][["unique_id", "ds", "AutoETS", "LightGBM"]]
    ins = fit[fit.cutoff == c][["unique_id", "ds", "y", "AutoETS"]].copy()
    ins["LightGBM"] = ins["AutoETS"]  # placeholder: methods used for LightGBM below do not use residuals
    hr = HierarchicalReconciliation(reconcilers=[BottomUp(), TopDown(method="forecast_proportions"), MinTrace(method="ols"),
                                                 MinTrace(method="wls_struct"), MinTrace(method="mint_shrink")])
    out = hr.reconcile(Y_hat_df=base, Y_df=ins, S_df=S_df, tags=tags)
    out["cutoff"] = c
    rec_rows.append(out)
REC = pd.concat(rec_rows).merge(cv[["unique_id", "ds", "cutoff", "y", "ae", "level", "h"]], on=["unique_id", "ds", "cutoff"])
keep = [c for c in REC.columns if c.startswith("AutoETS") or (c.startswith("LightGBM") and "mint_shrink" not in c and "TopDown" not in c)]
rec_sum = {}
for col in keep:
    a = ((REC[col] - REC.y).abs() / REC.ae).groupby(REC.level).mean()
    rec_sum[col] = {k2: round(float(v), 4) for k2, v in a.items()}
R["reconciliation_MASE"] = rec_sum
R["reconciliation_rel"] = {col: {lv: round(rec_sum[col][lv] / rec_sum[col.split("/")[0]][lv], 4) for lv in rec_sum[col]} for col in keep if "/" in col}
REC.to_csv("out/reconciled.csv", index=False)

# ---- series features and when global helps ----
feat = []
cut0 = cv.cutoff.min()
for u in tags["World/Region/Economy"]:
    s = Y[(Y.unique_id == u) & (Y.ds <= cut0)].set_index("ds").y
    if len(s) < 36 or (s <= 0).any():
        continue
    r = STL(np.log(s.values), period=12, robust=True).fit()
    vr = np.var(r.resid)
    FT = max(0.0, 1 - vr / np.var(r.trend + r.resid)); FS = max(0.0, 1 - vr / np.var(r.seasonal + r.resid))
    feat.append((u, np.log10(s.mean()), FT, FS, s.std() / s.mean()))
F = pd.DataFrame(feat, columns=["unique_id", "log10_level", "F_trend", "F_season", "cv"]).set_index("unique_id")
F["rel_lgbm_comb"] = np.log(ser["LightGBM"] / ser["Comb"]).reindex(F.index)
F = F.dropna()
X = sm.add_constant(F[["log10_level", "F_trend", "F_season"]])
ols = sm.OLS(F.rel_lgbm_comb, X).fit(cov_type="HC1")
R["feature_regression"] = {"n": int(ols.nobs), "r2": round(float(ols.rsquared), 3),
                           "coef": {k2: round(float(v), 4) for k2, v in ols.params.items()},
                           "p": {k2: float(v) for k2, v in ols.pvalues.items()}}
R["features_median"] = {c: round(float(F[c].median()), 3) for c in ["log10_level", "F_trend", "F_season", "cv"]}
F.to_csv("out/features.csv")

# ---- descriptive data facts ----
w = Y[Y.unique_id == "World"].set_index("ds").y
R["world"] = {"first": [str(w.index[0].date()), round(float(w.iloc[0] / 1e6), 2)], "last": [str(w.index[-1].date()), round(float(w.iloc[-1] / 1e6), 2)],
              "annual_2019": round(float(w["2019"].sum() / 1e6), 1), "annual_2025": round(float(w["2025"].sum() / 1e6), 1),
              "min_month": [str(w.idxmin().date()), round(float(w.min() / 1e6), 2)], "max_month": [str(w.idxmax().date()), round(float(w.max() / 1e6), 2)]}
econ = Y[Y.unique_id.isin(tags["World/Region/Economy"])]
tot = econ.groupby("unique_id").y.sum().sort_values(ascending=False)
R["top10_share"] = round(float(tot.head(10).sum() / w.sum()), 4)
R["top10"] = [(u.split("/")[-1], round(float(tot[u] / w.sum() * 100), 2)) for u in tot.head(10).index]
R["region_of"] = {u.split("/")[-1]: u.split("/")[1] for u in tags["World/Region/Economy"]}
R["region_counts"] = {u.split("/")[-1]: int(sum(1 for x in tags["World/Region/Economy"] if x.split("/")[1] == u.split("/")[-1])) for u in tags["World/Region"]}
R["tur"] = {"share_pct": round(float(Y[Y.unique_id.str.endswith("/TUR")].y.sum() / w.sum() * 100), 3), "mase": {m: round(float(E[(E.unique_id.str.endswith("/TUR"))&(E.model==m)].ase.mean()), 3) for m in MODELS}}
R["region_share_pct"] = {u.split("/")[-1]: round(float(Y[Y.unique_id == u].y.sum() / w.sum() * 100), 2) for u in tags["World/Region"]}
R["level_range_economy_Mt_month"] = [round(float(econ.groupby("unique_id").y.mean().min() / 1e6), 4), round(float(econ.groupby("unique_id").y.mean().max() / 1e6), 2)]
json.dump(R, open("out/results.json", "w"), indent=1, default=str)

# ================= FIGURES =================
import matplotlib.dates as mdates
UID = {u.split("/")[-1]: u for u in Y.unique_id.unique()}
def yfmt(a):
    a.xaxis.set_major_locator(mdates.YearLocator(2)); a.xaxis.set_major_formatter(mdates.DateFormatter("%Y"))
COL = {"Naive": "#9e9e9e", "SeasonalNaive": "#616161", "AutoETS": "#1f77b4", "AutoARIMA": "#2ca02c", "AutoTheta": "#9467bd",
       "Comb": "#000000", "LightGBM": "#d62728", "LightGBM_SD": "#e377c2", "GlobalAR": "#ff7f0e", "Hybrid": "#17becf"}
# Fig 1 data
fig, ax = plt.subplots(1, 3, figsize=(10, 3.1))
ax[0].plot(w.index, w / 1e6, color="k", lw=1); ax[0].plot(w.index, (w / 1e6).rolling(12).mean(), color="#d62728", lw=1.2)
ax[0].set_title("(a) World, monthly maritime CO$_2$"); ax[0].set_ylabel("Mt CO$_2$ per month"); yfmt(ax[0])
for code, c in [("GRC", "#1f77b4"), ("CHN", "#d62728"), ("SGP", "#2ca02c"), ("TUR", "#ff7f0e"), ("USA", "#9467bd")]:
    u = UID[code]; s = Y[Y.unique_id == u].set_index("ds").y
    ax[1].plot(s.index, s / s["2019"].mean() * 100, lw=1, color=c, label=u.split("/")[-1])
ax[1].axhline(100, color="grey", lw=0.6, ls=":"); ax[1].set_title("(b) Selected economies (2019 = 100)"); ax[1].legend(ncol=2, frameon=False); yfmt(ax[1])
ax[2].scatter(F.F_trend, F.F_season, s=10 + 6 * (F.log10_level - F.log10_level.min()), c="#1f77b4", alpha=0.6, edgecolor="none")
ax[2].set_xlabel("Trend strength $F_T$"); ax[2].set_ylabel("Seasonal strength $F_S$"); ax[2].set_title("(c) Features of 141 bottom-level series")
fig.tight_layout(); fig.savefig("fig/fig1_data.png"); plt.show()
# Fig 2 accuracy by horizon
fig, ax = plt.subplots(1, 2, figsize=(9, 3.3))
for m in MODELS:
    d = byh[byh.model == m].sort_values("h")
    ax[0].plot(d.h, d.MASE, marker="o", ms=2.5, lw=1.2 if m in ("Comb", "LightGBM") else 0.8, color=COL[m], label=NAMES[m])
    ax[1].plot(d.h, d.MSIS, marker="o", ms=2.5, lw=1.2 if m in ("Comb", "LightGBM") else 0.8, color=COL[m])
ax[0].set_xlabel("Forecast horizon (months)"); ax[0].set_ylabel("MASE"); ax[0].set_title("(a) Point accuracy, bottom level (141 series)")
ax[1].set_xlabel("Forecast horizon (months)"); ax[1].set_ylabel("MSIS (95%)"); ax[1].set_title("(b) Interval accuracy")
ax[0].legend(frameon=False, ncol=2, fontsize=7); fig.tight_layout(); fig.savefig("fig/fig2_horizon.png"); plt.show()
# Fig 3 ranks with Nemenyi CD
order = sorted(MODELS, key=lambda m: R["ranks"][m])
fig, ax = plt.subplots(figsize=(6.5, 3.2))
for i, m in enumerate(order):
    r = R["ranks"][m]; ax.errorbar(r, i, xerr=CD / 2, fmt="o", color=COL[m], capsize=3)
ax.axvspan(R["ranks"][order[0]] - CD / 2, R["ranks"][order[0]] + CD / 2, color="#cfe2f3", alpha=0.6, lw=0)
ax.set_yticks(range(len(order))); ax.set_yticklabels([NAMES[m] for m in order]); ax.invert_yaxis()
ax.set_xlabel(f"Mean rank of MASE across 141 bottom-level series (Nemenyi CD = {CD:.2f})")
fig.tight_layout(); fig.savefig("fig/fig3_ranks.png"); plt.show()
# Fig 4 calibration
fig, ax = plt.subplots(figsize=(6.5, 3.2))
xs = np.arange(len(MODELS))
c80 = [summ[(summ.level == "Economy") & (summ.model == m)].cov80.iloc[0] * 100 for m in MODELS]
c95 = [summ[(summ.level == "Economy") & (summ.model == m)].cov95.iloc[0] * 100 for m in MODELS]
ax.bar(xs - 0.2, c80, 0.4, color="#9ecae1", label="80% interval"); ax.bar(xs + 0.2, c95, 0.4, color="#3182bd", label="95% interval")
ax.axhline(80, color="#9ecae1", ls="--", lw=1); ax.axhline(95, color="#3182bd", ls="--", lw=1)
ax.set_xticks(xs); ax.set_xticklabels([NAMES[m].split(" (")[0] for m in MODELS], rotation=30, ha="right")
ax.set_ylabel("Empirical coverage (%)"); ax.set_ylim(40, 105); ax.legend(frameon=False, loc="upper center", ncol=2, bbox_to_anchor=(0.5, 1.12))
fig.tight_layout(); fig.savefig("fig/fig4_coverage.png"); plt.show()
# Fig 5 reconciliation
meths = [c for c in keep if c.startswith("AutoETS/")]
fig, ax = plt.subplots(figsize=(7, 3.2))
levels = ["World", "Region", "Economy"]
for j, lv in enumerate(levels):
    vals = [R["reconciliation_rel"][mm][lv] for mm in meths]
    ax.bar(np.arange(len(meths)) + (j - 1) * 0.27, vals, 0.27, label=lv, color=["#08519c", "#6baed6", "#c6dbef"][j])
ax.axhline(1, color="k", lw=0.8)
ax.set_xticks(range(len(meths))); ax.set_xticklabels([mm.split("/")[1].replace("_method-", " ").replace("MinTrace", "MinT") for mm in meths], rotation=25, ha="right")
ax.set_ylabel("MASE relative to base ETS"); ax.legend(frameon=False, ncol=3); fig.tight_layout(); fig.savefig("fig/fig5_reconciliation.png"); plt.show()
# Fig 6 example forecasts at last origin
last = cv.cutoff.max()
fig, axs = plt.subplots(2, 2, figsize=(9, 5.2))
for a, u in zip(axs.flat, ["World", UID["TUR"], UID["GRC"], UID["SGP"]]):
    hist = Y[(Y.unique_id == u) & (Y.ds <= last + pd.DateOffset(months=12))].set_index("ds").y / 1e6
    f = cv[(cv.unique_id == u) & (cv.cutoff == last)].set_index("ds")
    a.plot(hist.index, hist, color="k", lw=0.9, label="Observed")
    for m, c in [("Comb", "#1f77b4"), ("LightGBM", "#d62728")]:
        a.plot(f.index, f[m] / 1e6, color=c, lw=1.2, label=NAMES[m].split(" (")[0])
        a.fill_between(f.index, f[m + "-lo-95"] / 1e6, f[m + "-hi-95"] / 1e6, color=c, alpha=0.15, lw=0)
    a.axvline(last, color="grey", ls=":", lw=0.8); a.set_title({"TUR": "Türkiye", "GRC": "Greece", "SGP": "Singapore"}.get(u.split("/")[-1], "World total")); a.set_ylabel("Mt CO$_2$/month"); yfmt(a)
axs[0, 0].legend(frameon=False, fontsize=7); fig.tight_layout(); fig.savefig("fig/fig6_examples.png"); plt.show()
# Fig 7 when does global help
fig, ax = plt.subplots(1, 2, figsize=(9, 3.2))
ax[0].scatter(F.log10_level, F.rel_lgbm_comb, s=12, alpha=0.6, color="#d62728", edgecolor="none"); ax[0].axhline(0, color="k", lw=0.7)
ax[0].set_xlabel("log$_{10}$ mean monthly emissions (t CO$_2$)"); ax[0].set_ylabel("log(MASE LightGBM / MASE Comb)")
ax[1].scatter(F.F_season, F.rel_lgbm_comb, s=12, alpha=0.6, color="#1f77b4", edgecolor="none"); ax[1].axhline(0, color="k", lw=0.7)
ax[1].set_xlabel("Seasonal strength $F_S$"); ax[1].set_ylabel("log(MASE LightGBM / MASE Comb)")
fig.tight_layout(); fig.savefig("fig/fig7_features.png"); plt.show()
print(json.dumps({k: R[k] for k in ["ranks", "nemenyi_CD", "friedman", "best_rank_model", "not_sig_diff_from_best", "feature_regression"]}, 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/rp01_outputs", "zip", ".", ".")
files.download("/content/rp01_outputs.zip")