Climate Change Trend Analysis and Forecasting¶

raw_df = pd.read_csv(DATA_PATH)

print("Shape:", raw_df.shape)
print("\nFirst 10 rows:")
raw_df.head(10)

print("Column names ({} total):".format(len(raw_df.columns)))
print(list(raw_df.columns))
print("\nData types:")
raw_df.dtypes

null_pct = (raw_df.isna().sum() / len(raw_df) * 100).round(2).sort_values(ascending=False)
null_report = null_pct.to_frame("pct_null_of_total_rows")
print("Columns with the most missing data:")
null_report.head(15)

# Which countries/years have the most complete co2 coverage?
coverage = (
    raw_df.dropna(subset=["co2"])
    .groupby("country")["year"]
    .agg(first_year="min", last_year="max", n_years="count")
    .sort_values("n_years", ascending=False)
)
print("Most complete CO2 coverage (top 15):")
coverage.head(15)

# Rows per year, to see which years have the broadest country coverage
rows_per_year = raw_df.dropna(subset=["co2"]).groupby("year")["country"].nunique()
print("Number of countries reporting CO2, most recent years:")
rows_per_year.tail(10)

df = raw_df[(raw_df["year"] >= 1990) & (raw_df["iso_code"].notna())].copy()

print(f"Rows before filtering: {len(raw_df):,}")
print(f"Rows after filtering:  {len(df):,}")
print(f"Countries retained:    {df['country'].nunique()}")
print(f"Aggregate rows removed (no iso_code): {raw_df['iso_code'].isna().sum():,}")

global_co2 = df.groupby("year")["co2"].sum().reset_index()

fig = px.line(
    global_co2, x="year", y="co2",
    title="Global CO2 Emissions by Sovereign Countries, 1990-2024",
    labels={"co2": "CO2 Emissions (Mt)", "year": "Year"},
)
fig.update_traces(line_color="#1f77b4", line_width=3)
fig.update_layout(template="plotly_white")
fig.show()

top5_emitters = ["China", "United States", "India", "Russia", "Japan"]
top5_df = df[df["country"].isin(top5_emitters)]

fig = px.line(
    top5_df, x="year", y="co2", color="country",
    title="CO2 Emission Trends for the Top 5 Emitting Countries, 1990-2024",
    labels={"co2": "CO2 Emissions (Mt)", "year": "Year", "country": "Country"},
    color_discrete_map=COUNTRY_COLORS,
)
fig.update_traces(line_width=2.5)
fig.update_layout(template="plotly_white", legend_title_text="Country")
fig.show()

decade_df = df.copy()
decade_df["decade"] = (decade_df["year"] // 10 * 10).astype(str) + "s"
decade_df = decade_df[decade_df["decade"].isin(["1990s", "2000s", "2010s", "2020s"])]

gas_by_decade = (
    decade_df.groupby("decade")[["co2", "methane", "nitrous_oxide"]]
    .sum()
    .rename(columns={"co2": "CO2", "methane": "Methane (CH4)", "nitrous_oxide": "Nitrous Oxide (N2O)"})
)
gas_share = gas_by_decade.div(gas_by_decade.sum(axis=1), axis=0) * 100
gas_share_long = gas_share.reset_index().melt(id_vars="decade", var_name="Gas", value_name="Share (%)")

fig = px.bar(
    gas_share_long, x="decade", y="Share (%)", color="Gas",
    title="Share of Total Global GHG Emissions by Gas Type, per Decade",
    labels={"decade": "Decade"},
    color_discrete_sequence=px.colors.qualitative.Safe,
)
fig.update_layout(template="plotly_white", barmode="stack")
fig.show()

base = (
    df[df["country"].isin(COUNTRIES)]
    .sort_values(["country", "year"])
    .reset_index(drop=True)
    .copy()
)
print(f"Rows: {len(base)}  |  Countries: {base['country'].nunique()}  |  "
      f"Years: {base['year'].min()}-{base['year'].max()}")
base[["country", "year", "co2", "population", "gdp"]].head()

base["decade"] = (base["year"] // 10 * 10).astype(str) + "s"
base["years_since_1990"] = base["year"] - 1990

# 5-year rolling average of CO2, computed independently per country
base["co2_5yr_rolling_mean"] = (
    base.groupby("country")["co2"]
    .transform(lambda s: s.rolling(window=5, min_periods=1).mean())
)

base[["country", "year", "decade", "years_since_1990", "co2", "co2_5yr_rolling_mean"]].head(8)

for lag in (1, 2, 3):
    base[f"co2_lag{lag}"] = base.groupby("country")["co2"].shift(lag)

base[["country", "year", "co2", "co2_lag1", "co2_lag2", "co2_lag3"]].head(8)

# Cross-check co2_per_capita = co2 (Mt) * 1e6 / population, for 3 countries x 3 years
check_years = [1990, 2010, 2023]
check_countries = ["China", "United States", "India"]

check_rows = []
for country in check_countries:
    for y in check_years:
        row = base[(base["country"] == country) & (base["year"] == y)]
        if row.empty:
            continue
        computed = row["co2"].values[0] * 1e6 / row["population"].values[0]
        reported = row["co2_per_capita"].values[0]
        check_rows.append({
            "country": country, "year": y,
            "computed_per_capita": round(computed, 3),
            "reported_per_capita": round(reported, 3),
            "diff": round(computed - reported, 4),
        })

pd.DataFrame(check_rows)

base["ghg_intensity"] = np.where(
    base["gdp"].notna() & (base["gdp"] > 0),
    base["total_ghg"] / base["gdp"],
    np.nan,
)

missing_intensity = base[base["ghg_intensity"].isna()][["country", "year"]]
print(f"Country-years where ghg_intensity could not be computed (missing/zero GDP): {len(missing_intensity)}")
missing_intensity

base["co2_yoy_change"] = base.groupby("country")["co2"].diff()
base["co2_yoy_pct_change"] = base.groupby("country")["co2"].pct_change() * 100

avg_growth = (
    base[base["year"] >= 1991]
    .groupby("country")["co2_yoy_pct_change"]
    .mean()
    .sort_values(ascending=False)
)

print("Top 5 countries by highest average annual CO2 growth rate since 1990:")
print(avg_growth.head(5).round(2))

co2_1990 = base[base["year"] == 1990].set_index("country")["co2"]
co2_latest = base[base["year"] == base["year"].max()].set_index("country")["co2"]
total_reduction = (co2_1990 - co2_latest).sort_values(ascending=False)

print(f"Top 5 countries by largest absolute CO2 reduction since 1990 (1990 level minus "
      f"{base['year'].max()} level, Mt):")
print(total_reduction.head(5).round(1))

feature_cols = [
    "country", "year", "co2", "co2_per_capita", "co2_5yr_rolling_mean",
    "co2_lag1", "co2_lag2", "co2_lag3", "co2_yoy_pct_change", "ghg_intensity",
]
features_df = base[feature_cols].copy()

features_df.to_csv("../data/ghg_features.csv", index=False)
print(f"Saved ../data/ghg_features.csv with shape {features_df.shape}")
features_df.head()

model_df = base.copy()
model_df["target_co2_next"] = model_df.groupby("country")["co2"].shift(-1)

feature_names = [
    "years_since_1990", "co2", "co2_per_capita",
    "co2_5yr_rolling_mean", "co2_lag1", "co2_lag2", "co2_lag3",
]

# Drop rows without a full lag history or without a next-year target to predict
model_df = model_df.dropna(subset=feature_names + ["target_co2_next"]).reset_index(drop=True)
print(f"Modelling rows: {len(model_df)}  |  Year range: {model_df['year'].min()}-{model_df['year'].max()}")
model_df[["country", "year"] + feature_names + ["target_co2_next"]].head()

train = model_df[model_df["year"] <= 2018].copy()
test = model_df[model_df["year"] >= 2019].copy()

split_counts = pd.DataFrame({
    "train_rows": train.groupby("country").size(),
    "test_rows": test.groupby("country").size(),
}).fillna(0).astype(int)
split_counts

# Naive baseline: predict next year's CO2 = this year's CO2 (no-change model)
test_baseline = test.copy()
test_baseline["pred_naive"] = test_baseline["co2"]

baseline_results = []
for country in COUNTRIES:
    t = test_baseline[test_baseline["country"] == country]
    if t.empty:
        continue
    mae = mean_absolute_error(t["target_co2_next"], t["pred_naive"])
    rmse = mean_squared_error(t["target_co2_next"], t["pred_naive"]) ** 0.5
    baseline_results.append({"country": country, "baseline_mae": mae, "baseline_rmse": rmse})

baseline_df = pd.DataFrame(baseline_results)
baseline_df

fig = make_subplots(rows=1, cols=3, subplot_titles=["China", "United States", "Germany"])
for i, country in enumerate(["China", "United States", "Germany"], start=1):
    t = test_baseline[test_baseline["country"] == country].sort_values("year")
    fig.add_trace(go.Scatter(x=t["year"] + 1, y=t["target_co2_next"], name=f"{country} actual",
                              mode="lines+markers", line=dict(color="black")), row=1, col=i)
    fig.add_trace(go.Scatter(x=t["year"] + 1, y=t["pred_naive"], name=f"{country} naive pred",
                              mode="lines+markers", line=dict(color="crimson", dash="dash")), row=1, col=i)
fig.update_layout(title="Naive Baseline: Actual vs Predicted CO2 (Test Set)", template="plotly_white",
                   showlegend=False, height=400)
fig.show()

lr_models = {}
lr_results = []
lr_test_preds = []

for country in COUNTRIES:
    train_c = train[train["country"] == country]
    test_c = test[test["country"] == country]
    if len(train_c) < 3 or test_c.empty:
        continue

    X_train, y_train = train_c[feature_names], train_c["target_co2_next"]
    X_test, y_test = test_c[feature_names], test_c["target_co2_next"]

    lr = LinearRegression().fit(X_train, y_train)
    preds = lr.predict(X_test)

    mae = mean_absolute_error(y_test, preds)
    rmse = mean_squared_error(y_test, preds) ** 0.5

    lr_models[country] = lr
    lr_results.append({"country": country, "lr_mae": mae, "lr_rmse": rmse})
    lr_test_preds.append(test_c.assign(pred_lr=preds))

lr_df = pd.DataFrame(lr_results)
lr_test_preds_df = pd.concat(lr_test_preds, ignore_index=True)
lr_df

fig = make_subplots(rows=1, cols=3, subplot_titles=["China", "United States", "Germany"])
for i, country in enumerate(["China", "United States", "Germany"], start=1):
    t = lr_test_preds_df[lr_test_preds_df["country"] == country].sort_values("year")
    fig.add_trace(go.Scatter(x=t["year"] + 1, y=t["target_co2_next"], name=f"{country} actual",
                              mode="lines+markers", line=dict(color="black")), row=1, col=i)
    fig.add_trace(go.Scatter(x=t["year"] + 1, y=t["pred_lr"], name=f"{country} LR pred",
                              mode="lines+markers", line=dict(color="royalblue", dash="dash")), row=1, col=i)
fig.update_layout(title="Linear Regression: Actual vs Predicted CO2 (Test Set)", template="plotly_white",
                   showlegend=False, height=400)
fig.show()

coef_table = pd.DataFrame({
    country: model.coef_ for country, model in lr_models.items()
}, index=feature_names).T
coef_table["intercept"] = [model.intercept_ for model in lr_models.values()]
coef_table.round(3)

train_pool = train.copy()
test_pool = test.copy()

country_encoder = LabelEncoder().fit(model_df["country"])
train_pool["country_encoded"] = country_encoder.transform(train_pool["country"])
test_pool["country_encoded"] = country_encoder.transform(test_pool["country"])

rf_feature_names = feature_names + ["country_encoded"]

rf = RandomForestRegressor(n_estimators=100, random_state=42)
rf.fit(train_pool[rf_feature_names], train_pool["target_co2_next"])

print(f"Random Forest trained on {len(train_pool)} pooled rows across {train_pool['country'].nunique()} countries")

rf_results = []
rf_test_preds = []

for country in COUNTRIES:
    test_c = test_pool[test_pool["country"] == country]
    if test_c.empty:
        continue
    preds = rf.predict(test_c[rf_feature_names])
    mae = mean_absolute_error(test_c["target_co2_next"], preds)
    rmse = mean_squared_error(test_c["target_co2_next"], preds) ** 0.5
    rf_results.append({"country": country, "rf_mae": mae, "rf_rmse": rmse})
    rf_test_preds.append(test_c.assign(pred_rf=preds))

rf_df = pd.DataFrame(rf_results)
rf_test_preds_df = pd.concat(rf_test_preds, ignore_index=True)
rf_df

importance_df = pd.DataFrame({
    "feature": rf_feature_names,
    "importance": rf.feature_importances_,
}).sort_values("importance", ascending=True)

fig = px.bar(
    importance_df, x="importance", y="feature", orientation="h",
    title="Pooled Random Forest — Feature Importance (single model, all 10 countries)",
    labels={"importance": "Importance", "feature": "Feature"},
    color_discrete_sequence=["#2ca02c"],
)
fig.update_layout(template="plotly_white")
fig.show()

comparison_df = (
    baseline_df.merge(lr_df, on="country").merge(rf_df, on="country")
)
comparison_df = comparison_df[
    ["country", "baseline_mae", "lr_mae", "rf_mae", "baseline_rmse", "lr_rmse", "rf_rmse"]
]
comparison_df.columns = [
    "Country", "Baseline MAE", "LR MAE", "RF MAE", "Baseline RMSE", "LR RMSE", "RF RMSE",
]
comparison_df["Best Model"] = comparison_df[["Baseline MAE", "LR MAE", "RF MAE"]].idxmin(axis=1).str.replace(" MAE", "")

comparison_df.round(2)

train_years = range(1990, 2019)

ets_models = {}
for country in COUNTRIES:
    series = base[(base["country"] == country) & (base["year"].isin(train_years))].sort_values("year")
    y = pd.Series(series["co2"].values, index=series["year"].values)

    model = ExponentialSmoothing(y, trend="add", damped_trend=True, seasonal=None)
    ets_models[country] = model.fit(optimized=True)

print(f"Fitted ETS(A,Ad,N) models for {len(ets_models)} countries on the 1990-2018 training window.")

param_rows = []
for country in ["United Kingdom", "Germany", "India", "China"]:
    p = ets_models[country].params
    param_rows.append({
        "country": country,
        "alpha (level)": round(p["smoothing_level"], 3),
        "beta* (trend)": round(p["smoothing_trend"], 3),
        "phi (damping)": round(p["damping_trend"], 3),
    })
pd.DataFrame(param_rows)

FORECAST_YEARS = list(range(2019, 2044))  # 2019-2023 holdout + 2024-2043 out-of-sample
N_STEPS = len(FORECAST_YEARS)

ets_forecasts = {}
for country in COUNTRIES:
    fit = ets_models[country]
    point_forecast = fit.forecast(N_STEPS)
    point_forecast.index = FORECAST_YEARS

    # HoltWintersResults has no get_forecast(); build a 95% CI via bootstrap simulation instead
    sims = fit.simulate(nsimulations=N_STEPS, repetitions=1000, error="add",
                         random_errors="bootstrap", random_state=RANDOM_STATE)
    sims.index = FORECAST_YEARS

    ets_forecasts[country] = {
        "mean": point_forecast,
        "ci_lower": sims.quantile(0.025, axis=1),
        "ci_upper": sims.quantile(0.975, axis=1),
        "fitted": fit.fittedvalues,
    }

print("Generated 1990-2018 fitted values plus 2019-2043 forecasts (with 95% CI) for all 10 countries.")

def plot_country_forecast(country):
    hist = base[(base["country"] == country) & (base["year"] <= 2018)].sort_values("year")
    holdout = base[(base["country"] == country) & (base["year"].between(2019, 2023))].sort_values("year")
    fc = ets_forecasts[country]

    fig = go.Figure()
    fig.add_trace(go.Scatter(x=hist["year"], y=hist["co2"], name="Historical actual (1990-2018)",
                              mode="lines", line=dict(color="black", width=2)))
    fig.add_trace(go.Scatter(x=fc["fitted"].index, y=fc["fitted"].values, name="Fitted values",
                              mode="lines", line=dict(color="gray", dash="dot")))
    fig.add_trace(go.Scatter(x=holdout["year"], y=holdout["co2"], name="Holdout actual (2019-2023)",
                              mode="markers", marker=dict(color="black", size=8, symbol="diamond")))
    fig.add_trace(go.Scatter(x=fc["mean"].index, y=fc["mean"].values, name="Forecast to 2043",
                              mode="lines", line=dict(color="crimson", width=2)))
    fig.add_trace(go.Scatter(
        x=list(fc["ci_upper"].index) + list(fc["ci_lower"].index[::-1]),
        y=list(fc["ci_upper"].values) + list(fc["ci_lower"].values[::-1]),
        fill="toself", fillcolor="rgba(220,20,60,0.15)", line=dict(color="rgba(0,0,0,0)"),
        name="95% CI", hoverinfo="skip",
    ))
    fig.update_layout(title=f"ETS(A,Ad,N) Forecast: {country}", xaxis_title="Year", yaxis_title="CO2 Emissions (Mt)",
                       template="plotly_white", legend=dict(orientation="h", y=-0.2))
    return fig


for country in COUNTRIES:
    plot_country_forecast(country).show()

summary_rows = []
for country in COUNTRIES:
    mean_fc = ets_forecasts[country]["mean"]
    actual_2020 = base.loc[(base["country"] == country) & (base["year"] == 2020), "co2"].values[0]
    f2030, f2035, f2040 = mean_fc.loc[2030], mean_fc.loc[2035], mean_fc.loc[2040]
    pct_change = (f2040 - actual_2020) / actual_2020 * 100
    summary_rows.append({
        "Country": country,
        "2030 Forecast (MtCO2)": round(f2030, 1),
        "2035 Forecast (MtCO2)": round(f2035, 1),
        "2040 Forecast (MtCO2)": round(f2040, 1),
        "2020 Actual (MtCO2)": round(actual_2020, 1),
        "% Change 2020->2040": round(pct_change, 1),
    })

forecast_summary_df = pd.DataFrame(summary_rows)
forecast_summary_df

ets_eval = []
for country in COUNTRIES:
    actual_holdout = (
        base[(base["country"] == country) & (base["year"].between(2019, 2023))]
        .sort_values("year")["co2"].values
    )
    pred_holdout = ets_forecasts[country]["mean"].loc[2019:2023].values
    mae = mean_absolute_error(actual_holdout, pred_holdout)
    rmse = mean_squared_error(actual_holdout, pred_holdout) ** 0.5
    ets_eval.append({"country": country, "ets_mae": mae, "ets_rmse": rmse})

ets_eval_df = pd.DataFrame(ets_eval)
ets_eval_df.round(2)

four_model_df = comparison_df.merge(
    ets_eval_df.rename(columns={"country": "Country", "ets_mae": "ETS MAE", "ets_rmse": "ETS RMSE"}),
    on="Country",
)
four_model_df = four_model_df[
    ["Country", "Baseline MAE", "LR MAE", "RF MAE", "ETS MAE",
     "Baseline RMSE", "LR RMSE", "RF RMSE", "ETS RMSE", "Best Model"]
]
four_model_df.round(2)

scenario_rows = []
scenario_years = list(range(2025, 2041))

for country in COUNTRIES:
    base_2024 = base.loc[(base["country"] == country) & (base["year"] == 2024), "co2"].values[0]
    bau_mean = ets_forecasts[country]["mean"]

    for year in scenario_years:
        n = year - 2024
        scenario_rows.append({"country": country, "year": year, "scenario": "BAU",
                               "co2_projected": bau_mean.loc[year]})
        scenario_rows.append({"country": country, "year": year, "scenario": "Moderate Mitigation",
                               "co2_projected": base_2024 * (1 - 0.02) ** n})
        scenario_rows.append({"country": country, "year": year, "scenario": "Aggressive Mitigation",
                               "co2_projected": base_2024 * (1 - 0.05) ** n})

scenario_df = pd.DataFrame(scenario_rows)
scenario_df.to_csv("../data/scenario_projections.csv", index=False)
print(f"Saved ../data/scenario_projections.csv with shape {scenario_df.shape}")
scenario_df.head(9)

def plot_country_scenarios(country):
    hist = base[(base["country"] == country) & (base["year"].between(1990, 2024))].sort_values("year")
    level_1990 = base.loc[(base["country"] == country) & (base["year"] == 1990), "co2"].values[0]

    fig = go.Figure()
    fig.add_trace(go.Scatter(x=hist["year"], y=hist["co2"], name="Historical actual (1990-2024)",
                              mode="lines", line=dict(color="lightgray", width=3)))
    for scen, color in SCENARIO_COLORS.items():
        s = scenario_df[(scenario_df["country"] == country) & (scenario_df["scenario"] == scen)].sort_values("year")
        fig.add_trace(go.Scatter(x=s["year"], y=s["co2_projected"], name=scen,
                                  mode="lines", line=dict(color=color, width=2.5)))
    fig.add_hline(y=level_1990, line_dash="dash", line_color="black",
                  annotation_text="1990 level (policy benchmark)", annotation_position="top left")
    fig.update_layout(title=f"Emissions Scenarios, 2020-2040: {country}", xaxis_title="Year",
                       yaxis_title="CO2 Emissions (Mt)", template="plotly_white",
                       legend=dict(orientation="h", y=-0.2))
    fig.update_xaxes(range=[2020, 2040])
    return fig


for country in COUNTRIES:
    plot_country_scenarios(country).show()

global_scenario = scenario_df.groupby(["scenario", "year"])["co2_projected"].sum().reset_index()

fig = px.line(
    global_scenario, x="year", y="co2_projected", color="scenario",
    title="Global Aggregate Emissions Scenarios (Sum of All 10 Countries), 2025-2040",
    labels={"co2_projected": "CO2 Emissions (Mt)", "year": "Year", "scenario": "Scenario"},
    color_discrete_map=SCENARIO_COLORS,
)
fig.update_traces(line_width=3)
fig.update_layout(template="plotly_white")
fig.show()

cumulative_df = (
    scenario_df.groupby(["country", "scenario"])["co2_projected"]
    .sum()
    .reset_index()
    .rename(columns={"co2_projected": "cumulative_co2_2025_2040"})
)

fig = px.bar(
    cumulative_df, x="country", y="cumulative_co2_2025_2040", color="scenario",
    barmode="group",
    title="Cumulative CO2 Emissions by Country and Scenario, 2025-2040",
    labels={"cumulative_co2_2025_2040": "Cumulative CO2 (Mt)", "country": "Country", "scenario": "Scenario"},
    color_discrete_map=SCENARIO_COLORS,
)
fig.update_layout(template="plotly_white", xaxis_tickangle=-30)
fig.show()