# Supplementary Code for:
# Beyond participation: where carbon accounts cannot distinguish sustainable practice from development constraint
# Anonymous author(s)
# Reproduces every estimate, classification, table and figure in the article.
# Inputs: Global Carbon Budget 2025 (Global Carbon Project dataset v15, https://doi.org/10.18160/GCP-2025)
# via Our World in Data; World Bank World Development Indicators.
# All sources are publicly downloadable without registration.
# Requires: Python 3.10+, pandas, numpy, requests, statsmodels, scipy, matplotlib.

import io
import pandas as pd
import numpy as np
import requests
import statsmodels.api as sm
from scipy.stats import spearmanr, shapiro, skew
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import FancyBboxPatch, FancyArrowPatch

THRESH_T = 2.3          # 1.5 C-compatible operational benchmark (t CO2), frozen ex ante; sensitivity 2.1 / 2.5
POP_MIN = 300_000       # sample rule, frozen ex ante
ALPHA = 0.05            # two-tailed
WORLD_POP = 8.06e9

# Colourblind-safe palette (Okabe-Ito)
C_ZONE = "#7F7F7F"   # grey: blindness zone
C_ON   = "#0072B2"   # blue: on expectation
C_OVER = "#D55E00"   # vermillion: above expectation
C_TEAL = "#009E73"   # bluish green: realised-sustainability bracket

UA = {"User-Agent": "research-reproduction/1.0"}

# 1. Data: consumption-based CO2 per capita (Global Carbon Budget 2025 via Our World in Data)
r = requests.get("https://ourworldindata.org/grapher/consumption-co2-per-capita.csv", headers=UA, timeout=120)
r.raise_for_status()
owid = pd.read_csv(io.StringIO(r.text))
owid.columns = ["Entity", "Code", "Year", "co2_cons_pc"]
owid = owid[owid["Code"].str.len() == 3]
co2_now = owid[owid["Year"] == 2023].set_index("Code")[["Entity", "co2_cons_pc"]]

def wb(indicator, start=2020, end=2023):
    """World Development Indicators, most recent non-null observation per country in [start, end]."""
    url = (f"https://api.worldbank.org/v2/country/all/indicator/{indicator}"
           f"?format=json&date={start}:{end}&per_page=2000")
    d = pd.DataFrame(requests.get(url, headers=UA, timeout=120).json()[1])
    d = d[d["countryiso3code"].str.len() == 3]
    d["value"] = pd.to_numeric(d["value"]); d["date"] = d["date"].astype(int)
    return (d.dropna(subset=["value"]).sort_values("date")
             .groupby("countryiso3code").tail(1).set_index("countryiso3code")["value"])

df = co2_now.join([wb("NY.GDP.PCAP.PP.CD").rename("gdp_pc"), wb("EG.ELC.ACCS.ZS").rename("elec"),
                   wb("SP.URB.TOTL.IN.ZS").rename("urb"), wb("SP.POP.TOTL").rename("pop")],
                  how="inner").dropna()
df = df[df["pop"] > POP_MIN]
assert len(df) == 118

# 2. Constraint-consistent expectation model (frozen specification)
def design(d):
    X = pd.concat([np.log(d["gdp_pc"]).rename("ln_gdp"), d["elec"], d["urb"]], axis=1)
    X.insert(0, "const", 1.0)
    return X

def fit(d):
    return sm.OLS(np.log(d["co2_cons_pc"]), design(d)).fit()

def classify(model, d):
    pr = model.get_prediction(design(d)).summary_frame(ALPHA)
    exp_t = np.exp(pr["mean"])
    ln_y = np.log(d["co2_cons_pc"])
    cls = pd.Series(np.select([exp_t < THRESH_T,
                               (exp_t >= THRESH_T) & (ln_y < pr["obs_ci_lower"]),
                               (exp_t >= THRESH_T) & (ln_y > pr["obs_ci_upper"])],
                              ["Blindness zone", "Below expectation", "Above expectation"],
                              default="On expectation"), index=d.index)
    return cls, exp_t, pr

m = fit(df)
df["class"], df["expected_t"], pr = classify(m, df)
df["gap_pct"] = 100 * (df["co2_cons_pc"] / df["expected_t"] - 1)
df["total_mt"] = df["co2_cons_pc"] * df["pop"] / 1e6
print(m.summary())                                   # R2 = 0.861, n = 118
zone = df[df["class"] == "Blindness zone"]
print("Zone:", len(zone), "countries,", round(zone["pop"].sum() / 1e9, 2), "bn people")

# 2a. Normality check on model residuals (Nature statistical reporting requirement)
W, p_norm = shapiro(m.resid)
print("Shapiro-Wilk W = %.4f, P = %.4f; residual skew = %.3f" % (W, p_norm, skew(m.resid)))
# Expected: P = 0.0521, skew = 0.361

# 2b. Detection limit of the classification (how far below expectation registers as restraint)
ratio = np.exp(pr["obs_ci_lower"]) / np.exp(pr["mean"])
print("Lower 95%% prediction bound as share of expected: %.1f%% to %.1f%%" % (100*ratio.min(), 100*ratio.max()))
print("Implied shortfall needed to classify below expectation: %.0f%% to %.0f%%" % (100*(1-ratio.max()), 100*(1-ratio.min())))
# Expected: 36.8% to 39.6%; shortfall 60% to 63%

# 3. Persistence: identical model on 2000 and 2010 data
resid = {2023: pd.Series(m.resid, index=df.index)}
for yr in [2000, 2010]:
    c = owid[owid["Year"] == yr].set_index("Code")[["co2_cons_pc"]]
    d = c.join([wb("NY.GDP.PCAP.PP.CD", yr, yr + 1).rename("gdp_pc"),
                wb("EG.ELC.ACCS.ZS", yr, yr + 1).rename("elec"),
                wb("SP.URB.TOTL.IN.ZS", yr, yr + 1).rename("urb")], how="inner").dropna()
    resid[yr] = pd.Series(fit(d).resid, index=d.index)
P = pd.concat([resid[2000].rename("r2000"), resid[2010].rename("r2010"), resid[2023].rename("r2023")], axis=1).dropna()
print("Spearman 2000-2010: %.2f | 2010-2023: %.2f | 2000-2023: %.2f" % (
    spearmanr(P.r2000, P.r2010).statistic, spearmanr(P.r2010, P.r2023).statistic, spearmanr(P.r2000, P.r2023).statistic))
# Expected: 0.65 | 0.64 | 0.43 (n = 117)

# 4. Dual reporting
pc_rank = df["co2_cons_pc"].rank(ascending=False); tot_rank = df["total_mt"].rank(ascending=False)
print("Top-10 total emitters outside per-capita top 20:", int(((tot_rank <= 10) & (pc_rank > 20)).sum()))
print("India: pc rank %d (%.2f t) | total rank %d (%.2f Gt)" % (
    pc_rank["IND"], df.loc["IND", "co2_cons_pc"], tot_rank["IND"], df.loc["IND", "total_mt"] / 1000))
# Expected: 7; India 92nd (1.77 t) | 3rd (2.54 Gt)

# 5. Leave-one-country-out loop (Supplementary Table S1)
rows = []
for c in df.index:
    d_ = df.drop(index=c)
    ml = fit(d_); cls_, _, _ = classify(ml, d_)
    rows.append({"removed": df.loc[c, "Entity"],
                 "zone_bn": round(d_[cls_ == "Blindness zone"]["pop"].sum() / 1e9, 2),
                 "R2": round(ml.rsquared, 3),
                 "class_changes": int((cls_ != df["class"].drop(index=c)).sum())})
loo = pd.DataFrame(rows)
print(loo.to_string(index=False))
# Expected: 104 exclusions change 0 classifications, 14 change exactly 1; India exclusion: 0 changes, R2 0.860

# 6. Threshold sensitivity
for t in [2.1, 2.3, 2.5]:
    z = df[df["expected_t"] < t]["pop"].sum()
    print(f"Threshold {t} t: {z/1e9:.2f} bn = {100*z/WORLD_POP:.1f}% of world population")
# Expected: 36.7% / 38.0% / 39.2%

# 7. Post hoc out-of-sample check (added after pre-submission review; not part of the frozen design)
# Each country's 95% prediction interval is computed from a model fitted WITHOUT that country.
rows = []
for c in df.index:
    d_ = df.drop(index=c)
    ml = fit(d_)
    pr_c = ml.get_prediction(design(df.loc[[c]])).summary_frame(ALPHA)
    ln_y = np.log(df.loc[c, "co2_cons_pc"])
    rows.append({"Entity": df.loc[c, "Entity"],
                 "heldout_expected_t": round(float(np.exp(pr_c["mean"].iloc[0])), 3),
                 "heldout_lo_t": round(float(np.exp(pr_c["obs_ci_lower"].iloc[0])), 3),
                 "heldout_hi_t": round(float(np.exp(pr_c["obs_ci_upper"].iloc[0])), 3),
                 "realised_t": df.loc[c, "co2_cons_pc"],
                 "oos_below": bool(ln_y < pr_c["obs_ci_lower"].iloc[0]),
                 "oos_above": bool(ln_y > pr_c["obs_ci_upper"].iloc[0])})
oos = pd.DataFrame(rows)
print("Out-of-sample below band:", int(oos["oos_below"].sum()), "| above band:", int(oos["oos_above"].sum()))
# Expected: 0 below; 4 above (Kuwait, Mongolia, Namibia, Trinidad and Tobago);
# blindness-zone classification identical under held-out expectations (35 countries).

# 8. Tables
top15 = df.nlargest(15, "total_mt").sort_values("Entity")
table1 = top15[["Entity", "gdp_pc", "expected_t", "co2_cons_pc", "gap_pct", "total_mt", "class"]]
print("\nTable 1 (15 largest population-scale consumers, alphabetical):")
print(table1.to_string(index=False))
df.sort_values("Entity")[["Entity", "gdp_pc", "elec", "urb", "pop", "co2_cons_pc",
                          "expected_t", "gap_pct", "total_mt", "class"]].to_csv(
    "Supplementary_Data_S1_analysis_dataset.csv", index=False)
loo.to_csv("Table_S1_loo_loop.csv", index=False)
oos.to_csv("Table_S3_out_of_sample.csv", index=False)

# 9. Figures
plt.rcParams.update({"font.size": 9, "axes.spines.top": False, "axes.spines.right": False})

# Figure 1 (schematic): four-layer measurement architecture
fig, ax = plt.subplots(figsize=(9.6, 5.5)); ax.axis("off")
ax.set_xlim(0, 10); ax.set_ylim(0, 6)
layers = [("LAYER 1\nEnabling conditions", "Income, prices,\ninfrastructure, services"),
          ("LAYER 2\nPractice", "Realised behaviour\nunder constraint,\nnot self-report"),
          ("LAYER 3\nPersistence", "Behaviour sustained\nacross years,\nnot one-off"),
          ("LAYER 4\nEnvironmental consequence", "Emissions and impacts,\nper capita and\npopulation scale")]
for i, (title, sub) in enumerate(layers):
    x = 0.4 + i * 2.45
    ax.add_patch(FancyBboxPatch((x, 2.2), 2.0, 1.9, boxstyle="round,pad=0.08",
                                fc="#DCE9F5", ec="#1565C0", lw=1.4))
    ax.text(x + 1.0, 3.55, title, ha="center", va="center", fontsize=9, fontweight="bold", color="#0D3B66")
    ax.text(x + 1.0, 2.75, sub, ha="center", va="center", fontsize=7.5, color="#333333")
    if i < 3:
        ax.add_patch(FancyArrowPatch((x + 2.12, 3.15), (x + 2.42, 3.15),
                                     arrowstyle="-|>", mutation_scale=16, color="#555555"))
ax.plot([0.4, 4.85], [4.5, 4.5], color=C_OVER, lw=2)
ax.text(2.6, 4.75, "What participation-based indicators observe", ha="center", color=C_OVER, fontsize=9, fontweight="bold")
ax.plot([5.0, 9.85], [4.5, 4.5], color=C_TEAL, lw=2)
ax.text(7.4, 4.75, "Where realised sustainability actually lives", ha="center", color=C_TEAL, fontsize=9, fontweight="bold")
ax.annotate("constraint filter\n(expectation from means)", xy=(3.2, 2.05), xytext=(1.6, 1.0),
            fontsize=8, color="#555555", arrowprops=dict(arrowstyle="->", ls="--", color="#555555"))
ax.annotate("dual reporting\n(per capita + population scale)", xy=(7.9, 2.05), xytext=(6.2, 1.0),
            fontsize=8, color="#555555", arrowprops=dict(arrowstyle="->", ls="--", color="#555555"))
ax.text(5.0, 0.35, "Schematic of the study design. Country weights and thresholds are fixed before analysis "
                   "(design frozen); no population bonus and no deprivation credit are applied at any layer.",
        ha="center", fontsize=8, color="#555555")
fig.savefig("fig1_architecture.png", dpi=300, bbox_inches="tight"); plt.close(fig)

# Figure 2: the measurement space
fig, ax = plt.subplots(figsize=(7.2, 7.2))
xs = np.linspace(df["gdp_pc"].min(), df["gdp_pc"].max(), 200)
grid = pd.DataFrame({"gdp_pc": xs,
                     "elec": np.full_like(xs, df["elec"].mean()),
                     "urb": np.full_like(xs, df["urb"].mean())})
gp = m.get_prediction(design(grid)).summary_frame(ALPHA)
ax.plot(xs, np.exp(gp["mean"]), color="#333333", lw=1.6, label="Expectation (at mean covariates)")
ax.fill_between(xs, np.exp(gp["obs_ci_lower"]), np.exp(gp["obs_ci_upper"]), color="#BBBBBB", alpha=0.35, label="95% prediction band")
thresh_gdp = np.exp((np.log(THRESH_T) - m.params["const"] - m.params["elec"]*df["elec"].mean()
                     - m.params["urb"]*df["urb"].mean()) / m.params["ln_gdp"])
ax.axvline(thresh_gdp, color=C_OVER, ls="--", lw=1.4)
for cls, col in [("Blindness zone", C_ZONE), ("On expectation", C_ON), ("Above expectation", C_OVER)]:
    sub = df[df["class"] == cls]
    ax.scatter(sub["gdp_pc"], sub["co2_cons_pc"], s=26, color=col, alpha=0.85, edgecolor="white", lw=0.4,
               label=f"{cls} (n={len(sub)})", zorder=3)
ax.set_xscale("log"); ax.set_yscale("log")
ax.set_xlabel("GDP per capita, PPP (current international $, log scale)")
ax.set_ylabel("Consumption-based CO2 per capita, 2023 (t, log scale)")
ax.legend(loc="upper left", frameon=False, fontsize=8)
fig.savefig("fig2_measurement_space.png", dpi=300, bbox_inches="tight"); plt.close(fig)

# Figure 3: persistence of position
fig, ax = plt.subplots(figsize=(7.0, 7.0))
ax.scatter(P["r2000"], P["r2023"], s=24, color=C_ON, alpha=0.8, edgecolor="white", lw=0.4)
lim = [-1.4, 1.9]
ax.plot(lim, lim, ls="--", color="#555555", lw=1.2)
ax.axhline(0, color="#999999", lw=0.8); ax.axvline(0, color="#999999", lw=0.8)
ax.set_xlim(lim); ax.set_ylim(lim)
ax.set_xlabel("Deviation from expectation, 2000 (log points)")
ax.set_ylabel("Deviation from expectation, 2023 (log points)")
ax.set_title("Spearman rho = 0.43 across 23 years (0.65 for 2000-2010, 0.64 for 2010-2023), n = 117", fontsize=9)
fig.savefig("fig3_persistence.png", dpi=300, bbox_inches="tight"); plt.close(fig)

# Figure 4: dual reporting slope chart (log scales)
fig, (axl, axr) = plt.subplots(1, 2, figsize=(8.4, 8.0), sharey=False)
t15 = top15.sort_values("co2_cons_pc")
yl = t15["co2_cons_pc"].values; yr_ = t15["total_mt"].values
for k, (a, b) in enumerate(zip(yl, yr_)):
    col = C_OVER if tot_rank[t15.index[k]] <= 10 and pc_rank[t15.index[k]] > 20 else "#999999"
    axl.plot([0, 1], [np.log10(a), np.log10(b)], color=col, lw=1.4, alpha=0.9)
axl.scatter(np.zeros_like(yl), np.log10(yl), color=C_ON, s=30, zorder=3)
axl.scatter(np.ones_like(yl), np.log10(yr_), color=C_OVER, s=30, zorder=3)
for k, e in enumerate(t15["Entity"]):
    axl.text(-0.04, np.log10(yl[k]), f"{e}  {yl[k]:.2f} t", ha="right", va="center", fontsize=7.5)
    axl.text(1.04, np.log10(yr_[k]), f"{yr_[k]/1000:.2f} Gt", ha="left", va="center", fontsize=7.5)
axl.set_xlim(-0.55, 1.75); axl.set_xticks([0, 1])
axl.set_xticklabels(["Per-capita consumption CO2, 2023\n(t, log scale)",
                     "Population-scale consumption CO2, 2023\n(Gt, log scale)"])
yt = np.arange(np.log10(yl).min(), np.log10(yr_).max() + 0.1, 0.5)
axl.set_yticks(yt); axl.set_yticklabels([f"{10**v:.2g}" for v in yt])
axr.axis("off")
fig.savefig("fig4_dual_reporting.png", dpi=300, bbox_inches="tight"); plt.close(fig)

print("\nDone. All estimates, Tables 1/S1/S2/S3 and Figures 1-4 regenerated from public sources.")