Lake Mead 水危機：由水量守恆到情景預測

封存文章

本文源自 HiMCM 2021 Problem B。舊版本將不同模型的分歧當成 ensemble 權重，並以逐期 marginal interval 的上下界作累積和，不能形成正確的水位預測區間。本版本保留歷史探索，撤回 2025–2050 的點預測作為管理依據。

Lake Mead 的狀態同時受 Colorado River inflow、下游 release、蒸發、降雨、 bank storage、用水規則和上游水庫操作影響。水位只是儲量的一個非線性轉換；因此只把 elevation 當作單變量時間序列，最多可作短期 benchmark。

研究應分開三個問題：

歷史水位有甚麼趨勢、季節與 drought episode？
一個純統計模型在短期 holdout 上能預測多準？
在不同氣候、需求和管理情景下，reservoir storage 會如何演變？

一、物理水量帳

以月為時間步長：

V_{t+1} =V_t +Q_t^{\rm in} -Q_t^{\rm release} -W_t -E_tA(\bar h_t) +P_tA(\bar h_t) +R_t,

其中：

$V_t$ ：庫容；
$Q^{\rm in}$ ：河流及支流流入；
$Q^{\rm release}$ ：dam release；
$W_t$ ：直接取水；
$E_t$ 、 $P_t$ ：蒸發與降雨深度；
$A(h)$ ：水面面積；
$R_t$ ：bank storage、測量誤差和未建模項。

高程 $h_t$ 由 elevation–storage curve

h_t=g(V_t)

取得。 $g$ 不是線性：同樣一呎水位變化在不同高度對應不同水量。因此以 feet 直接相加 inflow／outflow 在單位上不成立。

二、資料整理

舊 Excel 以 year × month 寬表儲存。轉 long format 時應建立真正月度 index：

import numpy as np
import pandas as pd

MONTHS = {
    "Jan": 1, "Feb": 2, "Mar": 3, "Apr": 4,
    "May": 5, "Jun": 6, "Jul": 7, "Aug": 8,
    "Sep": 9, "Oct": 10, "Nov": 11, "Dec": 12,
}

monthly = (
    raw.replace("---", np.nan)
       .melt(id_vars="Year", var_name="month_name",
             value_name="elevation_ft")
)
monthly["month"] = monthly["month_name"].map(MONTHS)
monthly["date"] = pd.to_datetime(
    dict(year=monthly["Year"], month=monthly["month"], day=1)
)
monthly["elevation_ft"] = pd.to_numeric(
    monthly["elevation_ft"], errors="coerce"
)
monthly = monthly.sort_values("date").set_index("date")

不要一開始便線性插補所有 missing values。先檢查缺失原因；若缺口進入 training， interpolation uncertainty 應保留，若缺口在 validation period，不能用未來觀測填補。

三、歷史描述

可以畫：

月度 elevation；
60-month trailing mean；
每年高、低水位；
monthly climatology；
storage 而非只用 elevation；
inflow、release 和 evaporation 的同期變化。

「drought period」不應只以低於歷史平均判斷。可預先定義：

D_t=\mathbf1(h_t<h_{\rm threshold}),

或用 standardized inflow index／官方 shortage tiers，再報告 duration、depth 和 recovery。若 threshold 是看完圖後才選，episode comparison 會有 selection bias。

四、SARIMA benchmark

月度序列可考慮

\operatorname{SARIMA}(p,d,q)\times(P,D,Q)_{12}.

若有真正 exogenous inputs 才稱 SARIMAX：

\phi(B)\Phi(B^{12}) (1-B)^d(1-B^{12})^Dy_t = \theta(B)\Theta(B^{12})\varepsilon_t +\boldsymbol\beta^\mathsf T\mathbf x_t.

舊模型沒有加入 inflow、release 或 climate covariates，實質只是 SARIMA。而且部分程式先把月度資料按年平均，再仍使用 seasonal period 12；對 annual series，12 表示 12 年季節，不是月份，模型含義已改變。

from statsmodels.tsa.statespace.sarimax import SARIMAX

model = SARIMAX(
    train["elevation_ft"],
    order=(1, 1, 1),
    seasonal_order=(1, 0, 1, 12),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
fit = model.fit(disp=False)

order 不能只靠 auto_arima 在全資料選完再以同一資料評分；應放入每個 training fold。

五、Structural time series

Local linear trend 加 seasonality：

\begin{aligned} y_t&=\mu_t+\gamma_t+\varepsilon_t,\\ \mu_{t+1}&=\mu_t+\nu_t+\eta_t,\\ \nu_{t+1}&=\nu_t+\zeta_t. \end{aligned}

from statsmodels.tsa.statespace.structural import UnobservedComponents

ucm = UnobservedComponents(
    train["elevation_ft"],
    level="local linear trend",
    seasonal=12,
)
ucm_fit = ucm.fit(disp=False)

它容許 trend slope 隨時間改變，但長期外推仍由 state evolution 假設控制，不會自動理解 allocation rules 或 climate change。

六、差分還原

若模型直接預測一階差分

\Delta y_{T+h}=y_{T+h}-y_{T+h-1},

mean path 可還原為

\hat y_{T+h} =y_T+\sum_{j=1}^{h} \widehat{\Delta y}_{T+j}.

但 interval 不能把每一步 lower bounds 累加、upper bounds 累加。差分 forecast 在不同 horizon 之間相關：

\operatorname{Var}(y_{T+h}) = \sum_{i=1}^{h}\sum_{j=1}^{h} \operatorname{Cov} (\Delta y_{T+i},\Delta y_{T+j}).

正確方法是：

讓 state-space model 在 level scale 直接產生 forecast covariance；或
從 joint predictive distribution 模擬整條 future path，再對每個 horizon 取 quantile。

paths = np.stack([
    simulate_future_differences(fit, horizon=360, seed=s)
    for s in range(5000)
])
level_paths = last_level + np.cumsum(paths, axis=1)
median = np.quantile(level_paths, 0.5, axis=0)
lower = np.quantile(level_paths, 0.025, axis=0)
upper = np.quantile(level_paths, 0.975, axis=0)

七、Rolling-origin validation

不能用 1935–2021 全資料 fit，然後只比較 in-sample AIC 便稱 2050 reliable。設多個 forecast origins：

def rolling_origins(series, min_train, horizons):
    for origin in range(min_train, len(series) - max(horizons)):
        train = series.iloc[:origin]
        truth = {
            h: series.iloc[origin + h - 1]
            for h in horizons
        }
        yield train, truth

對每個模型和 horizon 報告：

\mathrm{MAE}_h =\frac1K\sum_k |\hat y_{k,h}-y_{k,h}|,

以及 RMSE、interval coverage、interval width 和 bias。比較 baselines：

last value；
seasonal naive；
drift；
SARIMA；
UCM；
water-balance model。

若複雜模型不能穩定勝過 naive forecast，就不應用於長期政策結論。

八、訓練窗口

全歷史資料涵蓋 reservoir filling、不同操作制度和氣候階段；最近十年則樣本極少且集中於一段 drought。兩者分歧不代表其中一個必然較好，而是顯示 nonstationarity。

可以比較 expanding window、fixed 10/20/30-year windows，但 window length 應以 out-of-sample performance 選擇。只用十年 annual data 卻估計 trend、seasonal 和多個 ARMA 參數，通常不可辨識。

九、舊點預測為何不能保留

舊程式報告 2050：

SARIMAX：1,062.70 ft；
UCM：1,088.30 ft；
recent SARIMAX：1,067.68 ft；
recent UCM：969.96 ft。

跨度超過 118 ft，反映模型結構與窗口主導結果。更重要的是：

annual／monthly frequency 混用；
interval 還原錯誤；
模型沒有 exogenous water drivers；
30 年遠超 validation horizon；
recent model 樣本不足。

因此這些數字是舊 notebook 的探索輸出，不是更新後的 forecast。

十、錯誤的 ensemble

舊 dynamic_weighted_average 用模型彼此的距離計算權重：

error1 = abs(forecast_1 - forecast_2)

模型 disagreement 不是 forecast error，因為未來真值未知。它甚至可能給離群預測不合理權重。Ensemble weights 應由歷史 holdout error 取得，例如：

w_m =\frac{(\mathrm{RMSE}_m+\epsilon)^{-1}} {\sum_\ell(\mathrm{RMSE}_\ell+\epsilon)^{-1}},

並在獨立 validation period 檢查。更穩健可用 constrained regression／stacking，但仍不能把純統計模型變成 policy scenario model。

十一、情景模型

管理問題應在 storage scale 模擬：

V_{t+1}^{(s)} =F\left( V_t^{(s)}, Q_{t}^{\rm in,(s)}, Q_{t}^{\rm release,(s)}, W_t^{(s)}, E_t^{(s)} \right),

情景 $s$ 可包括：

不同 runoff／temperature ensemble；
demand growth 或 conservation；
operating rules；
shortage tiers；
wastewater reuse；
agricultural efficiency；
upstream reservoir coordination。

每項 intervention 要避免 double counting。若回收水本來已間接回到河系，新增 reuse 的 net augmentation 不能等於處理量全數。

決策指標

比單一 2050 elevation 更有用的是：

\Pr(h_t<h_{\rm critical}),

首次跌破 threshold 的分布、shortage duration、expected delivery deficit，以及不同政策的 regret。

十二、模型檢查

elevation、area、volume curve 單調且單位一致；
每月 water balance residual 可追查；
residual ACF 和 diagnostics；
rolling forecast 勝過 naive baseline；
interval coverage 接近 nominal；
temporal aggregation 不混淆 seasonality；
scenario covariates 在 forecast horizon 有明確來源；
wastewater／conservation effect 有 system boundary；
不把 historical association 當作 causal policy effect。

結論

Lake Mead 預測需要把物理與統計放在正確位置。SARIMA 和 UCM 可作短期 benchmark，幫助辨識 trend、seasonality 和 residual structure；長期管理則必須回到庫容守恆、 elevation–storage curve、氣候 inflow ensemble 與 operation rules。

原 2025–2050 點預測和 1,075.98 ft ensemble 因窗口、頻率、區間與權重問題，不應作規劃依據。最可靠的輸出不是一條看似平滑的 30 年線，而是經 rolling validation 的短期預測，加上透明情景下跌破關鍵門檻的概率。

參考資料

COMAP, Inc. (2021). HiMCM 2021 Problem B and accompanying Lake Mead datasets.
U.S. Bureau of Reclamation. Colorado River Basin and Lake Mead operational data and planning studies.
Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
Harvey, A. C. (1989). Forecasting, Structural Time Series Models and the Kalman Filter. Cambridge University Press.
Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.). OTexts.