時間序列分析（進階）：條件異質變異、Kalman 濾波與 Granger 因果

當「平均數會漂移」：為什麼 ARIMA 不夠用，而波動率本身會說話

你已經學過自相關（autocorrelation）、定態（stationarity）與 ARIMA。現在請看一段真實的金融報酬序列：它的平均數幾乎是定態的（每天報酬都在 0 附近抖動），但有些時期波動劇烈、有些時期風平浪靜——而且波動會「成群結隊」地出現。ARIMA 對這種現象束手無策，因為 ARIMA 假設誤差項的變異數是常數（同質變異，homoskedasticity）。

這篇進階文章不再重複「怎麼判斷定態」或「怎麼挑 ARIMA 階數」。我們直接走進三個入門篇刻意略過的深水區：條件異質變異（conditional heteroskedasticity）、狀態空間模型與 Kalman 濾波（state-space & Kalman filter），以及向量自迴歸與 Granger 因果（VAR & Granger causality）。這三者分別回答了三個入門 ARIMA 答不出的問題：變異數會變怎麼辦、隱藏狀態看不見怎麼辦、多條序列互相影響怎麼辦。

一、條件異質變異：當變異數有記憶

ARIMA 把序列拆成「可預測的條件平均 + 白噪音誤差」。問題是它假設誤差 $\varepsilon_t$ 的變異數固定不變。但金融、能源、甚至生理訊號（如 HRV）常出現波動聚集（volatility clustering）：大波動後面常跟著大波動，小波動後面跟著小波動。

ARCH 模型：把變異數也建模

Engle（1982）的 ARCH（Autoregressive Conditional Heteroskedasticity）模型核心想法是：讓今天的條件變異數依賴過去的誤差平方。最簡單的 ARCH(1)：

$$ \varepsilon_t = \sigma_t z_t, \quad z_t \sim \text{i.i.d.}\,\mathcal{N}(0,1) $$

$$ \sigma_t^2 = \omega + \alpha_1 \varepsilon_{t-1}^2 $$

這裡 $z_t$ 是標準化的隨機衝擊，$\sigma_t^2$ 是條件變異數（conditional variance）——它隨時間改變，而且被昨天的誤差平方 $\varepsilon_{t-1}^2$ 推動。昨天衝擊大（$\varepsilon_{t-1}^2$ 大），今天變異數就大，於是今天也容易出現大波動。波動聚集就這樣被機制化地解釋了。

GARCH：給變異數加上自身的慣性

ARCH 要捕捉長記憶的波動，往往需要很高的階數。Bollerslev（1986）的 GARCH(1,1) 加入「變異數對自己的迴歸項」：

$$ \sigma_t^2 = \omega + \alpha_1 \varepsilon_{t-1}^2 + \beta_1 \sigma_{t-1}^2 $$

這個式子之所以在實務上稱霸三十年，是因為它用三個參數就刻畫了極為頑強的波動持續性。其中 $\alpha_1 + \beta_1$ 衡量波動的持續度（persistence）：越接近 1，波動回到長期均值的速度越慢。長期（無條件）變異數為：

$$ \bar{\sigma}^2 = \frac{\omega}{1 - \alpha_1 - \beta_1}, \quad \text{需要 } \alpha_1 + \beta_1 < 1 $$

當 $\alpha_1 + \beta_1 \to 1$，模型逼近 IGARCH（integrated GARCH），波動衝擊不會消散——這正是金融危機期間的特徵。

動手算一下：GARCH(1,1) 的一步波動預測

假設我們估出 $\omega = 0.00001$、$\alpha_1 = 0.08$、$\beta_1 = 0.90$。某天觀察到誤差 $\varepsilon_{t-1} = 0.03$（3% 的異常報酬），且前一天的條件變異數 $\sigma_{t-1}^2 = 0.0004$（對應約 2% 的日波動）。今天的條件變異數：

$$ \sigma_t^2 = 0.00001 + 0.08 \times (0.03)^2 + 0.90 \times 0.0004 $$

$$ = 0.00001 + 0.08 \times 0.0009 + 0.00036 = 0.00001 + 0.000072 + 0.00036 = 0.000442 $$

所以今天的條件標準差 $\sigma_t = \sqrt{0.000442} \approx 0.0210$，即 2.10%——比前一天的 2.00% 明顯上升。這就是「昨天的大波動推高了今天的風險預期」。

我們也可順手檢查持續度：$\alpha_1 + \beta_1 = 0.98$，非常高，波動半衰期約為 $\ln(0.5)/\ln(0.98) \approx 34$ 天，意味著一次衝擊的影響要一個多月才衰退一半。長期變異數 $\bar{\sigma}^2 = 0.00001/(1-0.98) = 0.0005$，對應長期日波動約 $\sqrt{0.0005} \approx 2.24\%$。我們今天的 2.10% 仍在向這個長期水準收斂的路上。

一個常見迷思

很多人以為「波動聚集 = 序列有自相關」。不是。 報酬序列 $\varepsilon_t$ 本身常常幾乎沒有自相關（這符合效率市場），但 $\varepsilon_t^2$（或 $|\varepsilon_t|$）卻有顯著自相關。判斷要不要上 GARCH，正確做法是對殘差平方做 Ljung-Box 檢定或 ARCH-LM 檢定，而不是對殘差本身。

二、狀態空間模型：把看不見的狀態解出來

ARIMA 把一切都寫成觀測值的函數。但很多真實系統有「隱藏狀態」：真實的趨勢、真實的位置、真實的學習能力——我們只能透過充滿雜訊的觀測去推斷它。狀態空間模型（state-space model）正是為此而生。

兩條方程式

狀態空間模型由一對方程式定義。狀態方程（state equation）描述隱藏狀態如何隨時間演化：

$$ x_t = F x_{t-1} + w_t, \quad w_t \sim \mathcal{N}(0, Q) $$

觀測方程（observation equation）描述我們看到的資料如何由狀態產生：

$$ y_t = H x_t + v_t, \quad v_t \sim \mathcal{N}(0, R) $$

$x_t$ 是隱藏狀態（可以是向量），$y_t$ 是觀測，$F$ 是狀態轉移矩陣，$H$ 是觀測矩陣，$Q$ 與 $R$ 分別是過程雜訊與觀測雜訊的共變異數。一個關鍵洞見是：ARIMA 其實是狀態空間模型的特例——任何 ARMA 模型都能寫成這種形式。狀態空間是更廣的語言。

Kalman 濾波：最佳的遞迴更新

給定上面的線性高斯模型，Kalman 濾波（Kalman filter）以遞迴方式給出隱藏狀態的最小均方誤差估計。它分兩步循環：

預測步（predict）——用昨天的估計推測今天：

$$ \hat{x}_{t|t-1} = F \hat{x}_{t-1|t-1}, \quad P_{t|t-1} = F P_{t-1|t-1} F^\top + Q $$

更新步（update）——用今天的觀測修正預測：

$$ K_t = \frac{P_{t|t-1} H^\top}{H P_{t|t-1} H^\top + R}, \quad \hat{x}_{t|t} = \hat{x}_{t|t-1} + K_t (y_t - H \hat{x}_{t|t-1}) $$

$$ P_{t|t} = (I - K_t H) P_{t|t-1} $$

這裡 $K_t$ 是 Kalman 增益（Kalman gain），它決定我們要多相信新觀測。直覺非常漂亮：當觀測雜訊 $R$ 大（資料不可信），$K_t$ 變小，我們就少修正、多依賴模型；當預測不確定性 $P_{t|t-1}$ 大（模型不可信），$K_t$ 變大，我們就多採信新資料。Kalman 增益是模型信念與資料證據之間的動態權衡。

看一個例子：用 Kalman 濾波追蹤一個漂移的均值

假設我們相信某個量測（例如某學習者每日的「真實專注度」）的真實值會緩慢隨機漂移：$x_t = x_{t-1} + w_t$（即 $F=1$），而我們每天的量測 $y_t = x_t + v_t$（即 $H=1$）有雜訊。設過程雜訊 $Q = 0.01$、觀測雜訊 $R = 0.25$。

從 $\hat{x}_0 = 50$、$P_0 = 1$ 開始，今天觀測到 $y_1 = 53$。

預測步：$\hat{x}_{1|0} = 50$，$P_{1|0} = 1 + 0.01 = 1.01$。

計算增益：

$$ K_1 = \frac{1.01}{1.01 + 0.25} = \frac{1.01}{1.26} \approx 0.802 $$

更新步：

$$ \hat{x}_{1|1} = 50 + 0.802 \times (53 - 50) = 50 + 2.405 = 52.405 $$

$$ P_{1|1} = (1 - 0.802) \times 1.01 \approx 0.200 $$

注意兩件事：第一，我們沒有完全跳到 53，而是停在 52.4——因為觀測有雜訊，我們只部分相信它。第二，更新後的不確定性從 1.01 降到 0.20——觀測讓我們對真實狀態更有把握了。隨著資料累積，$P_t$ 會收斂到一個穩態值，$K_t$ 也會穩定下來。

這正是入門 ARIMA 做不到的：ARIMA 給你點預測，但 Kalman 濾波同時給你狀態估計與其不確定性，而且能線上（online）即時更新——這是 GPS 定位、衛星追蹤、感測器融合背後的同一套數學。

三、向量自迴歸與 Granger 因果：多條序列的互動

入門 ARIMA 是單變量的——它只看一條序列的過去。但真實世界裡，序列彼此牽動：利率影響匯率、睡眠影響壓力、二氧化碳濃度影響注意力。向量自迴歸（VAR, Vector Autoregression）把多條序列放進同一個系統。

VAR(p) 的結構

兩變量 VAR(1) 長這樣（以 $y$ 與 $x$ 兩條序列為例）：

$$ \begin{aligned} y_t &= c_1 + \phi_{11} y_{t-1} + \phi_{12} x_{t-1} + \varepsilon_{1t} \\ x_t &= c_2 + \phi_{21} y_{t-1} + \phi_{22} x_{t-1} + \varepsilon_{2t} \end{aligned} $$

每條序列都同時被自己的過去和對方的過去解釋。寫成矩陣形式 $\mathbf{y}_t = \mathbf{c} + \Phi \mathbf{y}_{t-1} + \boldsymbol{\varepsilon}_t$，就是一個向量版的 AR 模型。

Granger 因果：可預測性，不是真因果

如果加入 $x$ 的過去值能顯著改善對 $y$ 的預測（相較於只用 $y$ 自己的過去），我們就說 $x$ Granger 導致（Granger-causes） $y$。檢定方式是對上面 $y_t$ 方程式中所有 $x$ 的落後項係數做聯合 F 檢定：

$$ H_0: \phi_{12}^{(1)} = \phi_{12}^{(2)} = \dots = \phi_{12}^{(p)} = 0 $$

若拒絕 $H_0$，則 $x$ 對 $y$ 有 Granger 因果。

這裡有一個必須講清楚的迷思：Granger 因果不是哲學意義上的真因果。它只是「時間上的可預測性增益」。$x$ Granger 導致 $y$，可能是因為兩者都被一個未觀測的共同因子驅動，也可能是純粹的領先—落後關係。一個經典反例：公雞啼叫 Granger 導致日出（啼叫總在日出前），但顯然不是公雞讓太陽升起。Granger 因果是必要的線索，不是充分的證明。

動手算一下：用 RSS 計算 Granger 檢定統計量

假設我們要檢定「$x$ 是否 Granger 導致 $y$」，落後階數 $p=2$，樣本 $T=100$。

• 受限模型（只用 $y$ 的落後項）的殘差平方和 $\text{RSS}_R = 52.0$
• 非受限模型（加入 $x$ 的落後項）的殘差平方和 $\text{RSS}_U = 44.0$

加入 $x$ 的 2 個落後項後，非受限模型估了 5 個參數（截距 + 2 個 $y$ 落後 + 2 個 $x$ 落後）。F 統計量：

$$ F = \frac{(\text{RSS}_R - \text{RSS}_U)/q}{\text{RSS}_U/(T-k)} = \frac{(52.0 - 44.0)/2}{44.0/(100-5)} $$

$$ = \frac{8.0/2}{44.0/95} = \frac{4.0}{0.4632} \approx 8.64 $$

在分子自由度 $q=2$、分母自由度 $95$ 下，$F=8.64$ 遠超過 5% 臨界值（約 3.09），因此拒絕 $H_0$——資料支持 $x$ 對 $y$ 有 Granger 因果（可預測性增益顯著）。但記住：這只說明 $x$ 的過去含有預測 $y$ 的資訊，不直接證明機制性的因果。

致命陷阱：偽迴歸

對多條序列做 VAR 前，定態檢查比單變量時更重要。若兩條都是非定態（含單位根）的隨機漫步，直接跑 VAR 或迴歸會出現偽迴歸（spurious regression）：$R^2$ 很高、t 值很顯著，但全是假的。正確做法是先做單位根檢定（如 ADF），若都非定態，再檢查它們是否共整合（cointegration）——若存在某個線性組合是定態的，就改用 VECM（向量誤差修正模型）；若無共整合，則對差分後的序列建 VAR。

重點回顧

• ARIMA 假設變異數固定；面對波動聚集要用 ARCH/GARCH 把條件變異數也建模。判斷依據是檢定殘差平方的自相關，不是殘差本身。
• GARCH(1,1) 的 $\alpha_1 + \beta_1$ 衡量波動持續度，越接近 1 衰退越慢；$\geq 1$ 時長期變異數不存在（IGARCH）。
• 狀態空間模型把隱藏狀態與觀測分離，ARIMA 是它的特例。Kalman 濾波以預測—更新遞迴給出狀態的最佳估計與其不確定性。
• Kalman 增益是模型信念與資料證據的動態權衡：觀測越可信，增益越大，越採信新資料。
• VAR 處理多序列互動，Granger 因果是「可預測性」不是「真因果」；多序列建模前務必檢查定態與共整合，否則落入偽迴歸。

深入探討（研究所視角）

GARCH 家族的非對稱延伸。 標準 GARCH 對「好消息」與「壞消息」一視同仁——正負誤差平方後相同。但金融市場有槓桿效應（leverage effect）：壞消息（價格下跌）對波動的衝擊大於同等幅度的好消息。EGARCH（Nelson, 1991）對 $\ln \sigma_t^2$ 建模並引入 $z_{t-1}$ 的符號項以捕捉非對稱；GJR-GARCH 則加入一個指示函數項 $\gamma\, \mathbb{1}(\varepsilon_{t-1}<0)\,\varepsilon_{t-1}^2$。此外，估計通常採用準最大概似（QMLE），即使 $z_t$ 並非真的常態，QMLE 仍給出一致估計；但實務上常假設 $z_t$ 為 Student's $t$ 分布以容納金融報酬的厚尾（heavy tails）。

Kalman 濾波與貝氏濾波的統一。 Kalman 濾波其實是貝氏遞迴濾波在線性高斯假設下的封閉解——預測步是先驗的時間傳播，更新步就是貝氏定理。一旦模型非線性或非高斯，封閉解不存在，要改用擴展 Kalman 濾波（EKF，對非線性做一階泰勒線性化）、無跡 Kalman 濾波（UKF，用 sigma points 傳播分布的前兩階動差），或粒子濾波（particle filter）——用蒙地卡羅樣本近似整個後驗分布。粒子濾波讓我們能處理任意非線性、非高斯的狀態空間，代價是計算量大。值得一提的是，狀態空間模型的參數（$F, H, Q, R$）本身可用 EM 演算法估計，其中 E 步用 Kalman 平滑器（RTS smoother）算出狀態的後驗期望，M 步再更新參數。

共整合的深層結構。 Engle 與 Granger 因「共整合」與「ARCH」分享 2003 年諾貝爾經濟學獎，這並非巧合——兩者都在處理「定態假設失效」的不同面向。共整合的核心是：若 $d$ 條 $I(1)$ 序列存在 $r$ 個線性獨立的定態組合，則它們被 $r$ 個長期均衡關係綁住，短期可偏離但會被「誤差修正項」拉回。Johansen 檢定透過 VAR 的係數矩陣 $\Pi$ 的秩來判定共整合向量的個數：$\text{rank}(\Pi) = r$。當 $r$ 介於 0 與 $d$ 之間，系統就有 VECM 表示。這把「短期動態」與「長期均衡」優雅地放進同一個框架——正是把入門時間序列推向計量經濟學前沿的關鍵橋樑。