二氧化碳與全球升溫：由曲線擬合到可靠外推

封存文章

本文源自 HiMCM 2022 的建模練習。舊版本包含單位誤讀、錯誤外推輸入，並以樣本內 $R^2$ 選擇百年預測模型；本版本保留原有比較問題，但把結果改寫為可檢查的統計流程。它不是新的氣候情景預測。

工業革命前，大氣二氧化碳濃度約為 280 ppm；現代連續觀測顯示濃度長期上升。 HiMCM 題目要求學生使用提供的 $\mathrm{CO_2}$ 與 land–ocean temperature anomaly 資料回答三類問題：

2004 年前十年的增幅是否歷來最大？
不同趨勢模型如何外推至 2050 或 2100？
$\mathrm{CO_2}$ 與溫度異常有何統計關係，而這種關係可外推多遠？

這些看似是 curve fitting 問題，實際上更考驗資料對齊、時間依賴、模型選擇與外推誠信。

一、資料與前處理

原始工作簿包含年度 $\mathrm{CO_2}$ 濃度與全球 land–ocean temperature anomaly。讀入後應先：

確認年份是否為 calendar year；
確認濃度是年度平均、月平均還是某站觀測；
確認溫度異常的基準期；
檢查缺失值與資料修訂；
只在共同年份比較兩組序列。

import numpy as np
import pandas as pd

co2 = pd.read_excel(
    "2022_HiMCM_Data.xlsx",
    sheet_name="CO2 Data Set 1",
).rename(columns={"PPM": "co2_ppm"})

temp = pd.read_excel(
    "2022_HiMCM_Data.xlsx",
    sheet_name="Temps Data Set 2",
).rename(columns={"Degrees C": "temp_anomaly_c"})

joined = (
    co2[["Year", "co2_ppm"]]
    .merge(temp[["Year", "temp_anomaly_c"]], on="Year", how="inner")
    .dropna()
    .sort_values("Year")
)

原資料開頭大致如下：

年份	$\mathrm{CO_2}$ （ppm）
1959	315.97
1960	316.91
1961	317.64
1962	318.45
1963	318.99

溫度欄位是相對基準期的異常值，不是全球平均絕對溫度。

二、十年增幅

「2004 年增幅最大」至少有三種解讀：

1994 至 2004 的端點差；
截至 2004 年的十年線性趨勢；
2004 年附近月度資料的十年增長率。

若題目資料是年度值，可先計算 rolling endpoint change：

co2 = co2.sort_values("Year").copy()
co2["increase_10y"] = co2["co2_ppm"] - co2["co2_ppm"].shift(10)
largest = co2.loc[co2["increase_10y"].idxmax()]

舊分析報告 1994–2004 約增加 13.82 ppm，而 1998–2008 約增加 18.99 ppm，因此「截至 2004 年比資料內任何其他十年都大」並不成立。不過若原聲稱是「截至當時」或使用月度序列，便要按同一觀測頻率和資料截止日重做，不能用後來十年直接改寫歷史語境。

三、三種 $\mathrm{CO_2}$ 趨勢

為減少年份數值過大造成的病態，令

\tau=t-t_0.

比較：

\begin{aligned} \text{Linear:}\quad&C(t)=\beta_0+\beta_1\tau,\\ \text{Exponential:}\quad&\log C(t)=\gamma_0+\gamma_1\tau,\\ \text{Quadratic:}\quad&C(t)=\alpha_0+\alpha_1\tau+\alpha_2\tau^2. \end{aligned}

from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import mean_squared_error

t0 = co2["Year"].min()
X = (co2["Year"] - t0).to_numpy().reshape(-1, 1)
y = co2["co2_ppm"].to_numpy()

linear = LinearRegression().fit(X, y)
quadratic = make_pipeline(
    PolynomialFeatures(2, include_bias=False),
    LinearRegression(),
).fit(X, y)
exponential = LinearRegression().fit(X, np.log(y))

舊程式得到的樣本內結果是：

模型	樣本內 RMSE（ppm）	樣本內 $R^2$	2100 外推（ppm）
線性	3.92	0.98	534.88
指數	2.89	0.99	583.70
二次	0.72	約 1.00	688.33

這張表只重現舊運算，不能據此宣布二次模型「最準確」。加入一個二次項幾乎必然改善樣本內 fit，但百年外推由最高次項支配，微小係數誤差會被放大。應使用 rolling-origin validation：

def rolling_rmse(model_factory, years, values, min_train=25):
    errors = []
    origin = years.min()
    for stop in range(min_train, len(values)):
        X_train = (years[:stop] - origin).reshape(-1, 1)
        X_test = np.array([[years[stop] - origin]])
        model = model_factory().fit(X_train, values[:stop])
        errors.append(values[stop] - model.predict(X_test)[0])
    return np.sqrt(np.mean(np.square(errors)))

還應檢查殘差自相關、參數穩定性，以及改變 training window 後的 forecast sensitivity。

四、何時到達 685 ppm？

對線性模型，若 $\beta_1>0$ ，

t_{685}=t_0+\frac{685-\beta_0}{\beta_1}.

對指數模型，

t_{685} =t_0+\frac{\log685-\gamma_0}{\gamma_1}.

對二次模型則解

\alpha_2\tau^2+\alpha_1\tau+\alpha_0-685=0

並只保留位於資料末年之後的實根。舊程式沒有正確求指數模型的 crossing，因而把它報成 inf；這是實作缺漏，不是模型永不到達 685 ppm 的證據。

更重要的是， $\mathrm{CO_2}$ 未來濃度取決於排放、碳循環與政策路徑。純時間曲線沒有情景變數，所以即使求出 crossing year，也只表示「若歷史數學形狀繼續」，而不是政策不變或物理系統必然如此。

五、溫度趨勢

可同樣擬合時間與 temperature anomaly：

T(t)=a_0+a_1\tau

以及

T(t)=b_0+b_1\tau+b_2\tau^2.

舊運算的樣本內數值為：

模型	RMSE（°C）	$R^2$
線性	0.11	0.89
二次	0.09	0.92

並外推不同 threshold year：

相對基準期異常	線性	二次
1.25 °C	2044	2030
1.50 °C	2060	2038
2.00 °C	2090	2052

兩種簡單模型給出差距很大的年份，正好顯示 extrapolation uncertainty，而不是讓我們挑一條看來最合理的曲線。年度溫度還受火山、ENSO、氣溶膠與自然變率影響；殘差也不是獨立同分布。

六、 $\mathrm{CO_2}$ 與溫度的統計關係

最簡單的 contemporaneous regression 是

T_t=\theta_0+\theta_1C_t+\varepsilon_t.

Xc = joined[["co2_ppm"]].to_numpy()
yt = joined["temp_anomaly_c"].to_numpy()
co2_temp = LinearRegression().fit(Xc, yt)

effect_per_100ppm = 100 * co2_temp.coef_[0]

舊程式直接打印 coef_[0]，卻把文字標成「每增加 100 ppm」；如果係數單位是 °C/ppm，必須乘 100。原文所稱「每 100 ppm 只升 0.01 °C」因此是單位錯誤。

更嚴重的是：

co2_temp_linear_model.predict([[2100]])

把 2100 當成 ppm 輸入，而不是年份，然後把輸出 18.61 °C 誤標為「2100 年預測」。二次版本的 27.51 °C 同樣無效，兩者應從文章結論移除。

即使修正程式，高 $R^2$ 仍不等於因果效應。兩條具有共同時間趨勢的非平穩序列可產生 spurious regression。最低限度應：

檢查差分或 detrended 關係；
考慮 ocean thermal inertia 與 distributed lag；
加入其他 forcing 或使用由物理模型支持的 energy-balance formulation；
對 residual autocorrelation 使用合適的不確定區間。

例如先探索年度變化：

\Delta T_t =\eta_0+\eta_1\Delta\log C_t+u_t,

但這也只是一個診斷模型，不取代 attribution study。

七、模型可靠範圍

這個練習可以支持：

描述資料期間內的長期上升；
比較短期 holdout 下不同曲線的預測誤差；
顯示 2050／2100 外推對函數形式極敏感；
解釋共同趨勢與因果推論的分別。

它不能單獨支持：

確定某個 warming threshold 的到達年份；
由相關性估計 equilibrium climate sensitivity；
對政策、排放或碳循環改變作反事實推論；
把二次曲線無限延伸。

若要作更完整預測，應由 emission scenario 經 carbon-cycle model 得到濃度，再以 radiative forcing 與 energy-balance／Earth-system model 連接溫度，並報告 ensemble uncertainty。

結論

原始分析的價值在於提出三個值得比較的趨勢模型；它的主要問題則是把樣本內擬合當成長期預測能力，並在 $\mathrm{CO_2}$ –溫度回歸中出現兩個單位／輸入錯誤。

經修正後，最合理的結論是：

2004 的十年增幅聲稱要按相同資料截止日與定義檢查；
二次模型在歷史資料內 fit 最緊，不代表 2100 外推最可信；
不同模型的 threshold year 差異本身就是 model-form uncertainty；
$\mathrm{CO_2}$ 與溫度異常共同上升，但簡單同期回歸不能建立完整因果機制；
所有預測都應附上情景、validation window 和不確定區間。

參考資料

COMAP, Inc. (2022). HiMCM 2022 Problem B: CO₂ and Global Warming and accompanying datasets.
NOAA Global Monitoring Laboratory. Trends in Atmospheric Carbon Dioxide.
NASA Goddard Institute for Space Studies. GISTEMP Surface Temperature Analysis.
IPCC. (2021). Climate Change 2021: The Physical Science Basis. Working Group I contribution to the Sixth Assessment Report.