ARK · 金融機器學習CHAPTER 07 / 12

CHAPTER 07 / 12 · PART 2 · 時序學習

隱馬可夫模型與粒子濾波

Probabilistic Sequence Modeling

市場狀態看不見、只看得見報酬——HMM、Kalman 與粒子濾波教你從觀測倒推隱狀態。

§01學習重點

§02課程內容

一、看不見的狀態，看得見的報酬

打開任何一張大盤日 K 線圖，你看到的是一串報酬：漲、跌、漲漲、跌。你看不到的，是市場「處在什麼狀態」。多頭市場、空頭市場、橫盤整理——這些詞每天出現在財經新聞裡，但沒有任何一個交易所會掛牌公告「本日起進入熊市」。狀態是事後追認的、爭論不休的、本質上觀測不到的。你手上有的，永遠只是狀態投影出來的影子：一天一筆的報酬。

這個「影子與本體分離」的結構，值得用一個正式的框架來裝。設一條看不見的狀態序列 \(s_1, s_2, \dots, s_T\)（每天市場處於牛市或熊市），它自己按某種規律隨時間演化；再設一條看得見的觀測序列 \(y_1, y_2, \dots, y_T\)（每天的報酬方向），每一天的觀測由當天的狀態決定機率分布。狀態驅動觀測、觀測透露狀態，但兩者永遠不相等——這就是隱馬可夫模型（Hidden Markov Model，HMM）的世界觀，也是更廣的狀態空間模型（state-space model）家族的共同骨架。

把它跟你已經會的東西接上。第 2 章的高斯混合模型（GMM）已經有「隱變數」了：每個樣本背後有一個看不見的成分歸屬，先抽成分、再從該成分的分布抽觀測。但 GMM 的隱變數是獨立同分布的——今天抽到哪個成分，跟昨天抽到哪個毫無關係。金融資料顯然不是這樣：熊市不會一天結束，波動高的日子傾向連著來（第 6 章講 GARCH 時你已經見過波動叢聚）。HMM 對混合模型只做了一個修改，卻是關鍵的修改：讓隱變數自己形成一條馬可夫鏈——今天的狀態依賴昨天的狀態。一句話總結：

HMM＝混合模型＋時間記憶。

同時它也補上了第 6 章工具箱的一個缺口。自迴歸、移動平均這些模型直接對「觀測本身」建方程，沒有隱藏層；GARCH 雖然有一條看不見的波動度序列，但那條序列是由過去觀測確定性地算出來的，不是自己帶隨機性的過程。HMM 與狀態空間模型讓隱狀態擁有自己的隨機動態，於是「市場換了一個習性」這件事可以被建模成「隱狀態跳到另一格」，而不必假設整段歷史服從同一組參數——這正是它們面對非平穩資料時比前一章的頻率派時序模型多出來的本錢。

比喻： 你和室友隔著一道牆。你聽得見聲音——鍵盤聲、翻書聲、鼾聲——但看不見他在做什麼。

你心裡其實一直在跑一個 HMM：他的「狀態」（趕報告、耍廢、睡著）看不見，你聽到的聲音是狀態

「發射」出來的觀測——趕報告時多半是鍵盤聲，但偶爾也翻書；睡著時多半安靜，但偶爾打鼾。而且

狀態有黏性：睡著的人下一分鐘大概率還睡著，不會每分鐘重新擲骰子決定人生。聽到連續十分鐘

的鍵盤聲，你敢打賭他在趕報告——你剛剛就做了一次隱狀態推斷。

這一章的路線圖：先把 HMM 的三要素與三個經典問題講清楚（離散狀態），再把狀態從離散推廣到連續——線性高斯世界的 Kalman 濾波，以及非線性非高斯世界的粒子濾波——最後落地到金融應用與濾波器的校準。

二、HMM 三要素：初始、轉移、發射

一個 HMM 由三組參數完全決定。假設隱狀態有 \(K\) 種（牛、熊就是 \(K=2\)），觀測有 \(N\) 種符號（漲、跌就是 \(N=2\)）：

第一，初始分布 \(\pi\)：一個長度 \(K\) 的機率向量，\(\pi_k\) 是第一天處於狀態 \(k\) 的機率。故事的開頭從哪裡抽。

第二，轉移矩陣 \(A\)：一個 \(K\times K\) 矩陣，\(A_{ij}\) 是「今天在狀態 \(i\)、明天跳到狀態 \(j\)」的機率，每一列加總為 1。它就是狀態的「黏性表」——對角線元素大，代表狀態持久；對角線小，代表市場性格善變。這個矩陣正是 HMM 與混合模型的全部差別：把 \(A\) 的每一列都改成同一個向量（明天去哪跟今天無關），HMM 立刻退化回第 2 章的混合模型。

第三，發射分布 \(B\)：狀態決定觀測的機率。離散觀測時是一個 \(K\times N\) 矩陣，\(B_{kn}\) 是「狀態 \(k\) 之下觀測到符號 \(n\)」的機率；連續觀測時每個狀態配一個密度（例如牛市配均值略正的高斯、熊市配均值略負且更胖尾的高斯）。狀態不直接露臉，它透過發射分布「說話」。

三要素到齊，整條序列的聯合機率就可以一步一步拆開。狀態鏈是馬可夫的（明天只看今天），觀測只看當天狀態，所以隱狀態序列 \(\mathbf{s}\) 與觀測序列 \(\mathbf{y}\) 的聯合機率是一條規整的連乘：

$$ p(\mathbf{s}, \mathbf{y}) = \pi_{s_1}\, B_{s_1}(y_1) \prod_{t=2}^{T} A_{s_{t-1} s_t}\, B_{s_t}(y_t) $$

逐項拆解：\(\pi_{s_1}\) 是開局落在狀態 \(s_1\) 的機率；\(B_{s_1}(y_1)\) 是該狀態發射出第一筆觀測的機率；之後每往前走一天，就乘上一次「狀態轉移」\(A_{s_{t-1}s_t}\) 與一次「觀測發射」\(B_{s_t}(y_t)\)。整條式子讀起來就是故事本身：抽開局、發射、轉移、發射、轉移、發射……任何一條給定的狀態路徑配上觀測序列，其機率都能沿著這條連乘直接算出來。

本章通篇使用同一個自創的玩具市場，把數字先立好：狀態 0＝牛市、狀態 1＝熊市；觀測 0＝當日漲、1＝當日跌。

$$ \pi = (0.75,\; 0.25),\qquad A = \begin{pmatrix} 0.85 & 0.15 \\ 0.10 & 0.90 \end{pmatrix},\qquad B = \begin{pmatrix} 0.65 & 0.35 \\ 0.40 & 0.60 \end{pmatrix} $$

逐項拆解：開局七成五機率是牛市；牛市有 85% 機率延續、熊市更黏（90% 延續）——這對「熊市又長又難熬」的市場直覺點了頭；牛市裡單日上漲機率 65%，熊市裡下跌機率 60%。注意發射分布故意設計得很不乾脆：牛市也有 35% 的日子在跌。單看一天的漲跌幾乎說明不了什麼，必須把整條序列的證據累積起來——這正是接下來三個經典問題的舞台。

三、三個經典問題與前向算法

拿到一個 HMM，你會想問的事幾乎都落在三個經典問題裡：

評估（evaluation）——給定參數與觀測序列，這條序列出現的機率 \(p(\mathbf{y})\) 是多少？這是模型比較的基礎：兩組參數誰更能解釋眼前的資料，比的就是這個數。

解碼（decoding）——給定參數與觀測序列，最有可能的隱狀態路徑是哪一條？「過去這兩個月，市場到底哪幾天在熊市？」問的就是解碼。

學習（learning）——只給觀測序列，參數 \(\pi, A, B\) 本身要從資料裡估出來。

先攻評估。直覺的算法是暴力枚舉：把所有可能的狀態路徑列出來，每條用聯合機率公式算一遍再加總。問題是路徑有 \(K^T\) 條——兩個狀態、250 個交易日就是 \(2^{250}\) 條，宇宙的原子數都不夠用。出路是前向算法（forward algorithm）。定義前向變數：

$$ \alpha_t(k) = p(y_1, \dots, y_t,\; s_t = k) $$

逐項拆解：\(\alpha_t(k)\) 是「前 \(t\) 天的觀測全都發生、且第 \(t\) 天恰好處於狀態 \(k\)」的聯合機率。它把「到目前為止的所有路徑證據」壓縮成 \(K\) 個數字。遞迴關係一行寫完：

$$ \alpha_1(k) = \pi_k B_k(y_1), \qquad \alpha_{t}(k) = \Big(\sum_{i=1}^{K} \alpha_{t-1}(i)\, A_{ik}\Big)\, B_k(y_t) $$

逐項拆解：第一天，開局機率乘上發射機率。之後每一天分兩步——轉移：把昨天各狀態的前向變數沿轉移矩陣匯流到今天的狀態 \(k\)（括號裡的加總）；發射：乘上狀態 \(k\) 產生今日觀測的機率。走到最後一天，把 \(K\) 個前向變數加總就是 \(p(\mathbf{y})\)。計算量從 \(K^T\) 條路徑塌縮成 \(T\times K^2\) 次乘加——它能這麼省，靠的是馬可夫性質：所有通往「今天狀態 \(k\)」的路徑，未來長得一模一樣，所以可以先加總、再前進，不必各走各的。

用玩具市場手算兩步。觀測序列取「漲、跌、跌、漲、漲」。第一天觀測到漲：\(\alpha_1(\text{牛}) = 0.75\times 0.65 = 0.4875\)，\(\alpha_1(\text{熊}) = 0.25\times 0.40 = 0.10\)。第二天觀測到跌，先轉移再發射：\(\alpha_2(\text{牛}) = (0.4875\times 0.85 + 0.10\times 0.10)\times 0.35 = 0.14853125\)；\(\alpha_2(\text{熊}) = (0.4875\times 0.15 + 0.10\times 0.90)\times 0.60 = 0.097875\)。剩下交給程式，並用暴力枚舉對答案：

PYTHON
import numpy as np

# 兩狀態 HMM：0=牛市、1=熊市；觀測：0=漲、1=跌
pi = np.array([0.75, 0.25])                  # 初始分布
A  = np.array([[0.85, 0.15],                 # 轉移矩陣：列=今天、行=明天
               [0.10, 0.90]])
B  = np.array([[0.65, 0.35],                 # 發射矩陣：列=狀態、行=觀測
               [0.40, 0.60]])
obs = [0, 1, 1, 0, 0]                        # 觀測序列：漲 跌 跌 漲 漲

alpha = pi * B[:, obs[0]]                    # 前向變數 α₁
for t in range(1, len(obs)):
    alpha = (alpha @ A) * B[:, obs[t]]       # 遞迴：先轉移、再發射
print("前向算法 P(y) =", alpha.sum())

# 暴力法：枚舉全部 2^5 = 32 條隱狀態路徑，逐條算聯合機率後加總
from itertools import product
total = 0.0
for path in product([0, 1], repeat=len(obs)):
    p = pi[path[0]] * B[path[0], obs[0]]
    for t in range(1, len(obs)):
        p *= A[path[t-1], path[t]] * B[path[t], obs[t]]
    total += p
print("暴力枚舉 P(y) =", total)
# 輸出：
# 前向算法 P(y) = 0.03009501678466797
# 暴力枚舉 P(y) = 0.03009501678466796

兩個數字到小數第十六位才分家——那是浮點加法順序造成的精度殘差，數學上完全相同。五天的序列暴力法還跑得動；換成一年的資料，前向算法照樣一眨眼，暴力法則已經不在可行宇宙裡。

順帶把學習問題的直覺交代掉。若隱狀態路徑有標籤（哪天牛哪天熊都知道），估參數就是數頻率的統計作業。麻煩在於狀態看不見——這正是第 2 章 EM 演算法登場過的局面。Baum-Welch 演算法就是 EM 的時序版：E 步用前向變數搭配一個對稱的「後向變數」（從序列尾端往回累積證據），算出每一天處於各狀態的後驗機率、以及相鄰兩天各種轉移的後驗機率；M 步拿這些「軟歸屬」當加權樣本，重新估 \(\pi, A, B\)；反覆迭代到收斂。跟 GMM 的 EM 一樣，它保證概似單調上升、也一樣可能停在局部最優，實務上要多組初始值起跑。長推導本課不展開，機制與完整式子見原書 pp.224–226。

四、Viterbi：動態規劃找最佳路徑

解碼問題要的不是機率總和，而是單獨一條最好的路徑：在所有 \(K^T\) 條隱狀態路徑中，找出與觀測聯合機率最大的那條。暴力枚舉同樣不可行，而救兵幾乎是免費的——把前向算法裡的「加總」換成「取最大值」，就得到 Viterbi 演算法。定義：

$$ V_t(k) = \max_{s_1,\dots,s_{t-1}} p(s_1,\dots,s_{t-1},\, s_t = k,\; y_1,\dots,y_t) $$

逐項拆解：\(V_t(k)\) 是「所有在第 \(t\) 天抵達狀態 \(k\) 的路徑中，最好那條」的聯合機率。遞迴與前向算法平行：

$$ V_t(k) = \Big(\max_{i}\; V_{t-1}(i)\, A_{ik}\Big)\, B_k(y_t) $$

逐項拆解：要在今天最優地抵達狀態 \(k\)，只需要知道「昨天各狀態的最優抵達分數」乘上轉移機率後誰最大——因為最佳路徑的任何前段，必然也是通往該中繼點的最佳路徑（不然把前段換掉，整條會更好，矛盾）。這個性質讓我們每天只需替每個狀態記兩件事：最優分數 \(V_t(k)\)，以及它從哪個前狀態來（回溯指標）。走到最後一天，挑分數最高的終點狀態，沿回溯指標一路倒著走回第一天，整條最佳路徑就出土了。

比喻： 規劃台北到高雄的火車轉乘。你不必列出所有可能的轉乘組合——只要在每個中繼站

（台中、嘉義、台南）記住「到這一站最省時的走法是哪條、它從哪一站接過來」就夠了。因為

不管你最後怎麼走到高雄，路過台南時，你在台南以前的走法一定是「到台南的最快走法」。

到了終點站，沿著各站記下的「從哪來」倒著讀，最佳路線自動浮現。Viterbi 的回溯指標，

就是每一站牆上那張「本站最快抵達方式」的小抄。

跑一個會出戲的例子。觀測序列取「漲、漲、跌、跌、跌、漲」——前兩天像牛市、中間三天連跌、最後一天又漲。直覺會說中段是熊市吧?看程式怎麼判：

PYTHON
import numpy as np

pi = np.array([0.75, 0.25])
A  = np.array([[0.85, 0.15], [0.10, 0.90]])
B  = np.array([[0.65, 0.35], [0.40, 0.60]])
obs = [0, 0, 1, 1, 1, 0]                     # 漲 漲 跌 跌 跌 漲
K, T = 2, len(obs)

V = np.zeros((T, K)); back = np.zeros((T, K), dtype=int)
V[0] = pi * B[:, obs[0]]
for t in range(1, T):
    for k in range(K):
        cand = V[t-1] * A[:, k]              # 從各前狀態走到 k 的候選分數
        back[t, k] = np.argmax(cand)         # 記下最佳前手
        V[t, k] = cand.max() * B[k, obs[t]]
    print(f"t={t+1}: V(牛)={V[t,0]:.6f}  V(熊)={V[t,1]:.6f}")

path = [int(np.argmax(V[-1]))]               # 終點取分數最高的狀態
for t in range(T - 1, 0, -1):
    path.append(int(back[t, path[-1]]))      # 沿 back 指標回溯
path.reverse()
print("Viterbi 最佳路徑:", ["牛" if s == 0 else "熊" for s in path])

# 對照組：每個時點各自取「濾波機率最高」的狀態
alpha = pi * B[:, obs[0]]
greedy = [int(np.argmax(alpha))]
for t in range(1, T):
    alpha = (alpha @ A) * B[:, obs[t]]
    alpha /= alpha.sum()
    greedy.append(int(np.argmax(alpha)))
print("逐時點濾波 argmax:", ["牛" if s == 0 else "熊" for s in greedy])
# 輸出：
# t=2: V(牛)=0.269344  V(熊)=0.036000
# t=3: V(牛)=0.080130  V(熊)=0.024241
# t=4: V(牛)=0.023839  V(熊)=0.013090
# t=5: V(牛)=0.007092  V(熊)=0.007069
# t=6: V(牛)=0.003918  V(熊)=0.002545
# Viterbi 最佳路徑: ['牛', '牛', '牛', '牛', '牛', '牛']
# 逐時點濾波 argmax: ['牛', '牛', '牛', '熊', '熊', '熊']

兩個答案吵起來了。Viterbi 說整整六天都是牛市：連跌三天雖然可疑，但要解釋成熊市，路徑得付出「牛轉熊」的轉移代價（0.15 的過路費），還得在最後一天用熊市解釋上漲（0.40 對 0.65 的發射劣勢）——加加減減之後，「牛市裡的一段回檔」整條路徑的機率反而更高（注意 t=5 兩邊分數 0.007092 對 0.007069，只差千分之三，最後一天的上漲才讓天平倒向牛）。而逐時點的濾波機率在第 4 天起已經偏向熊——它每一天只回答「此刻最像什麼」，不需要對整條路徑的一致性負責。

兩者沒有誰對誰錯，它們回答不同的問題。濾波機率的逐點 argmax 甚至可能拼出一條「轉移機率為零、根本不可能發生」的狀態序列；Viterbi 保證整條路徑合法且整體最優，卻可能在個別時點上與「當下最可能的狀態」相左。做市場狀態報告時引用哪一個，取決於你要的是「逐日的即時判讀」還是「事後一致的歷史分段」——把這個區別講清楚，是用 HMM 做狀態分析的第一堂職業素養課。

五、連續隱狀態：狀態空間模型與 Kalman 濾波

牛熊二分是個好起點，但很多金融隱狀態天生是連續的：資產的公允價值、潛在的波動度水準、殖利率曲線的潛在因子。把 HMM 的離散狀態換成連續向量，就得到一般的狀態空間模型，寫成兩條方程：

$$ s_t = F(s_{t-1}) + \varepsilon_t, \qquad y_t = G(s_t) + \xi_t $$

逐項拆解：第一條是狀態方程——隱狀態沿某個函數 \(F\) 演化，加上過程噪聲 \(\varepsilon_t\)（世界本身的隨機性）；第二條是觀測方程——你看到的 \(y_t\) 是隱狀態經過函數 \(G\) 的投影，再疊上觀測噪聲 \(\xi_t\)（量測的隨機性）。HMM 是它的離散特例：\(F\) 換成轉移矩陣抽籤、\(G\) 換成發射分布抽籤。

當 \(F\)、\(G\) 都是線性、兩個噪聲都是高斯，這個模型有一個著名的精確解法：Kalman 濾波。它維護的不是一團樣本，而是隱狀態的完整後驗分布——因為高斯分布經過線性變換還是高斯、兩個高斯的貝葉斯更新還是高斯，所以後驗永遠是高斯，用兩個數（均值與變異數）就完整刻畫。每個時間步走一個兩拍迴圈：

預測：把信念沿狀態方程往前推。均值跟著動態走，變異數膨脹（往前猜一步，把握變小）。修正：新觀測進來，把預測值往觀測方向拉，拉多少由 Kalman 增益決定。一維情形的增益一目瞭然：

$$ K_t = \frac{P_t}{P_t + r^2} $$

逐項拆解：\(P_t\) 是預測的不確定性（你對自己猜測的不放心程度），\(r^2\) 是觀測噪聲的變異數（儀器有多不可靠）。兩者的比例決定信任分配——自己很不確定、儀器很準時，增益趨近 1（幾乎全信觀測）；自己很有把握、儀器很吵時，增益趨近 0（幾乎不理觀測）。Kalman 濾波的全部智慧就在這個分數裡：它是一台自動調參的信任分配機。

比喻： 山路開車的導航。進隧道，GPS 斷了，導航靠「剛才的位置＋車速」推算你到哪

（預測），而且推得越久、藍點的光暈畫得越大——它誠實地承認越來越不確定。出隧道，

衛星訊號回來，藍點被往量到的位置拉回去（修正），光暈瞬間收小。拉多少？訊號滿格就

大膽拉，訊號只有一格就只拉一點點——這個「該信推算幾分、信衛星幾分」的權重，就是

Kalman 增益。你的手機每一秒都在跑這個預測、修正迴圈。

拿一個微結構味道的玩具驗證：看不見的「公允價值偏離」慢慢隨機遊走（持續性 0.98），你只看得到疊了買賣價差噪聲的成交價。比較「直接拿成交價當估計」與「Kalman 濾波」誰貼得近：

PYTHON
import numpy as np

rng = np.random.default_rng(21)
T = 250
phi, q, r = 0.98, 0.10, 0.60          # 狀態持續性、狀態噪聲sd、觀測噪聲sd

x = np.zeros(T); y = np.zeros(T)      # x=看不見的公允價值偏離, y=成交價偏離
y[0] = x[0] + rng.normal(0, r)
for t in range(1, T):
    x[t] = phi * x[t-1] + rng.normal(0, q)
    y[t] = x[t] + rng.normal(0, r)

m, P = 0.0, 1.0                       # 狀態的均值與變異數（我們的「信念」）
est = np.zeros(T)
for t in range(T):
    if t > 0:
        m = phi * m                   # 預測：信念沿動態方程前推
        P = phi**2 * P + q**2         #        不確定性隨之膨脹
    K = P / (P + r**2)                # Kalman 增益：該信新觀測幾分
    m = m + K * (y[t] - m)            # 修正：往觀測方向挪 K 的比例
    P = (1 - K) * P                   #        不確定性收縮
    est[t] = m
print(f"直接拿成交價當估計 RMSE: {np.sqrt(np.mean((y - x)**2)):.4f}")
print(f"Kalman 濾波估計   RMSE: {np.sqrt(np.mean((est - x)**2)):.4f}")
print(f"收斂後的 Kalman 增益: {K:.4f}")
# 輸出：
# 直接拿成交價當估計 RMSE: 0.5779
# Kalman 濾波估計   RMSE: 0.1822
# 收斂後的 Kalman 增益: 0.1386

誤差從 0.578 壓到 0.182——三倍多的改善，來源沒有任何魔法：模型知道真值移動得很慢（狀態噪聲 0.10 遠小於觀測噪聲 0.60），所以增益收斂在 0.139，每天只讓新成交價把估計拉動 14% 左右，其餘靠歷史累積的信念撐著。濾波的本質是用「狀態動態的先驗知識」去換「觀測噪聲的壓制」。

還有一個值得記下的對照：第 8 章會遇到循環神經網路（RNN），它的隱狀態更新式長得跟狀態空間模型神似，但有兩個本質差異——RNN 的狀態轉移是非線性的、且同時吃前一狀態與當期觀測；以及它的狀態方程沒有噪聲項，隱狀態是確定性計算的結果。可以粗略地說：RNN 是把「濾波器」直接用資料學出來的確定性親戚。

六、粒子濾波：當世界不再線性高斯

Kalman 濾波的精確性建立在一個嬌貴的前提上：線性動態、高斯噪聲。前提一破——狀態方程有指數、觀測分布是胖尾、後驗長出兩個峰——「兩個數刻畫整個後驗」就不再成立。這時需要一個不挑食的近似框架：粒子濾波（particle filtering）。

想法樸素到令人意外：既然後驗分布寫不出解析式，那就用一群樣本點代表它。維護 \(M\) 顆「粒子」，每顆是隱狀態的一個候選值、帶一個權重；粒子雲的加權直方圖就是後驗分布的近似，粒子越多、近似越好，而且對分布形狀零假設——雙峰、偏斜、胖尾統統照單全收。最常用的序貫重要性重抽樣（Sequential Importance Resampling，SIR）演算法，每個時間步走三拍：

第一拍，抽樣：每顆粒子各自沿狀態方程往前走一步（把過程噪聲抽出來加上去）。粒子雲代表「預測分布」。第二拍，加權：新觀測進來，每顆粒子問一句「如果隱狀態真是我，看到這筆觀測的概似多高?」概似高的加權、概似低的降權，權重歸一化。加權後的粒子雲就代表「後驗分布」，此刻要估什麼（後驗均值、分位數）都從它算。第三拍，重抽樣：按權重比例把粒子重新抽一輪——高權重的粒子被複製多份，低權重的被淘汰，權重歸零重設為均等。

第三拍為什麼必要?因為沒有它，演算法會得一種職業病：粒子退化（degeneracy）。權重一輪一輪連乘下去，運氣稍差的粒子權重指數式萎縮，幾十步之後幾乎全部權重集中在一兩顆粒子上——名義上有兩千顆粒子，實際上只剩一顆在幹活，後驗近似徹底失真。診斷指標是有效樣本數：

$$ \mathrm{ESS} = \frac{1}{\sum_{i=1}^{M} (w^{(i)})^2} $$

逐項拆解：\(w^{(i)}\) 是歸一化後的粒子權重。權重完全均等時 \(\mathrm{ESS}=M\)（人人有貢獻）；權重全壓在一顆粒子時 \(\mathrm{ESS}=1\)（只剩獨腳戲）。實務上盯著 ESS 走勢，跌破門檻（常用 \(M/2\)）就重抽樣。重抽樣不是免費午餐——它引入額外的蒙地卡羅噪聲、且複製出的粒子彼此相同（多樣性下降）——但相比放任退化，這個代價非常划算。

比喻： 山區搜救。失蹤者的位置是隱狀態，指揮官撒出兩千名志工（粒子）沿「他可能移動

的方向」散開（抽樣）。新線索進來——溪邊發現一只鞋——每個搜索區按「線索與此區的吻合度」

打分（加權），然後指揮官果斷調度：把人力從低分區抽走、加倍投入高分區（重抽樣）。若從不

重新調度，兩天後多數志工還在最初隨機分配的山頭空轉，等於全隊只剩幾個人真正在搜——

這就是粒子退化。重抽樣，就是那位不斷把人力壓到證據最熱處的指揮官。

粒子濾波的代價同樣要攤開講：計算量是 Kalman 的百倍千倍（每步 \(M\) 次狀態演化與概似評估）；高維狀態空間裡粒子數需求爆炸；而在線性高斯世界裡它只是 Kalman 濾波的吵雜近似——能用 Kalman 就用 Kalman，粒子濾波是留給非線性非高斯場合的重武器。作業三會讓你親眼看到粒子濾波隨 \(M\) 增大收斂到 Kalman 精確解的過程。

七、金融應用：隨機波動度與市場狀態切換

粒子濾波在金融裡最經典的用武之地是隨機波動度模型（stochastic volatility，SV）。第 6 章的 GARCH 把波動度寫成過去觀測的確定性函數；SV 模型更進一步，讓波動度自己是一條隨機隱過程。標準寫法用對數變異數 \(x_t\) 當隱狀態：

$$ y_t = e^{x_t / 2}\, \epsilon_t, \qquad x_{t+1} = \mu(1-\phi) + \phi\, x_t + \sigma_v\, \eta_t $$

逐項拆解：第一條是觀測方程——日報酬 \(y_t\) 等於「當日波動度」\(e^{x_t/2}\)（對數變異數的一半取指數，即標準差）乘上標準高斯噪聲 \(\epsilon_t\)；第二條是狀態方程——對數變異數走一條 AR(1)，\(\mu\) 是長期水準、\(\phi\) 是持續性（接近 1 表示高波動期綿延）、\(\sigma_v\) 控制波動度自身的抖動。這個模型對 Kalman 是硬骨頭：觀測方程裡隱狀態藏在指數裡、與噪聲相乘，非線性得徹底。文獻中的完整版本還會讓兩條噪聲負相關以捕捉槓桿效應（跌市推升波動），並疊加跳躍項——原書以 SVLJ 模型為例完整推導了對應的粒子濾波（pp.230–231），本課用最素的版本把機制打通。

bootstrap 粒子濾波（用狀態方程本身當提議分布的 SIR）在這裡三拍照走：粒子沿 AR(1) 前進；概似是「以 \(e^{x/2}\) 為標準差的高斯在觀測值 \(y_t\) 的密度」；然後重抽樣。跑一段 500 天的模擬——參數設在日報酬的真實量級（\(\mu=-8.4\) 對應長期波動度約 1.5%，\(\phi=0.97\)，\(\sigma_v=0.25\)），對照組是天真估計「拿 \(\ln y_t^2\) 直接當對數變異數」：

PYTHON
import numpy as np

rng = np.random.default_rng(77)
T = 500
mu, phi_v, sig_v = -8.4, 0.97, 0.25   # 對數變異數的長期均值、持續性、波動

xt = np.zeros(T); ret = np.zeros(T)   # xt=隱藏的對數變異數, ret=日報酬
xt[0] = mu + rng.normal(0, sig_v / np.sqrt(1 - phi_v**2))
ret[0] = np.exp(xt[0] / 2) * rng.normal()
for t in range(1, T):
    xt[t] = mu * (1 - phi_v) + phi_v * xt[t-1] + rng.normal(0, sig_v)
    ret[t] = np.exp(xt[t] / 2) * rng.normal()

M = 2000
parts = mu + rng.normal(0, sig_v / np.sqrt(1 - phi_v**2), M)   # 初始粒子雲
est = np.zeros(T); ess_min = M
for t in range(T):
    if t > 0:                          # (1) 抽樣：每顆粒子沿狀態方程前進
        parts = mu * (1 - phi_v) + phi_v * parts + rng.normal(0, sig_v, M)
    var = np.exp(parts)                # (2) 加權：報酬在該粒子下的概似
    logw = -0.5 * (np.log(2 * np.pi * var) + ret[t]**2 / var)
    w = np.exp(logw - logw.max()); w /= w.sum()
    est[t] = np.sum(w * parts)         #     後驗均值 = 粒子的加權平均
    ess_min = min(ess_min, 1.0 / np.sum(w**2))
    idx = rng.choice(M, size=M, p=w)   # (3) 重抽樣：複製高權重、淘汰低權重
    parts = parts[idx]
print(f"粒子濾波 RMSE(對數變異數): {np.sqrt(np.mean((est - xt)**2)):.4f}")
naive = np.log(ret**2)                 # 天真估計：ln(單日報酬平方)
print(f"天真估計 RMSE(對數變異數): {np.sqrt(np.mean((naive - xt)**2)):.4f}")
print(f"全程最低有效樣本數 ESS: {ess_min:.1f} / {M}")
# 輸出：
# 粒子濾波 RMSE(對數變異數): 0.5591
# 天真估計 RMSE(對數變異數): 2.7917
# 全程最低有效樣本數 ESS: 128.1 / 2000

粒子濾波把追蹤誤差壓到天真估計的五分之一。差距的來源看一眼就懂：單日報酬平方是波動度的極度嘈雜的代理（\(\ln \epsilon_t^2\) 這個噪聲項的標準差超過 2，還帶著約 \(-1.27\) 的系統性偏移），天真估計整天被單日運氣甩來甩去；粒子濾波則透過狀態方程知道「波動度是黏的」，把 500 天的證據沿時間軸整合起來。與第五節的 Kalman 實驗遙相呼應：濾波的價值永遠來自動態先驗對觀測噪聲的壓制。另外注意 ESS 全程最低掉到 128——即使每步重抽樣，遇到出乎意料的觀測（平靜期突然一根大棒子）時權重仍會瞬間集中，這正是實務上要監控 ESS 的原因。

實務代碼寫作有兩個不起眼但會咬人的細節，上面的程式都示範了：概似在對數空間計算、減去最大值後再取指數（不然 \(e^{x_t}\) 量級的變異數會讓權重下溢成零）；重抽樣用 rng.choice 按權重抽 \(M\) 次，這就是原書所稱的多項式重抽樣（multinomial resampling）。

離散狀態這邊的金融應用同樣經典：市場狀態切換偵測。把第二節的玩具模型換上真實日報酬（發射分布改成兩個高斯：牛市均值略正、波動低；熊市均值略負、波動高），用 Baum-Welch 估參數、再用濾波機率或 Viterbi 標記歷史——這套流程是計量經濟學裡 Hamilton 狀態切換模型的近親，實務上被拿來做風險開關（濾波出的熊市機率超過門檻就降部位）。三個職業級提醒：其一，狀態數 \(K\) 是超參數，加狀態永遠能提高樣本內概似，跟第 1 章的過擬合是同一隻怪物，要用樣本外驗證管住；其二，估出來的狀態沒有天生的名字，「狀態 0」是牛是熊要靠事後檢視發射參數來命名，且重新訓練後標籤可能互換；其三，濾波機率在狀態交界處必然猶豫（回想 t=5 那個 0.007092 對 0.007069 的瞬間），拿它當交易訊號要配上遲滯門檻，否則會在邊界上反覆進出。

八、濾波器校準：點校準與貝葉斯校準

到目前為止，模型參數（\(\phi\)、\(\sigma_v\)、轉移矩陣……）都是給定的。真實世界裡它們得從資料估出來——這一步叫校準（calibration），有兩條路。

點校準走最大概似：把參數固定，跑一遍濾波器；粒子濾波跑的過程中，每步未歸一化權重的平均恰好近似「當期觀測的預測概似」，全程取對數加總就是這組參數的對數概似——於是「跑一次濾波」等於「評一次概似」，外面再套一個數值優化器（如 BFGS 擬牛頓法）搜參數即可。一個實務暗礁：重抽樣是抽籤操作，會讓概似對參數的曲面帶鋸齒，優化器容易被卡住，文獻裡因此發展了連續化的重抽樣方案（原書 pp.231–233 有完整討論）。

貝葉斯校準則不滿足於一個點估計：給參數設先驗，用 MCMC（Metropolis-Hastings、Gibbs 抽樣這一族）從「參數與整條隱狀態路徑的聯合後驗」裡抽樣，最後拿到的是參數的整個分布——你不只知道 \(\phi\) 大概是 0.97，還知道你對這個數有多不確定。這與第 2、3 章的貝葉斯精神一脈相承：參數與隱狀態在後驗裡地位平等，一起被推斷。代價是計算重、收斂診斷麻煩；回報是不確定性被誠實地量化——對風險管理場景，這常常不是奢侈品而是必需品。

一句話收束本章：混合模型加上時間記憶得到 HMM；隱狀態連續化、保持線性高斯得到 Kalman；放棄線性高斯、改用樣本近似後驗得到粒子濾波。三件武器共享同一個世界觀——你看到的不是世界本身，是世界透過噪聲投下的影子；濾波，就是從影子倒推本體的紀律。

§03原書對照

本課以自創例子重組了原書第七章的骨架，以下是原書有、但本課未展開的內容，附 Springer 印刷頁碼供深入。其一，牛熊市場的完整數值例：原書給定特定的發射與轉移矩陣，示範如何沿聯合分解式逐項算出「一條指定隱狀態路徑配上指定觀測序列」的機率（pp.223–224，例 7.1、式 7.1–7.3）。其二，前向與後向量的正式定義、兩者乘積給出任意時點狀態後驗的恆等式，以及相鄰兩期聯合後驗的分解——這組量正是 Baum-Welch 的 E 步引擎，原書並將其定位為 EM 家族的非監督式演算法（p.224，式 7.4–7.6）。其三，作弊莊家例：莊家在公平硬幣與灌鉛硬幣間依轉移機率切換，觀測擲幣結果解碼莊家用的是哪枚硬幣——與本課牛熊例同構但更輕巧，適合當第二個練習場（pp.225–226，例 7.2）。其四，Viterbi 遞迴的原書記法與回溯函數的形式化定義，以及「轉移發射矩陣已知時不需要 Baum-Welch」的提醒（p.226）。其五，濾波與平滑的訊號處理視角：把兩者放進「從噪聲中萃取訊號」的一般框架，並比較逐點最大化與整體路徑最優的差異（p.226，§2.1.1）。其六，Kalman 濾波與單變量 RNN 的逐式對照：狀態方程非線性且同時吃前狀態與觀測、以及無噪聲這兩點差異的精確陳述（p.227，式 7.7–7.12）。其七，多項式重抽樣的完整演算法：累積權重和、均勻亂數落區間選粒子、複製因子的多項分布詮釋（p.229，§3.2）。其八，帶槓桿與跳躍的隨機波動度模型（SVLJ）：報酬與波動噪聲的相關結構捕捉槓桿效應、Bernoulli 跳躍計數器與跳躍強度參數，及其專屬的 SIR 變體——處理非加性、相關噪聲的重要性抽樣步驟（pp.230–231，式 7.13–7.15）。其九，粒子濾波的概似估計式：由未歸一化權重的逐期平均連乘出邊際概似，據以做頻率派 MLE 校準；含重抽樣造成目標函數不連續的實務討論與連續重抽樣的解法（pp.231–233，式 7.16）。其十，貝葉斯校準的完整管線：SVL 模型的先驗選取、聯合後驗的條件分解、BUGS 的有向無環圖表示、Metropolis 到 Metropolis-Hastings 到 Gibbs 的方法系譜，以及線上更新參數的適應性濾波文獻指引（pp.233–234）。其十一，四道習題：ARMA 模型的 Kalman 狀態空間表示、Ornstein-Uhlenbeck 過程的 Kalman 建模、SVLJ 粒子濾波的逐步推導、以及不誠實賭場的 Viterbi 解碼（pp.235–236）。原書第七章對應印刷頁 pp.221–238（含參考文獻）。

§04作業和解答

作業一：前向算法手算＋程式驗證

沿用課文的玩具市場（\(\pi=(0.75,0.25)\)、\(A\)、\(B\) 同第二節），觀測序列改為「跌、漲、跌」。(1) 手算 \(\alpha_1\) 與 \(\alpha_2\)（保留全部小數位）；(2) 寫程式算出 \(p(\mathbf{y})\)，並用暴力枚舉全部 \(2^3=8\) 條路徑驗證；(3) 說明前向算法的計算量為什麼是 \(T\times K^2\) 而不是 \(K^T\)。

解答 SOLUTION

(1) 第一天觀測到跌：\(\alpha_1(\text{牛}) = 0.75\times 0.35 = 0.2625\)；\(\alpha_1(\text{熊}) = 0.25\times 0.60 = 0.15\)。第二天觀測到漲，先轉移再發射：\(\alpha_2(\text{牛}) = (0.2625\times 0.85 + 0.15\times 0.10)\times 0.65 = 0.238125\times 0.65 = 0.15478125\)；\(\alpha_2(\text{熊}) = (0.2625\times 0.15 + 0.15\times 0.90)\times 0.40 = 0.174375\times 0.40 = 0.06975\)。

(2) 程式驗證：

PYTHON
import numpy as np
from itertools import product

pi = np.array([0.75, 0.25])
A  = np.array([[0.85, 0.15], [0.10, 0.90]])
B  = np.array([[0.65, 0.35], [0.40, 0.60]])
obs = [1, 0, 1]                              # 跌 漲 跌

alpha = pi * B[:, obs[0]]
print("alpha1 =", alpha)
for t in range(1, 3):
    alpha = (alpha @ A) * B[:, obs[t]]
    print(f"alpha{t+1} =", alpha)
print("P(y) =", alpha.sum())

total = sum(
    pi[p[0]] * B[p[0], obs[0]]
    * A[p[0], p[1]] * B[p[1], obs[1]]
    * A[p[1], p[2]] * B[p[2], obs[2]]
    for p in product([0, 1], repeat=3))
print("暴力枚舉 =", total)
# 輸出：
# alpha1 = [0.2625 0.15  ]
# alpha2 = [0.15478125 0.06975   ]
# alpha3 = [0.04848867 0.05159531]
# P(y) = 0.10008398437499999
# 暴力枚舉 = 0.100083984375

手算的 \(\alpha_1, \alpha_2\) 與程式完全一致，\(p(\mathbf{y}) \approx 0.1001\) 兩法吻合。有趣的細節：第三天結束時 \(\alpha_3(\text{熊}) = 0.0516 > \alpha_3(\text{牛}) = 0.0485\)——「跌漲跌」的證據已讓熊市小勝，即使開局先驗七成五偏牛。

(3) 關鍵是馬可夫性質帶來的「證據可合併」：所有在第 \(t\) 天抵達狀態 \(k\) 的路徑，之後的演化規律完全相同，所以它們的機率可以先加總成一個數 \(\alpha_t(k)\) 再一起前進，不必各走各的。每個時間步只需把 \(K\) 個舊值透過 \(K\times K\) 的轉移矩陣匯流成 \(K\) 個新值（\(K^2\) 次乘加），共 \(T\) 步，總量 \(T\times K^2\)。指數爆炸被動態規劃拆成了線性成長。

作業二：多看一天，歷史改寫——Viterbi 的非線上性格

課文第四節對「漲、漲、跌、跌、跌、漲」解碼得到整條牛市。現在市場多走了一天：第 7 天收跌，序列變成「漲、漲、跌、跌、跌、漲、跌」。(1) 對七天序列重新跑 Viterbi，比較前六天的解碼與原答案；(2) 解釋為什麼多一筆觀測能改寫「過去」的判讀；(3) 這對「拿 Viterbi 路徑當每日交易訊號」的做法是什麼警訊?

解答 SOLUTION

(1) 用課文的 Viterbi 程式改觀測序列重跑（obs = [0, 0, 1, 1, 1, 0, 1]），輸出：六天版本解碼為「牛牛牛牛牛牛」（最佳路徑機率 0.003918）；七天版本解碼為「牛牛熊熊熊熊熊」（最佳路徑機率 0.001374，終點分數 V(熊)=0.001374 > V(牛)=0.001166）。同樣的前六天資料，第 3 到第 6 天的標籤全部從牛翻成熊——包括第 6 天那個上漲日，如今被解讀成「熊市裡的反彈」。

(2) Viterbi 最優化的對象是整條路徑的聯合機率，而非任何單日的判讀。六天版本裡，「牛市回檔」路徑勝出的關鍵是最後一天的上漲替牛市扳回一城（課文已指出 t=5 時兩邊只差千分之三，天平隨時會倒）。第 7 天再跌一根，證據結構改變：「三跌、一漲、一跌」用「熊市＋一天反彈」解釋，整體代價低於「牛市＋四次逆發射」——於是全域最優路徑換邊，連帶改寫了它途經的每一天。回溯指標是從終點往回讀的，終點變了，讀出來的歷史就變了。

(3) 警訊：Viterbi 是離線（offline）演算法——它的每一段判讀都以「看過整條序列」為前提，新資料到來時歷史標籤沒有穩定性保證。若把它當每日訊號，等於讓今天的部位決策引用一個明天就可能被改寫的歷史敘事，回測時更會不小心引入前視偏差（第 1 章的老怪物）。線上決策該用的是濾波機率 \(p(s_t \mid y_{1:t})\)——它只用截至當下的資訊、且不會回頭改寫；Viterbi 留給事後的歷史分段報告。

作業三：粒子濾波收斂到 Kalman——用精確解當裁判

在線性高斯模型 \(x_t = 0.9\,x_{t-1} + \varepsilon_t\)（\(\varepsilon_t\sim\mathcal{N}(0, 0.5^2)\)）、\(y_t = x_t + \xi_t\)（\(\xi_t\sim\mathcal{N}(0, 1)\)）上，Kalman 濾波是精確解。請：(1) 模擬 200 期資料並跑 Kalman 濾波；(2) 在同一資料上跑 bootstrap 粒子濾波，粒子數 \(M \in \{100, 1000, 10000\}\)；(3) 以 Kalman 為裁判，計算各 \(M\) 之下粒子濾波估計與 Kalman 估計的平均絕對差，觀察收斂速度。

解答 SOLUTION
PYTHON
import numpy as np

rng = np.random.default_rng(3)
T = 200
phi, q, r = 0.9, 0.5, 1.0
x = np.zeros(T); y = np.zeros(T)
y[0] = x[0] + rng.normal(0, r)
for t in range(1, T):
    x[t] = phi * x[t-1] + rng.normal(0, q)
    y[t] = x[t] + rng.normal(0, r)

m, P = 0.0, q**2 / (1 - phi**2)       # Kalman：精確解
kal = np.zeros(T)
for t in range(T):
    if t > 0:
        m = phi * m; P = phi**2 * P + q**2
    K = P / (P + r**2)
    m = m + K * (y[t] - m); P = (1 - K) * P
    kal[t] = m

for M in (100, 1000, 10000):          # bootstrap 粒子濾波
    r2 = np.random.default_rng(9)
    parts = r2.normal(0, np.sqrt(q**2 / (1 - phi**2)), M)
    pf = np.zeros(T)
    for t in range(T):
        if t > 0:
            parts = phi * parts + r2.normal(0, q, M)
        logw = -0.5 * (y[t] - parts)**2 / r**2
        w = np.exp(logw - logw.max()); w /= w.sum()
        pf[t] = np.sum(w * parts)
        parts = parts[r2.choice(M, size=M, p=w)]
    print(f"M={M:6d}: 與 Kalman 的平均絕對差 {np.mean(np.abs(pf - kal)):.4f}")
# 輸出：
# M=   100: 與 Kalman 的平均絕對差 0.0799
# M=  1000: 與 Kalman 的平均絕對差 0.0231
# M= 10000: 與 Kalman 的平均絕對差 0.0079

粒子數每乘 10，平均絕對差約縮成三分之一——正是蒙地卡羅方法 \(1/\sqrt{M}\) 收斂律的指紋（\(\sqrt{10}\approx 3.16\)）。三個讀法：其一，粒子濾波確實在逼近同一個貝葉斯後驗，Kalman 是它的極限；其二，收斂是「根號慢」的——想把誤差再砍半，粒子數得翻四倍，這就是粒子濾波計算成本的來源；其三，反過來說，在這個 Kalman 能精確處理的問題上花一萬顆粒子是純浪費——工具的選擇次序永遠是「先問線性高斯成不成立」。

作業四：關掉重抽樣——親眼看粒子退化

把課文第七節的隨機波動度實驗改成全程不重抽樣：權重在對數空間逐期累加（logw += ...），不再執行 rng.choice 那一步。(1) 追蹤 ESS 在 t = 1、5、20、100、500 的變化；(2) 比較追蹤誤差與有重抽樣版本的差距；(3) 解釋 ESS 崩潰的機制。

解答 SOLUTION

把課文程式的權重更新改為累積式（logw = np.zeros(M) 初始化，每期 logw += -0.5 * (np.log(2 * np.pi * var) + ret[t]**2 / var)，歸一化只為了計算 ESS 與估計值，粒子永不重抽），其餘同課文（同 seed、同資料）。輸出：

PYTHON
# t=  1: ESS =   1893.3
# t=  5: ESS =    893.8
# t= 20: ESS =    275.8
# t=100: ESS =      4.7
# t=500: ESS =      1.0
# 不重抽樣 RMSE: 0.7537   （課文有重抽樣版本：0.5591）

(1) ESS 從近乎滿員的 1893 一路崩塌：20 期剩 276、100 期剩 4.7、500 期只剩 1.0——兩千顆粒子名存實亡，整個後驗由一顆粒子獨演。(2) 追蹤誤差從 0.5591 惡化到 0.7537；而且這還是運氣好的結果——那顆倖存粒子恰巧走了一條不太離譜的軌跡，換個 seed 可以慘得多，因為單粒子近似的變異數毫無控制。(3) 機制：不重抽樣時，粒子 \(i\) 的累積權重是 500 個單期概似的連乘，等於 500 期「運氣」的複利。各粒子軌跡的對數概似差距隨時間線性擴大，取指數後權重差距指數式擴大——任何微小的持續劣勢都會把權重壓到機器零。重抽樣的作用正是在差距失控前「結算」：淘汰落後者、複製領先者，讓所有算力回到當前後驗的高機率區重新出發。這也解釋了為什麼課文版本的 ESS 能全程維持在 128 以上。

§05參考資料