第 12b 章 随机波动率与 Markov 转换模型

本章对应 Tsay 原书第 12 章 12.7–12.11 节。第 03a 章的 GARCH 模型让波动率是过去收益的确定函数；随机波动率（stochastic volatility, SV）模型则让波动率有自己的随机冲击，更接近连续时间期权定价中的波动率过程，但似然要对整条潜在波动率路径积分，极大似然很难做。Markov 转换（Markov switching）模型让市场在几个"区制"之间随机切换，似然是对所有状态路径的混合。这两类模型正是 MCMC 与数据增广最能发挥作用的地方：把每天的波动率、每天的状态都当作参数抽样。本章先讲单变量与多元 SV 的 Gibbs 估计，再讲基于 Kalman 滤波的前向滤波后向抽样（FFBS）新方法与杠杆效应，然后是 Markov 转换 GARCH-M 模型，最后是 MCMC 预测与 VaR。

学习目标

1. 写出单变量 SV 模型，说明为什么它的极大似然困难，而 Gibbs 抽样容易。
2. 推导单个 \(h_t\) 的条件后验 (12.23)，会用 Griddy Gibbs 抽样；了解多元 SV 的 Cholesky 构造与时变相关。
3. 掌握 Kim–Shephard–Chib 七正态混合近似与前向滤波后向抽样（FFBS），理解为什么联合抽取整条波动率路径比逐点抽样高效；了解杠杆效应的处理方法。
4. 会写两状态 Markov 转换模型的 Gibbs 步骤（状态、转移概率、各状态参数），理解期望持续期 \(1/e_i\)。
5. 会在 MCMC 迭代中模拟未来路径得到预测分布，并据此计算纳入参数不确定性的 VaR；知道滤波概率与平滑概率在回测中的区别。

读前导读

这一章在解决什么问题。 你在 FRM/CFA 里见过 EWMA 和 GARCH：今天的方差由昨天的方差和昨天的收益平方决定，是一个确定的公式。随机波动率（SV）模型多加了一个假设：波动率自己也会受到随机冲击，哪怕今天市场风平浪静，明天的波动率也可能突然上升。这更像期权交易员眼中的世界（Heston 模型、VIX 自身也在波动），但代价是波动率完全看不见，估计时必须把整条波动率路径当作未知量。Markov 转换模型则说市场有几种"天气"——平静与动荡——在它们之间随机切换，每种天气下收益的均值和波动不同。对风控来说，它就是一个有统计基础的"风险开关"。

两类模型都有大量看不见的量（每天的波动率、每天的区制），所以都用上一章的工具：数据增广 + Gibbs 抽样。本章的技术亮点是 FFBS：借用第 11 章的 Kalman 滤波，一次性抽出整条波动率路径，而不是一天一天地抽。最后用 MCMC 做预测和 VaR，这是与风险管理工作直接相关的部分：预测分布里同时包含"未来冲击""今天波动率到底多大""参数估得准不准"三种不确定性，所以 VaR 会比只代入点估计的做法更保守。

需要先想起来的数学。

• 变量替换与雅可比。若 \(x=\ln h\) 的密度是 \(g(x)\)，则 \(h\) 的密度是 \(g(\ln h)\cdot\frac{d\ln h}{dh}=g(\ln h)/h\)。乘上的这个导数（雅可比）负责把"\(x\) 轴上的一小段"换算成"\(h\) 轴上的一小段"。对数正态分布的密度里那个 \(1/x\) 就是这么来的。见 第 00 册第 03 章 积分 和 第 00 册第 05 章 多元微积分与优化。
• 配方与两个正态的乘积。两个关于 \(x\) 的正态密度相乘仍是正态，新精度 = 两个精度之和，新均值 = 按精度加权。这是第 12a 章 Result 12.1 的核心，本章 (12.23)、(12.32) 都在用。见 第 00 册第 07 章 概率中的分析工具。
• Cholesky 分解。正定的协方差矩阵可以写成"下三角 × 对角 × 上三角"，相当于把两资产的协方差拆成"资产 1 的方差""资产 2 对资产 1 的回归系数""资产 2 的残差方差"。（正定：对任意非零向量 \(w\)，\(w'\Sigma w>0\)，即任何组合的方差都为正。）见 第 00 册第 06 章 线性代数速成。
• 一阶泰勒展开（线性化）。\(G(z)\approx G(z_0)+G'(z_0)(z-z_0)\)。久期就是债券价格对收益率的一阶泰勒近似；本章用同样的思路把非线性的波动率转移方程近似成线性，以便用 Kalman 滤波。见 第 00 册第 02 章 导数与泰勒展开。
• 几何分布与马尔可夫链的平稳分布。每期以概率 \(e\) 离开，停留期的期望是 \(1/e\)。见 第 00 册第 04 章 级数与收敛（几何级数求和）。

记号提示：\(\mathrm{tr}(A)\) 是矩阵的迹，即对角元之和；\(IG\) 是逆 Gamma 分布（若 \(1/\omega\) 服从 Gamma，则 \(\omega\) 服从逆 Gamma），用来描述方差；\(\ln\chi^2_1\) 指"一个标准正态的平方再取对数"的分布。

怎么读这一章。 必读：12.8.1–12.8.3（SV 模型与 Gibbs 框架）、12.9.3–12.9.4（KSC 混合与 FFBS，本章方法论核心）、12.10.1–12.10.2（Markov 转换的 Gibbs 步骤）、12.11（MCMC 预测与 VaR）以及两个 Python 示例。12.8.5–12.8.6 多元 SV 第一次可以只看 (12.28) 的 Cholesky 构造和"\(q_{21,t}\) 就是时变 β"这句话；12.9.2 第 4 条 \((\rho,\sigma_\eta^2)\) 的重参数化和 12.9.5 杠杆处理可以先跳过，需要估计杠杆时再回来。各例题的数值表格浏览即可，重点看结论段。

12.8 随机波动率模型

12.8.1 模型

第 03b 章 3.12 节已给出 SV 模型的形式并指出其似然难以计算，本节给出完整的估计方法。MCMC 在金融中的一个重要应用，是 Jacquier, Polson & Rossi (1994)（JPR）对 SV 模型的估计。单变量 SV 模型为

\(x_{it}\) 在 \(t-1\) 时已知（可以包含滞后收益），\(\epsilon_t\sim N(0,1)\)，\(v_t\sim N(0,\sigma_v^2)\)，两者独立。取对数保证 \(h_t>0\)；\(|\alpha_1|<1\) 保证对数波动率平稳；也可以用更高阶的 AR。

与 GARCH 的区别在于 (12.21) 里有一个独立的随机冲击 \(v_t\)：即使今天收益为零，明天的波动率也可能跳升。这使 SV 模型能产生更厚的尾部，并更贴近 Heston 一类连续时间随机波动率模型（第 06 章、第 08 册）。

金融直觉：GARCH 像一位只看历史成交的风险经理：昨天跌得多，今天就上调波动率估计，规则完全确定。SV 像是承认还有看不见的消息流：即便价格没动，某个新闻、流动性变化也可能让明天的波动率跳升，所以波动率本身是一个带噪声的随机过程。 厚尾的来源：给定 \(h_t\) 时收益是正态的，但 \(h_t\) 本身随机，收益就是"方差各不相同的正态分布的混合"。混合后的分布比同方差的单一正态尖峰、厚尾，这与股票收益的经验特征一致。

12.8.2 为什么用 Gibbs

参数分两组：\(\beta=(\beta_0,\dots,\beta_p)'\) 和 \(\omega=(\alpha_0,\alpha_1,\sigma_v^2)'\)。不可观测的波动率 \(H=(h_1,\dots,h_n)'\) 是辅助变量。极大似然困难，因为似然要对 \(n\) 维的 \(H\) 积分：

贝叶斯框架下，把 \(H\) 当作增广参数，先验取 \(p(\beta,\omega)=p(\beta)p(\omega)\)，Gibbs 依次抽取三组条件后验

给定 \(H\)，均值方程是已知异方差的回归，波动方程是一个 AR(1)，都是第 12a 章的老问题；难点只在抽 \(H\)。

白话解释：这个积分为什么算不动？\(dH\) 是对 \(h_1,\dots,h_n\) 这 \(n\) 个变量同时积分。月度数据 \(n=575\)，就是 575 重积分。如果每个维度只取 10 个格点做数值积分，也要 \(10^{575}\) 次计算。GARCH 没有这个问题，因为它的 \(h_t\) 由过去数据算出来，不是未知量。 Gibbs 的绕法：不去求边际似然 \(f(R|\beta,\omega)\)，而是把 \(H\) 和参数一起抽样；最后只看参数那几列样本，就等于"把 \(H\) 积掉"了。这是 MCMC 的通用技巧：对样本"忽略某些列"就是对它们求边际分布。

12.8.3 单变量模型的条件后验

\(\beta\) 的条件后验。给定 \(H\)，两边除以 \(\sqrt{h_t}\) 得同方差回归：

先验 \(\beta\sim N(\beta_o,A_o)\)，后验 \(N(\beta_*,A_*)\)：\(A_*^{-1}=\sum x_{o,t}x_{o,t}'+A_o^{-1}\)，\(\beta_*=A_*(\sum x_{o,t}r_{o,t}+A_o^{-1}\beta_o)\)。这就是加权最小二乘加正态先验。

\(h_t\) 的条件后验（逐点抽取）。\(h_t\) 只出现在三处：观测方程 \(a_t|h_t\)，以及波动方程中的 \(\ln h_t|\ln h_{t-1}\) 与 \(\ln h_{t+1}|\ln h_t\)。因此对 \(1<t<n\)，

推导要点：

1. \(a_t|h_t\sim N(0,h_t)\) 贡献 \(h_t^{-1/2}\exp[-a_t^2/(2h_t)]\)；
2. 把 \(\ln h_t\) 看作变量，两个正态因子的指数为 \(-\frac{1}{2\sigma_v^2}[(\ln h_t-a)^2+\alpha_1^2(\ln h_t-b)^2]\)，其中 \(a=\alpha_0+\alpha_1\ln h_{t-1}\)，\(b=(\ln h_{t+1}-\alpha_0)/\alpha_1\)；
3. 用配方恒等式 \((x-a)^2A+(x-b)^2C=(x-c)^2(A+C)+(a-b)^2\frac{AC}{A+C}\)，\(c=\frac{Aa+Cb}{A+C}\)（Box & Tiao 1973 引理的标量版），取 \(A=1\)，\(C=\alpha_1^2\)，得 \(c=\mu_t\)，方差 \(\sigma_v^2/(1+\alpha_1^2)\)；
4. 从 \(\ln h_t\) 换回 \(h_t\) 的雅可比 \(d\ln h_t=h_t^{-1}dh_t\) 再贡献一个 \(h_t^{-1}\)，合计 \(h_t^{-1.5}\)。

（精读笔记记录原书第 2 步中 \(a\) 漏写了 \(\alpha_1\)，这里已更正。）

推导拆解：第 2 步里 \(\alpha_1^2(\ln h_t-b)^2\) 是怎么来的。波动方程在 \(t+1\) 期写作 \(\ln h_{t+1}=\alpha_0+\alpha_1\ln h_t+v_{t+1}\)，它对 \(\ln h_t\) 贡献的指数是 \(-(\ln h_{t+1}-\alpha_0-\alpha_1\ln h_t)^2/(2\sigma_v^2)\)。把括号里的 \(\alpha_1\) 提出来：\(\ln h_{t+1}-\alpha_0-\alpha_1\ln h_t=-\alpha_1\big(\ln h_t-\frac{\ln h_{t+1}-\alpha_0}{\alpha_1}\big)=-\alpha_1(\ln h_t-b)\)，平方后得到 \(\alpha_1^2(\ln h_t-b)^2\)。 直观上，前一天的 \(h_{t-1}\) 给出一个对 \(\ln h_t\) 的"预测" \(a\)（精度 \(1/\sigma_v^2\)），后一天的 \(h_{t+1}\) 给出一个"倒推" \(b\)（精度 \(\alpha_1^2/\sigma_v^2\)），两者按精度加权就是 \(\mu_t\)。再乘上当天收益 \(a_t\) 提供的信息，就是完整的条件后验。

白话解释：第 4 步的雅可比为什么必要。前三步得到的是 \(\ln h_t\) 的密度，而我们要对 \(h_t\) 本身做 Griddy Gibbs。同样宽度的区间，在 \(h\) 很小时对应的 \(\ln h\) 区间很宽，在 \(h\) 很大时对应的 \(\ln h\) 区间很窄，所以从"\(\ln h\) 的密度"换成"\(h\) 的密度"必须乘以 \(d\ln h/dh=1/h\) 来校正。漏掉它，抽到的 \(h_t\) 会系统性偏大。

(12.23) 不是标准分布。JPR 用 Metropolis 算法抽 \(h_t\)；原书用 Griddy Gibbs，\(h_t\) 的格点范围取样本无条件方差的若干倍。

另一个视角：(12.23) 也可以由第 12a 章 AR(1) 缺失值的结果得到。把 \(\ln h_t\) 当作缺失值，构造两观测的回归：\(y_t=\ln h_t+b_t\)（\(y_t=\alpha_0+\alpha_1\ln h_{t-1}\)）与 \(y_{t+1}=\alpha_1\ln h_t+b_{t+1}\)（\(y_{t+1}=\ln h_{t+1}-\alpha_0\)）。其 LS 估计正是 \(\mu_t\)，方差 \(\sigma_v^2/(1+\alpha_1^2)\)。于是 (12.23) 就是"\(a_t\sim N(0,h_t)\)"与"\(\widehat{\ln h_t}\sim N(\ln h_t,\sigma_v^2/(1+\alpha_1^2))\)"两部分信息的乘积再做变量变换。这样看也容易推广到 AR(\(p\))（假设前 \(p\) 个 \(h_t\) 固定）。\(h_t\) 的初值可以用第 03a 章 GARCH 模型的拟合值。

端点。(12.23) 只适用于 \(1<t<n\)。简单处理：固定 \(h_1\)；\(t=n\) 时只用 \(\ln h_n\sim N(\alpha_0+\alpha_1\ln h_{n-1},\sigma_v^2)\)。也可以用 \(h_{n+1}\) 的预测（原点 \(n-1\) 的 2 步预测 \(\alpha_0+\alpha_1(\alpha_0+\alpha_1\ln h_{n-1})\)）和 \(h_0\) 的后向预测：高斯 AR(1) 时间可逆，\(\ln h_t-\eta=\alpha_1(\ln h_{t+1}-\eta)+v_t^*\)，\(\eta=\alpha_0/(1-\alpha_1)\)，所以 2 步后向预测为 \(\eta+\alpha_1^2(\ln h_2-\eta)\)。

\(\omega\) 的条件后验。分成 \(\alpha=(\alpha_0,\alpha_1)'\) 与 \(\sigma_v^2\) 两块。给定 \(H\)，\(\ln h_t\) 就是一个 AR(1)。先验 \(\alpha\sim N(\alpha_o,C_o)\)，后验

\(v_t=\ln h_t-\alpha_0-\alpha_1\ln h_{t-1}\)，先验 \(m\lambda/\sigma_v^2\sim\chi^2_m\)，后验 \((m\lambda+\sum_{t=2}^nv_t^2)/\sigma_v^2\sim\chi^2_{m+n-1}\)。

12.8.4 例 12.3：S&P 500 月度收益的 GARCH 与 SV

数据为 S&P 500 月度对数收益（%），1962-01 至 2009-12，575 个观测（用每月第一个交易日的收盘指数）。

高斯 GARCH(1,1)：

t 值都大于 2.56；残差 \(Q(12)=10.04\)（p=0.61），平方残差 \(Q(12)=6.14\)（p=0.91）。

SV 模型：\(r_t=\mu+a_t\)，\(a_t=\sqrt{h_t}\epsilon_t\)，\(\ln h_t=\alpha_0+\alpha_1\ln h_{t-1}+v_t\)（12.27）。先验 \(\mu\sim N(0,4)\)，\(\alpha\sim N[(0,0.6)',\mathrm{diag}(0.25,0.04)]\)，\(10\times0.1/\sigma_v^2\sim\chi^2_{10}\)。初值：\(h_{0t}\) 取 GARCH 拟合值，\(\alpha,\sigma_v^2\) 取 \(\ln h_{0t}\) 的 LS 估计，\(\mu\) 取样本均值。\(h_t\) 用 400 个格点的 Griddy Gibbs，第 \(j\) 次迭代的格点范围为 \([\eta_{1t},\eta_{2t}]\)。精读笔记记录原文为 \(\eta_{1t}=0.6\max(h_{j-1,t},h_{0t})\)，\(\eta_{2t}=1.4\min(h_{j-1,t},h_{0t})\)；当两者相差超过 2 倍多时下界会大于上界，按意图应为 \(\eta_{1t}=0.6\min(\cdot)\)、\(\eta_{2t}=1.4\max(\cdot)\)，即覆盖两个参考值的较宽区间。迭代 2500 次，丢弃前 500 次。

后验均值（后验标准差）：

\(\alpha_1=0.837\) 说明波动率持续性强，但小于 JPR 用日数据所得的值（月度数据的持续性本来就比日度低）。原书图 12.5 比较了先验（虚线，较无信息）与后验（实线）密度，\(\mu\) 和 \(\sigma_v^2\) 的后验很集中；图 12.6 显示 SV 的后验均值波动率与 GARCH 拟合波动率形态相似。换初值、先验、迭代次数结果都稳定。Griddy Gibbs 的结果和效率依赖 \(h_t\) 范围的设定，这是逐点抽样的一个弱点。

12.8.5 多元随机波动率模型

用第 10 章的 Cholesky 分解构造二元 SV。令 \(b_{1t}=a_{1t}\)，\(b_{2t}=a_{2t}-q_{21,t}b_{1t}\)（即回归 \(a_{2t}=q_{21,t}a_{1t}+b_{2t}\)），则

\(g_{ii,t}=\mathrm{Var}(b_{it}|F_{t-1})\)，\(b_{1t}\) 与 \(b_{2t}\) 不相关。

推导拆解：把 (12.28) 乘出来（省略下标 \(t\)），

\[\Sigma=\begin{bmatrix}g_{11}&q_{21}g_{11}\\q_{21}g_{11}&q_{21}^2g_{11}+g_{22}\end{bmatrix}.\]

逐项看：资产 1 的方差是 \(g_{11}\)；协方差 \(=q_{21}g_{11}\)，所以 \(q_{21}=\mathrm{Cov}/\mathrm{Var}_1\)，正是回归系数 β；资产 2 的方差 = 系统性部分 \(\beta^2\sigma_M^2\) + 特质部分 \(g_{22}\)，就是单指数模型的方差分解。相关系数 \(\rho_{21}=q_{21}\sqrt{g_{11}}/\sqrt{q_{21}^2g_{11}+g_{22}}\)。 好处是：只要 \(g_{11},g_{22}>0\)，不管 \(q_{21}\) 取什么实数，\(\Sigma\) 都自动正定，所以三者可以各自独立地建 AR(1) 模型，而不用担心组合方差算出负数。

模型为

\(v_{1t},v_{2t},u_t\) 为独立高斯白噪声，方差 \(\sigma_{iv}^2\) 与 \(\sigma_u^2\)。传统参数为 \(\beta,\alpha_i,\sigma_{iv}^2,\gamma,\sigma_u^2\)；增广参数为 \(Q=\{q_{21,t}\}\)、\(G_1=\{g_{11,t}\}\)、\(G_2=\{g_{22,t}\}\)。注意 \(q_{21,t}\) 正是"资产 2 对资产 1 的时变回归系数"，若资产 1 是市场，它就是时变 β。

Gibbs 步骤：(1) 用 (12.22) 逐行抽 \(\beta_0,\beta_1\)；(2) 用 (12.23)（\(a_t\) 换成 \(a_{1t}\)）抽 \(g_{11,t}\)，并像单变量那样抽 \(\alpha_1,\sigma_{1v}^2\)；(3) 先算 \(b_{2t}=a_{2t}-q_{21,t}a_{1t}\sim N(0,g_{22,t})\)，再同法抽 \(g_{22,t}\)、\(\alpha_2\)、\(\sigma_{2v}^2\)。另需：

• \(\gamma\)：给定 \(Q\) 是 AR(1)，先验 \(N(\gamma_o,D_o)\)，\(D_*^{-1}=\sum z_tz_t'/\sigma_u^2+D_o^{-1}\)，\(\gamma_*=D_*(\sum z_tq_{21,t}/\sigma_u^2+D_o^{-1}\gamma_o)\)，\(z_t=(1,q_{21,t-1})'\)；
• \(\sigma_u^2\)：\((m\lambda+\sum u_t^2)/\sigma_u^2\sim\chi^2_{m+n-1}\)，\(u_t=q_{21,t}-\gamma_0-\gamma_1q_{21,t-1}\)；
• \(q_{21,t}\)：

\(\mu_t=\frac{\gamma_0(1-\gamma_1)+\gamma_1(q_{21,t-1}+q_{21,t+1})}{1+\gamma_1^2}\)，\(\sigma^2=\sigma_u^2/(1+\gamma_1^2)\)。第一个因子作为 \(q_{21,t}\) 的函数，是均值 \(a_{2t}/a_{1t}\)、方差 \(g_{22,t}/a_{1t}^2\) 的正态；第二个因子是 \(N(\mu_t,\sigma^2)\)。由 Result 12.1 得闭式正态后验：

直观上，\(a_{2t}/a_{1t}\) 是"只用今天一天数据算的 β"，它的精度与 \(a_{1t}^2\) 成正比：市场今天波动越大，这一天对 β 的信息就越多。

12.8.6 例 12.4：IBM 与 S&P 500

数据为 IBM 与 S&P 500 月度对数收益，1962-01 至 2009-12，\(r_t=(IBM_t,SP_t)'\)。比较三个模型。

Cholesky 时变相关 GARCH：\(r_t=\beta_0+a_t\)（12.33），\(g_{11,t}=\alpha_{10}+\alpha_{11}g_{11,t-1}+\alpha_{12}a_{1,t-1}^2\)（12.34），\(g_{22,t}=\alpha_{20}+\alpha_{22}b_{2,t-1}^2\)（12.35），\(q_{21,t}=\gamma_0\)（12.36）。估计值（标准误）：\(\beta_{01}=0.69(0.30)\)，\(\beta_{02}=0.49(0.18)\)，\(\alpha_{10}=3.98(1.22)\)，\(\alpha_{11}=0.80(0.04)\)，\(\alpha_{12}=0.12(0.03)\)，\(\alpha_{20}=10.67(0.53)\)，\(\alpha_{22}=0.12(0.04)\)，\(\gamma_0=0.37(0.01)\)。

BEKK(1,1)（Matlab 估计）：\(\beta_0=(0.70,0.54)'\)，\(A=\begin{bmatrix}0.80&0\\0.83&0.01\end{bmatrix}\)，\(A_1=\begin{bmatrix}0.07&0.33\\-0.06&0.43\end{bmatrix}\)，\(B_1=\begin{bmatrix}1.00&-0.12\\0.01&0.90\end{bmatrix}\)。

SV 模型：均值方程同上，\(\ln g_{ii,t}=\alpha_{i0}+\alpha_{i1}\ln g_{ii,t-1}+v_{it}\)（12.37–12.38），\(q_{21,t}=\gamma_0+u_t\)（12.39）。先验 \(\beta_{i0}\sim N(0,4)\)，\(\alpha_i\sim N[(0,0.7)',\mathrm{diag}(0.25,0.04)]\)，\(\gamma_0\sim N(0,1)\)，\(10\times0.1/\sigma_{iv}^2\sim\chi^2_{10}\)，\(5\times0.2/\sigma_u^2\sim\chi^2_5\)。初值取自 BEKK，\(t=1\) 的值固定，迭代 2500 次丢弃 500 次，\(g_{ii,t}\) 用 500 格点的 Griddy Gibbs。后验均值（标准差）：

收敛检查：用不同初值与迭代次数得到的两个 Gibbs 样本（500+2000 与 500+1000），\(g_{11,t},g_{22,t},q_{21,t},\sigma_{22,t},\sigma_{21,t},\rho_{21,t}\) 的后验均值散点都贴近 \(y=x\)（原书图 12.8）。

三模型比较：

1. 均值方程基本相同。
2. IBM 的条件方差（图 12.9）三者都显示波动聚集和上升趋势，但 GARCH 的峰值更高，且在 1993 年多出一个峰。
3. S&P 500 的条件方差（图 12.10）：二元 GARCH 在 1993 年附近出现一个额外的峰，而单变量分析（图 12.6）没有。这是二元 GARCH 通过 Cholesky 结构依赖 IBM 收益而产生的虚假波动峰；SV 与 BEKK 都没有，SV 的 S&P 波动与单变量结果相似。
4. 条件相关（图 12.11）差异很大：Cholesky-GARCH 的相关平滑且恒为正，均值 0.59、标准差 0.07、范围 (0.411, 0.849)；BEKK 在 1993 年附近出现小的负值，均值 0.59、标准差 0.13、范围 (−0.020, 0.877)；SV 逐月剧烈变化，均值 0.60、标准差 0.14、范围 (−0.161, 0.839)，若干孤立时期为负——因为 SV 中 \(q_{21,t}\) 含随机冲击 \(u_t\)。

这些 Gibbs 方法可以推广到其他二元 SV 模型，条件后验需要扩展，但思路不变。

12.9 SV 估计的新方法：FFBS 与杠杆效应

逐点抽 \(h_t\) 有两个缺点：相邻的 \(h_t\) 高度相关，逐点抽样混合很慢；每个 \(h_t\) 都要做一次 Griddy Gibbs，计算量大。这一节的方法在 Kalman 滤波框架内用前向滤波后向抽样（forward filtering and backward sampling, FFBS）联合抽取整条波动率路径，大幅提高效率，还能推广到带杠杆效应与跳跃的随机扩散模型。

12.9.1 重参数化与杠杆效应

\(z_t\) 是零均值的对数波动率，\((\epsilon_t,\eta_t)'\) 服从二元正态，协方差 \(\Sigma=\begin{bmatrix}1&\rho\sigma_\eta\\\rho\sigma_\eta&\sigma_\eta^2\end{bmatrix}\)。\(\rho\) 刻画杠杆效应（leverage effect），通常为负：负收益推高未来波动。与 12.8 节的关系是 \(z_t=\ln h_t-\ln\sigma_0^2\)，\(\sigma_0^2=\exp\{E[\ln h_t]\}\)。

关键是时间下标：\(\eta_t\) 是 \(z_{t+1}\) 的新息，与 \(z_t\) 独立。正是这一错位让杠杆效应可以处理——今天的收益冲击 \(\epsilon_t\) 与明天的波动冲击 \(\eta_t\) 相关。若写成 \(z_t=\alpha z_{t-1}+\eta_t\) 且 \(\eta_t\) 与 \(\epsilon_t\) 相关，\(z_t\) 就与 \(\epsilon_t\) 相关，会产生识别问题。等价的写法有 \(r_t=x_t'\beta+\sigma_0\exp(z_{t-1}/2)\epsilon_t\)，\(z_t=\alpha z_{t-1}+\eta_t\)；或 \(r_t=x_t'\beta+\exp(z_{t-1}^*/2)\epsilon_t\)，\(z_t^*=\alpha_0+\alpha z_{t-1}^*+\eta_t\)，此时 \(E(z_t^*)=\alpha_0/(1-\alpha)\ne0\)。

12.9.2 条件后验

参数为 \(\beta,\sigma_0,\alpha,\rho,\sigma_\eta\) 与 \(z=(z_1,\dots,z_n)'\)（设 \(z_1\) 已知）。

1. \(\beta\)：同 12.8.3，把 \(\sqrt{h_t}\) 换成 \(\sigma_0\exp(z_t/2)\)。
2. \(\alpha\)：给定 \(z\) 与 \(\sigma_\eta^2\)，是 AR(1) 系数，正态先验下后验为正态。
3. \(\sigma_0^2\)：\(v_t=(r_t-x_t'\beta)\exp(-z_t/2)=\sigma_0\epsilon_t\) 是 iid \(N(0,\sigma_0^2)\)。先验 \(m\lambda/\sigma_0^2\sim\chi^2_m\)，后验 \((m\lambda+\sum_{t=1}^nv_t^2)/\sigma_0^2\sim\chi^2_{m+n}\)。
4. \((\rho,\sigma_\eta^2)\)：给定其余，可算出 \(b_t=(\epsilon_t,\eta_t)'\)，似然 \(\propto|\Sigma|^{-N/2}\exp[-\frac12\mathrm{tr}(\Sigma^{-1}\sum b_tb_t')]\)，\(N\) 为 \(b_t\) 的个数。\(\rho\) 与 \(\sigma_\eta^2\) 在这个形式里纠缠在一起。Jacquier, Polson & Rossi (2004) 的重参数化是

于是 \(|\Sigma|=\omega\)，\(\Sigma^{-1}=\frac1\omega S+\begin{bmatrix}1&0\\0&0\end{bmatrix}\)，\(S=\begin{bmatrix}\phi^2&-\phi\\-\phi&1\end{bmatrix}\)。令 \(e=(\epsilon_2,\dots,\epsilon_n)'\)，\(\eta=(\eta_2,\dots,\eta_n)'\)，\(R=(e,\eta)'(e,\eta)\)，似然化为

取共轭先验 \(\omega\sim IG(\gamma_0/2,\gamma_1/2)\)，\(\phi|\omega\sim N(0,\omega/2)\)，配方得

其中 \(\omega\) 的后验已把 \(\phi\) 积掉（先抽 \(\omega\)、再抽 \(\phi|\omega\)，是一次分块抽样）。原书（PDF p.673）印的形状参数是 \((n+1+\gamma_0)/2\)，与上面按 \(N\) 个配对推出的 \((N+\gamma_0)/2\) 不同。按 \(R=\sum_{t=2}^n\) 只有 \(n-1\) 对、\(\phi\) 的正态先验再贡献 \(1/2\) 计算，应得 \((n+\gamma_0)/2\)，原书多出的 1 很可能是计数约定或排印问题；样本较大时影响可以忽略。抽完后令 \(\sigma_\eta^2=\omega+\phi^2\)，\(\rho=\phi/\sigma_\eta\)。这里 \(IG(\alpha,\beta)\) 的密度为 \(f(\omega)=\frac{\beta^\alpha}{\Gamma(\alpha)}\omega^{-(\alpha+1)}e^{-\beta/\omega}\)。

5. \(z\) 的联合抽取：见下两小节。

12.9.3 线性化与 KSC 混合近似

由 (12.40)，\((r_t-x_t'\beta)^2/\sigma_0^2=\exp(z_t)\epsilon_t^2\)。取对数，令 \(y_t=\ln[(r_t-x_t'\beta)^2/\sigma_0^2]\)，

以 (12.42) 为观测方程、(12.41) 为状态方程，就是一个线性状态空间模型（精读笔记指出原文此处误把状态方程写成 12.40）。但 \(\epsilon_t^*\sim\ln\chi^2_1\) 不是正态：它左偏，均值约 −1.27，方差 \(\pi^2/2\approx4.93\)。Kim, Shephard & Chib (1998)（KSC）用 7 个正态的混合来近似：

（表 12.3；另见 Chib, Nardari & Shephard 2002。原书图 12.12 用 100,000 个观测比较，\(\ln\chi^2_1\) 的密度与混合密度几乎重合。）

白话解释：这一步的逻辑是"先变成线性，再变成高斯"。

1. 原模型里 \(z_t\) 以 \(\exp(z_t/2)\) 乘在噪声上，是非线性的。对收益平方取对数，乘法变加法：\(\ln[\exp(z_t)\epsilon_t^2]=z_t+\ln\epsilon_t^2\)，观测方程变成线性。
2. 代价是噪声 \(\ln\epsilon_t^2\) 不是正态：\(\epsilon_t\) 接近 0 时 \(\ln\epsilon_t^2\) 会趋向负无穷，所以左尾很长。
3. 任何形状的分布都能用几个正态的加权和近似，KSC 选了 7 个。表里权重大的几个分量（第 5、6、7 个）描述主体，权重很小、均值很负的分量（第 1、3 个）专门负责那条长左尾。
4. 再给每天贴一个"属于第几个分量"的标签 \(I_t\)。标签一旦给定，当天的噪声就是普通正态，整个模型就是第 11 章的线性高斯状态空间模型。 金融类比：这就像把一个复杂的损失分布拆成几个情景（正常、轻度压力、极端压力），每个情景内部用正态近似，再按情景概率加权。

再引入独立的混合指示变量 \(I_t\in\{1,\dots,7\}\)（又一次数据增广）。给定 \(I_t=i\)，\(\epsilon_t^*\sim N(\mu_i,\varpi_i^2)\)，模型就成了线性高斯状态空间模型：

（这里先考虑无杠杆的情形，\(\eta_t\) 与 \(e_t\) 不相关。）其 Kalman 滤波为

混合指示变量的条件后验：给定 \(z_t\)，\(I_t=i\) 的似然是正态密度 \(q_{it}=\varpi_i^{-1}\phi[(y_t-z_t-\mu_i)/\varpi_i]\)（\(\phi\) 为标准正态密度；精读笔记记录原文写作标准正态 CDF，按似然的含义应为密度，且要带 \(\varpi_i^{-1}\)），以 \(p_i\) 为先验，后验 \(p_{it}=p_iq_{it}/\sum_jp_jq_{jt}\)。抽到 \(I_t=j\) 就令 \(c_t=\mu_j\)，\(H_t=\varpi_j^2\)。

12.9.4 前向滤波后向抽样（FFBS）

为什么一定要化成高斯状态空间模型？因为这样可以联合、高效地抽取整条 \(z\)。利用马尔可夫性——给定 \(z_{t+1}\)，\(z_t\) 与 \(z_{t+j}\)（\(j>1\)）独立——把联合后验分解为

第一个因子 \(p(z_n|F_n)=N(z_{n|n},\Sigma_{n|n})\) 由 Kalman 滤波给出。关键性质是

因为给定 \(z_n\) 后，\(v_n=y_n-c_n-z_{n|n-1}\) 只依赖 \(z_n\) 与 \(e_n\)，不再带有 \(z_{n-1}\) 的信息。由 Kalman 滤波，

由定理 11.1，\(p(z_{n-1}|z_n,F_n)=N(\mu_{n-1}^*,\Sigma_{n-1}^*)\)（12.50），

一般地，\(p(z_t|z_{t+1},F_n)=p(z_t|z_{t+1},F_t)=N(\mu_t^*,\Sigma_t^*)\)（12.51），

算法：给定 \(z_{1|0},\Sigma_{1|0}\)，前向运行 Kalman 滤波并保存 \(z_{t|t},\Sigma_{t|t}\)；从 \(N(z_{n|n},\Sigma_{n|n})\) 抽 \(z_n\)，再依次按 (12.51) 抽 \(z_{n-1},\dots,z_1\)，得到 \(z\) 的一个联合实现。这就是 FFBS（Carter & Kohn 1994；Frühwirth-Schnatter 1994）。因为 \(z_t\) 高度序列相关，联合抽取比逐点抽取高效得多。

白话解释：(12.46) 用的是概率的乘法规则：\(p(z_1,\dots,z_n)=p(z_n)\,p(z_{n-1}|z_n)\,p(z_{n-2}|z_{n-1},z_n)\cdots\)，再用马尔可夫性把 \(p(z_{n-2}|z_{n-1},z_n)\) 简化成 \(p(z_{n-2}|z_{n-1})\)（全部都在 \(F_n\) 条件下）。于是"抽一条 \(n\) 维路径"变成"从最后一天开始，一天一天往回抽一维正态"，每一步只需要前向滤波存下的 \(z_{t|t},\Sigma_{t|t}\)。 与第 11b 章平滑器的区别：平滑器给出每天的后验均值 \(z_{t|n}\) 和方差，是"最可能的那条路径"附近的描述；FFBS 抽出的是一条随机的完整路径，各天之间的相关性被正确保留。Gibbs 需要的正是后者。 为什么逐点抽样慢：\(\alpha=0.95\) 时，给定前后两天，\(z_t\) 的条件方差只有 \(\sigma_\eta^2/(1+\alpha^2)\)，几乎被邻居"钉死"。整条路径要整体上移，得靠每个点一次挪一小步，需要成千上万次迭代。FFBS 每次迭代都直接抽一条全新的路径，没有这个问题。

推导拆解：(12.51) 怎么来。在 \(F_t\) 下，\(z_t\sim N(z_{t|t},\Sigma_{t|t})\)，\(z_{t+1}=\alpha z_t+\eta_t\)，所以 \(z_{t+1}\) 的均值 \(\alpha z_{t|t}=z_{t+1|t}\)，方差 \(\Sigma_{t+1|t}\)，协方差 \(\mathrm{Cov}(z_t,\alpha z_t+\eta_t)=\alpha\Sigma_{t|t}\)。套条件正态公式（均值 + 协方差/方差 × 偏离；方差 − 协方差²/方差）即得 \(\mu_t^*\) 与 \(\Sigma_t^*\)。这和 CFA 里"已知市场收益后，个股收益的条件期望 = 均值 + β × 市场偏离"是同一个公式。

FFBS 适用于任何线性高斯状态空间模型，核心恒等式是 \(p(S_t|S_{t+1},F_n)=p(S_t|S_{t+1},F_t,v_{t+1},\dots,v_n)=p(S_t|S_{t+1},F_t)\)。它也是第 11a、11b 章 Kalman 滤波在贝叶斯计算中最重要的用途之一。

12.9.5 杠杆效应的处理

(12.42) 的平方变换丢掉了 \(\epsilon_t\) 的符号，也就丢掉了它与 \(\eta_t\) 的相关，因此上面的方法估计不了杠杆。Artigas & Tsay (2004) 的做法是：当 \(\rho\ne0\)，把 \(\eta_t\) 分解为 \(\eta_t=\rho\sigma_\eta\epsilon_t+\eta_t^*\)，\(\eta_t^*\) 与 \(\epsilon_t\) 独立，\(\mathrm{Var}(\eta_t^*)=\sigma_\eta^2(1-\rho^2)\)。代入 \(\epsilon_t=(r_t-x_t'\beta)\exp(-z_t/2)/\sigma_0\)：

转移方程变成非线性，Kalman 滤波不再直接适用。用时变线性化（类似扩展 Kalman 滤波，extended Kalman filter）修改 (12.45) 的最后两式：

\(g(z_{t|t})=\partial G(x)/\partial x|_{x=z_{t|t}}\)，在滤波状态处求导（精读笔记指出原文称之为"平滑状态"，实际是滤波状态）。

推导拆解：

1. \(\eta_t\) 的分解就是回归：\(\eta_t\) 对 \(\epsilon_t\) 回归，系数 \(=\mathrm{Cov}(\eta_t,\epsilon_t)/\mathrm{Var}(\epsilon_t)=\rho\sigma_\eta\)，残差 \(\eta_t^*\) 方差 \(=\sigma_\eta^2-(\rho\sigma_\eta)^2=\sigma_\eta^2(1-\rho^2)\)。
2. 记 \(k_t=\rho\sigma_\eta(r_t-x_t'\beta)/\sigma_0\)（给定数据和参数是已知数），则 \(G(z)=\alpha z+k_te^{-z/2}\)。求导：\(g(z)=\alpha-\frac{k_t}{2}e^{-z/2}\)（\(e^{-z/2}\) 的导数用链式法则得 \(-\frac12e^{-z/2}\)）。
3. 在 \(z_{t|t}\) 附近一阶泰勒展开 \(G(z_t)\approx G(z_{t|t})+g(z_{t|t})(z_t-z_{t|t})\)。取期望得预测均值 \(G(z_{t|t})\)；方差按"线性变换的方差 = 系数² × 原方差"得 \(g^2\Sigma_{t|t}\)，再加上 \(\eta_t^*\) 的方差。这就是 (12.53)。 金融含义：\(\rho<0\) 时，若今天收益为负（\(r_t-x_t'\beta<0\)），\(k_t>0\)，明天的对数波动率被推高，这就是杠杆效应。用一阶近似处理非线性，与用久期近似债券价格变化同理：变动小时很准，变动大时有误差（凸性被忽略）。

12.9.6 例 12.5：带杠杆的 SV

数据为 S&P 500 月度对数收益，1962-01 至 2004-11，515 个观测（例 12.3 的子序列）。模型

模型 1 无杠杆，模型 2 有 \(\mathrm{corr}(\epsilon_t,\eta_t)=\rho\)。用 Matlab 程序做 FFBS–Gibbs，共 10000 次迭代，前 2000 次为 burn-in。原书表 12.4 的后验均值（标准差）：

\(\hat\rho=-0.39\) 与文献中常见的杠杆系数接近。原书图 12.13 显示两个模型的后验均值波动率非常接近，形态和量级都与例 12.3 相似（注意图 12.6 是百分比收益的条件方差，图 12.13 是对数收益的条件标准差）。

阅读提示：表中 \(\alpha\approx-0.06\)、\(\sigma_\eta\approx2.56\) 意味着对数波动率几乎没有持续性、冲击极大，这与例 12.3 中 \(\alpha_1=0.837\) 的高持续性差异很大，也与"波动聚集"的经验事实不符。已对照原书（PDF p.679 表 12.4）：两个模型的 \(\alpha\) 确实印为 \(-0.0616\) 和 \(-0.0613\)、标准误约 0.12，\(\sigma_\eta\) 为 2.56 和 2.58，并非抽取错误。它与常见的高持续性 SV 估计差异很大，原书未解释；可能与月度数据、\(\sigma_o\) 的标度约定有关。把它当作原书的一个结果记录即可，不宜作为 SV 参数的典型值。

12.10 Markov 转换模型

12.10.1 模型

MCMC 相对传统似然方法的另一个优势场景是 Markov 转换模型。McCulloch & Tsay (1994b) 用 Gibbs 抽样估计各状态内波动率不变的 Markov 转换模型，应用于美国季调实际 GNP 季度增长率，发现"收缩"与"扩张"期的动态显著不同（第 04a 章 4.5 节从极大似然角度介绍过这类模型，转移概率记作 \(w_1,w_2\)，与下文的 \(e_1,e_2\) 含义相同）。本节关注波动率的转换。两状态、风险溢价和 GARCH 动态都随状态变化的模型为

\(a_t=\sqrt{h_t}\epsilon_t\)，\(\epsilon_t\sim N(0,1)\)，\(\alpha_{ij}\) 满足无条件方差存在的条件。状态 \(s_t\) 是两状态马尔可夫链，转移概率

\(e_i\) 小表示倾向于停留在状态 \(i\)。停留期服从几何分布，期望持续期为 \(1/e_i\)。

推导拆解：刚进入状态 \(i\) 后，每期以概率 \(e_i\) 离开。恰好停留 \(k\) 期的概率是 \((1-e_i)^{k-1}e_i\)（前 \(k-1\) 期都留下，第 \(k\) 期离开）。期望 \(\sum_{k\ge1}k(1-e_i)^{k-1}e_i=e_i\cdot\frac{1}{e_i^2}=\frac1{e_i}\)，这里用了几何级数求导的公式 \(\sum_kkx^{k-1}=1/(1-x)^2\)。 数值例：月度数据 \(e_1=0.089\)，期望持续 \(1/0.089\approx11\) 个月。长期处于各状态的时间比例（平稳分布）是 \(\pi_1=e_2/(e_1+e_2)\)，\(\pi_2=e_1/(e_1+e_2)\)：从 1 流向 2 的"流量" \(\pi_1e_1\) 必须等于反向流量 \(\pi_2e_2\)。

• 识别约束 \(\beta_2>\beta_1\)：状态 2 的风险溢价更高。它只是为了给状态唯一地贴标签（否则两个状态互换后似然不变，称为标签切换，label switching）。
• \(\alpha_{1j}=\alpha_{2j}\) 时，所有状态共用一个 GARCH；若把 \(\beta_i\sqrt{h_t}\) 换成 \(\beta_i\)，就是简单的 Markov 转换 GARCH。(12.55) 是 Markov 转换 GARCH-M 模型。
• 收益依赖波动率意味着收益有序列相关、有一定可预测性；但未来状态未知，预测是各种状态配置的混合，点预测的不确定性很高。

设 \(h_1\) 固定为样本方差（也可以当参数估计，大样本下影响可忽略）。传统参数为 \(\beta=(\beta_1,\beta_2)'\)、\(\alpha_i=(\alpha_{i0},\alpha_{i1},\alpha_{i2})'\)、\(e=(e_1,e_2)'\)；状态向量 \(S=(s_1,\dots,s_n)'\) 是增广参数。给定 \(h_1\)、\(\alpha_i\) 与 \(S\)，\(H\) 可以递推算出。

12.10.2 Gibbs 步骤

似然是对 \(2^n\) 种状态配置的混合，而且因为 GARCH 的路径依赖，Hamilton 滤波不能直接用。Gibbs 只需要 \(f(\beta|R,S,H,\alpha_1,\alpha_2)\)、\(f(\alpha_i|R,S,H,\alpha_{j\ne i})\)、\(P(S|R,h_1,\alpha_1,\alpha_2)\) 和 \(f(e_i|S)\)。先验：\(\beta_i\sim N(\beta_{io},\sigma_{io}^2)\)，\(e_i\sim\mathrm{Beta}(\gamma_{i1},\gamma_{i2})\)，\(\alpha_{ij}\) 在适当区间上均匀分布（非线性参数，用 Griddy Gibbs）。

1. \(\beta_i\)：只依赖状态 \(i\) 的数据。令 \(r_{it}=r_t/\sqrt{h_t}\)（\(s_t=i\)），则 \(r_{it}=\beta_i+\epsilon_t\)。记 \(n_i\) 为状态 \(i\) 的数据点数，\(\bar r_i=\sum_{s_t=i}r_{it}/n_i\)，后验正态：\(1/\sigma_{i*}^2=n_i+1/\sigma_{io}^2\)，\(\beta_i^*=\sigma_{i*}^2(n_i\bar r_i+\beta_{io}/\sigma_{io}^2)\)。
2. \(\alpha_{ij}\)：逐个用 Griddy Gibbs。对数条件后验 \(\propto-\frac12\sum_{s_t=i}[\ln h_t+(r_t-\beta_i\sqrt{h_t})^2/h_t]\)，\(h_t\) 中含 \(\alpha_{ij}\)；在满足约束的区间上取格点（如 \(0\le\alpha_{11}<1-\alpha_{12}\)）。
3. \(e_i\)：\(\ell_1\) 为 1→2 的转换次数，\(\ell_2\) 为 2→1 的转换次数，由 Result 12.3，后验 Beta(\(\gamma_{i1}+\ell_i\), \(\gamma_{i2}+n_i-\ell_i\))。
4. \(s_j\)：逐个抽取。\(P(s_j|\cdot)\propto\prod_{t=j}^nf(a_t|H)\cdot P(s_j|S_{-j})\)，其中 \(P(s_j=i|S_{-j})=P(s_j=i|s_{j-1},s_{j+1})\propto P(s_j=i|s_{j-1})P(s_{j+1}|s_j=i)\) 由转移概率算出。设 \(s_j=i\) 后从 \(t=j\) 起递推 \(h_t\)，似然 \(L(s_j=i)\propto\exp(f_{ji})\)，\(f_{ji}=\sum_{t=j}^n-\frac12[\ln h_t+a_t^2/h_t]\)，\(a_t=r_t-\beta_{s_t}\sqrt{h_t}\)。于是

再用一个 \(U[0,1]\) 随机数决定 \(s_j\)。

白话解释：这又是"两个假设的贝叶斯定理"，结构与第 12a 章 (12.19) 的异常值概率完全相同。先验部分 \(P(s_j=i|s_{j-1},s_{j+1})\) 来自邻居：前后两天都是动荡，今天也很可能是动荡；似然部分 \(L(s_j=i)\) 来自数据：把今天设为状态 \(i\) 后，收益看起来有多"正常"。例如前一天平静、后一天平静，\(e_1=0.05\) 时先验上今天是动荡的可能性约为 \(\frac{0.05\times0.05}{0.95\times0.95+0.05\times0.05}\approx0.3\%\)（\(e_2\) 也取 0.05 时）；只有当天的收益在平静状态下极不寻常，似然比才可能把它翻过来。

备注：\(e_1,e_2\) 很小时，\(s_j\) 与 \(s_{j+1}\) 高度相关，联合抽取若干个 \(s_j\) 更有效，但要枚举的状态配置数随联合个数指数增长。另外，GARCH 的路径依赖意味着改变 \(s_j\) 会影响其后所有 \(h_t\)，所以似然要从 \(j\) 累积到 \(n\)，单次迭代的计算量是 \(O(n^2)\)。若各状态内没有 GARCH（如本章实战代码），\(s_j\) 的似然只涉及 \(a_j\) 一项，计算量降为 \(O(n)\)；此时也可以像 FFBS 那样对离散状态做前向滤波后向抽样（Chib 1996）。

12.10.3 例 12.6：GE 月度收益

数据为 GE 月度对数收益（%），1926-01 至 1999-12，888 个观测。

GARCH-M 基准：

p 值都小于 0.0006，Ljung–Box 检验没有问题，风险溢价为正且显著。这是一个 GARCH(1,2) 型的方差方程（\(h_t\) 依赖两期滞后）。令 \(\eta_t=a_t^2-h_t\)，把 \(h_t=a_t^2-\eta_t\) 代入，得到 \(a_t^2\) 的 ARMA 表示

AR 多项式可以分解为 \((1-0.945B)(1-0.820B)\)，两个实根的倒数模都小于 1，方差平稳。（精读笔记把右边的 MA 部分记为 \((1-0.025B)\eta_t\)，按上面的代入应为 \((1-1.740B+0.775B^2)\eta_t\)。）无条件方差为 \(0.546/(1-1.765+0.775)\)。注意分母只有约 0.01，对系数的舍入极其敏感：用表中三位小数算得约 55，原书给出 49.64（对应分母 0.011，即未舍入的系数）。

Markov 转换 GARCH-M：先验 \(\beta_1\sim N(0.3,0.09)\)，\(\beta_2\sim N(1.3,0.09)\)，\(e_i\sim\mathrm{Beta}(5,95)\)。初值 \(e_i=0.1\)，\(s_1\) 取等概率 Bernoulli，其余按初始转移概率依次生成，\(\alpha_1=(1.0,0.6,0.2)'\)，\(\alpha_2=(2,0.7,0.1)'\)。Griddy Gibbs 用 400 个格点，范围 \(\alpha_{i0}\in[0,6]\)，\(\alpha_{i1}\in[0,1]\)，\(\alpha_{i2}\in[0,0.5]\)，约束 \(\alpha_{i1}+\alpha_{i2}<1\)。共迭代 7000 次，用后 2000 次。

原书表 12.5 的后验均值（标准差）：

结论：

• 风险溢价的差异显著（\(\beta_2-\beta_1=0.135\)，在 5% 水平显著）。注意差值的后验是直接对每次抽样求差得到的，这正是第 12a 章 12.2.2 节说的便利。
• 波动率参数的后验均值差异不显著，但后验分布的形状不同（图 12.15、12.16）。图 12.17 显示状态 1 的持续性 \(\alpha_{11}+\alpha_{12}\) 常常触及上界 1.0，状态 2 不会。
• 两个状态的期望持续期约为 11 个月（\(1/0.089\)）和 9 个月（\(1/0.112\)）。图 12.14(b) 是每个观测处于状态 2 的后验概率。
• 图 12.18 显示两个模型拟合的波动率形态相似，都与平方收益一致；简单 GARCH-M 更平滑，估计的波动率也更低。

12.11 MCMC 预测

12.11.1 方法

在 MCMC 框架下做预测非常简单：在每次 Gibbs 迭代中，用这一次抽到的参数和状态模拟一条未来路径。以单变量 SV（12.20–12.21）为例，第 \(j\) 次迭代的参数为 \(\beta_j\)、\(\alpha_j\)、\(\sigma_{v,j}^2\)：

步骤：抽 \(v_{n+1}\sim N(0,\sigma_{v,j}^2)\)，由 (12.59) 算出 \(h_{n+1,j}\)；抽 \(\epsilon_{n+1}\sim N(0,1)\)，令 \(a_{n+1,j}=\sqrt{h_{n+1,j}}\,\epsilon_{n+1}\)（精读笔记指出原文漏了开方），由 (12.58) 算出 \(r_{n+1,j}\)；对 \(i=2,\dots,\ell\) 依次重复。解释变量需要已知或能依次预测。

运行 \(M+N\) 次迭代时，只需对后 \(N\) 次做预测，得到随机样本 \(\{r_{n+1,j},\dots,r_{n+\ell,j}\}_{j=1}^N\) 与 \(\{h_{n+1,j},\dots,h_{n+\ell,j}\}_{j=1}^N\)。样本均值是点预测，样本标准差衡量预测误差。用重要性抽样还可以提高波动率预测的效率（Gelman et al. 2003）。

这样得到的是预测分布（predictive distribution），它同时包含三种不确定性：未来冲击、当前潜在波动率 \(h_n\) 的不确定性、参数的不确定性。VaR 所需的分位数可以直接从样本中读出。

金融直觉：为什么纳入参数不确定性后 VaR 会变大？把每次迭代当作一个"情景"：有的情景里 \(\sigma_v\) 偏大、今天的 \(h_n\) 偏高，有的偏低。把所有情景的模拟收益合在一起，就是许多不同方差的分布的混合。混合分布的方差等于"各情景方差的平均 + 各情景均值的离散"，而且尾部比用平均参数算出的单一分布更厚。99% 分位数这种尾部指标对此最敏感，所以 VaR 上升。 这和信用风险里"PD 本身估不准，应在压力测试中考虑 PD 的区间"是同一种审慎原则；监管资本模型对短样本、新产品要求附加的模型风险缓冲，背后也是这个逻辑。MCMC 的好处是不必另外拍一个缓冲系数，参数不确定性由后验分布自动给出。

12.11.2 例 12.7：S&P 500 的波动率预测

接例 12.3，用 1962–1999 的 S&P 500 月度对数收益，预测原点为 1999-12。原书表 12.6，1–5 步：

GARCH 的波动率预测逐渐升向无条件方差 \(3.349/(1-0.086-0.735)\approx18.7\)（原书写作 18.78；这里的 GARCH 参数是该样本期的估计，与 (12.26) 不同）。SV 用 2000+2000 次迭代。

SV 的预测高于 GARCH。原书的解释是：SV 的 MCMC 预测纳入了参数不确定性，而 GARCH 把参数当作固定已知；这也是 GARCH 预测相对期权隐含波动率倾向于低估波动的原因之一。MCMC 给出的是波动率的整个预测分布，比点预测信息更多，可直接得到 VaR 所需的分位数。

12.11.3 其他应用

Zhang, Russell & Tsay (2008) 用 MCMC 分析买卖报价的信息决定因素；McCulloch & Tsay (2001) 估计 IBM 交易数据的分层模型；Eraker (2001)、Elerian, Chib & Shephard (2001) 估计扩散方程；VaR 计算中可以自然地评估预测分布。原书的结论是：关键问题不是 MCMC 能不能用，而是效率。

量化实战

应用场景

1. 波动率预测与仓位：SV 模型的后验波动率可用于波动率目标策略、风险平价的协方差输入。FFBS + KSC 混合是 SV 估计的标准算法，单变量模型在普通笔记本上几秒就能跑完。
2. VaR 与参数不确定性：MCMC 预测分布天然包含参数和潜在波动率的不确定性，得到的 VaR 比插入式（plug-in）估计更保守。对样本短、参数估计不稳的新品种尤其重要。
3. 期权：带杠杆的 SV 是 Heston 模型的离散版本；\(\rho\approx-0.4\) 的杠杆系数对应期权市场的负偏斜（skew）。用历史收益估计的 SV 参数可以作为期权定价模型的先验或校准起点。
4. 区制识别与风控开关：Markov 转换模型识别"平静/动荡"区制，状态概率可以作为降杠杆、切换风险预算或换用不同因子模型的开关。实盘只能用滤波概率 \(P(s_t|F_t)\)；Gibbs 后验概率和 Hamilton 平滑概率都用了全样本，放进回测就是前视。
5. 时变相关与对冲：多元 SV 的 \(q_{21,t}\) 就是时变 β（原书习题 12.6），可与第 11b 章 Kalman 滤波的时变 β 对照。

Python 示例一：FFBS–Gibbs 估计 SV 模型并计算 VaR

模拟 1000 个日收益（%），真值 \(\mu=0.03\)，\(\sigma_0=1.2\)，\(\alpha=0.95\)，\(\sigma_\eta=0.25\)。Gibbs 每次迭代依次抽：混合指示 \(I_t\)、用 FFBS 联合抽整条 \(z\)、\(\alpha\)、\(\sigma_\eta^2\)、\(\mu\)、\(\sigma_0^2\)，并在 burn-in 之后每次迭代模拟 50 条未来 5 天的路径。最后比较 MCMC 预测分布、插入式估计和 GARCH(1,1) 的次日 99% VaR。

import numpy as np
from arch import arch_model

rng = np.random.default_rng(2024)

# ================= 模拟 SV 数据：r_t = mu + sigma0*exp(z_t/2)*eps_t, z_{t+1} = alpha z_t + eta_t =================
n = 1000
mu_t, s0_t, a_t, se_t = 0.03, 1.2, 0.95, 0.25          # 日收益(%)
z = np.zeros(n); z[0] = rng.normal(0, se_t / np.sqrt(1 - a_t ** 2))
for t in range(n - 1):
    z[t + 1] = a_t * z[t] + rng.normal(0, se_t)
r = mu_t + s0_t * np.exp(z / 2) * rng.normal(size=n)

# KSC (1998) 七正态混合近似 ln(chi2_1)（原书表 12.3）
pk = np.array([0.00730, 0.10556, 0.00002, 0.04395, 0.34001, 0.24566, 0.25750])
mk = np.array([-11.4004, -5.2432, -9.8373, 1.5075, -0.6510, 0.5248, -2.3586])
vk = np.array([5.7960, 2.6137, 5.1795, 0.1674, 0.6401, 0.3402, 1.2626])

def ffbs(y, c, H, alpha, s2eta):
    """线性高斯模型 y_t = c_t + z_t + e_t, z_{t+1}=alpha z_t + eta_t 的前向滤波后向抽样 (12.45)-(12.50)"""
    n_ = len(y)
    zf, Pf = np.empty(n_), np.empty(n_)
    zp, Pp = 0.0, s2eta / (1 - alpha ** 2)              # z_1 取平稳分布
    for t in range(n_):
        V = Pp + H[t]
        K = Pp / V
        zf[t] = zp + K * (y[t] - c[t] - zp)
        Pf[t] = Pp * (1 - K)
        zp, Pp = alpha * zf[t], alpha ** 2 * Pf[t] + s2eta
    zs = np.empty(n_)
    zs[-1] = rng.normal(zf[-1], np.sqrt(Pf[-1]))
    for t in range(n_ - 2, -1, -1):                     # p(z_t | z_{t+1}, F_t)
        Pnext = alpha ** 2 * Pf[t] + s2eta
        m = zf[t] + alpha * Pf[t] / Pnext * (zs[t + 1] - alpha * zf[t])
        v = Pf[t] - alpha ** 2 * Pf[t] ** 2 / Pnext
        zs[t] = rng.normal(m, np.sqrt(v))
    return zs

# 先验：mu~N(0,1)；alpha~N(0.9,0.1^2)；m*lam/sigma0^2 ~ chi2_m；m*lam/s2eta ~ chi2_m
def sv_gibbs(r, n_iter=3000, burn=1000, horizon=5, n_paths=50):
    n_ = len(r)
    mu, s02, alpha, s2eta = r.mean(), r.var(), 0.9, 0.1
    zz = np.zeros(n_)
    out, fc_r, fc_h, vol = [], [], [], np.zeros(n_)
    for it in range(n_iter):
        # (1) 混合指示 I_t：后验 p_i * N(y_t; z_t+m_i, v_i)（原书写作 CDF，应为密度）
        y = np.log((r - mu) ** 2 / s02 + 1e-8)
        resid = y[:, None] - zz[:, None] - mk[None, :]
        logw = np.log(pk) - 0.5 * np.log(vk) - 0.5 * resid ** 2 / vk
        w = np.exp(logw - logw.max(1, keepdims=True)); w /= w.sum(1, keepdims=True)
        I = (w.cumsum(1) > rng.random((n_, 1))).argmax(1)
        # (2) FFBS 联合抽取整条对数波动路径 z
        zz = ffbs(y, mk[I], vk[I], alpha, s2eta)
        # (3) alpha | z：AR(1) 系数，正态先验
        prec = zz[:-1] @ zz[:-1] / s2eta + 1 / 0.1 ** 2
        mean = (zz[:-1] @ zz[1:] / s2eta + 0.9 / 0.1 ** 2) / prec
        alpha = np.clip(rng.normal(mean, np.sqrt(1 / prec)), -0.999, 0.999)
        # (4) s2eta | z, alpha：逆卡方
        e = zz[1:] - alpha * zz[:-1]
        s2eta = (5 * 0.01 + e @ e) / rng.chisquare(5 + n_ - 1)
        # (5) mu | 其余：异方差回归 (12.22)
        sd = np.sqrt(s02) * np.exp(zz / 2)
        prec = np.sum(1 / sd ** 2) + 1.0
        mu = rng.normal(np.sum(r / sd ** 2) / prec, np.sqrt(1 / prec))
        # (6) sigma0^2 | 其余：v_t=(r_t-mu)exp(-z_t/2) iid N(0,sigma0^2)
        vv = (r - mu) * np.exp(-zz / 2)
        s02 = (5 * 1.0 + vv @ vv) / rng.chisquare(5 + n_)
        if it >= burn:
            out.append([mu, np.sqrt(s02), alpha, np.sqrt(s2eta)])
            vol += np.sqrt(s02) * np.exp(zz / 2)
            # (7) 预测：用本次迭代的参数模拟 n_paths 条未来路径 (12.58)(12.59)
            zf_ = np.full(n_paths, zz[-1]); rr, hh = [], []
            for h in range(horizon):
                zf_ = alpha * zf_ + rng.normal(0, np.sqrt(s2eta), n_paths)
                sig = np.sqrt(s02) * np.exp(zf_ / 2)
                hh.append(sig ** 2); rr.append(mu + sig * rng.normal(size=n_paths))
            fc_r.append(np.array(rr).T); fc_h.append(np.array(hh).T)
    return (np.array(out), np.concatenate(fc_r), np.concatenate(fc_h), vol / (n_iter - burn))

draws, fc_r, fc_h, vol_post = sv_gibbs(r)
print("参数        真值    后验均值  后验标准差")
for k, (nm, tv) in enumerate(zip(["mu", "sigma0", "alpha", "sigma_eta"], [mu_t, s0_t, a_t, se_t])):
    print(f"{nm:10s} {tv:6.3f}   {draws[:, k].mean():7.3f}   {draws[:, k].std():7.3f}")
true_vol = s0_t * np.exp(z / 2)
g = arch_model(r, mean="Constant", vol="GARCH", p=1, q=1).fit(disp="off")
print(f"真实波动率与估计的相关系数: SV 后验均值 {np.corrcoef(true_vol, vol_post)[0,1]:.3f}, "
      f"GARCH(1,1) {np.corrcoef(true_vol, g.conditional_volatility)[0,1]:.3f}")

# 预测分布与 VaR
print("\n未来 1-5 天 波动率预测(后验均值 sqrt(E h)):", np.round(np.sqrt(fc_h.mean(0)), 3))
gf = g.forecast(horizon=5)
print("GARCH(1,1) 波动率预测                     :", np.round(np.sqrt(gf.variance.values[-1]), 3))
q01 = np.quantile(fc_r[:, 0], 0.01)
# 插入式(plug-in)：参数固定在后验均值、最后一天 z 也固定在后验均值，只模拟未来冲击
m_, s0_, a_, se_ = draws.mean(0)
z_last = 2 * np.log(vol_post[-1] / s0_)
sim = m_ + s0_ * np.exp((a_ * z_last + rng.normal(0, se_, 200000)) / 2) * rng.normal(size=200000)
print(f"次日 99% VaR(收益率%)：MCMC 预测分布 {-q01:.3f} | 插入式(忽略参数与状态不确定性) {-np.quantile(sim, 0.01):.3f}"
      f" | GARCH 正态 {-(gf.mean.values[-1,0] - 2.326 * np.sqrt(gf.variance.values[-1,0])):.3f}")
print(f"100 万元多头次日 99% VaR ≈ {1e6 * -q01 / 100:,.0f} 元")

关键输出：

参数        真值    后验均值  后验标准差
mu          0.030    -0.024     0.037
sigma0      1.200     1.260     0.093
alpha       0.950     0.944     0.017
sigma_eta   0.250     0.264     0.038
真实波动率与估计的相关系数: SV 后验均值 0.833, GARCH(1,1) 0.682

未来 1-5 天 波动率预测(后验均值 sqrt(E h)): [1.763 1.755 1.748 1.739 1.73 ]
GARCH(1,1) 波动率预测                     : [1.801 1.796 1.791 1.786 1.781]
次日 99% VaR(收益率%)：MCMC 预测分布 4.405 | 插入式(忽略参数与状态不确定性) 4.087 | GARCH 正态 4.215
100 万元多头次日 99% VaR ≈ 44,051 元

读输出：

• 3000 次迭代（前 1000 次 burn-in）在几秒内完成，\(\sigma_0\)、\(\alpha\)、\(\sigma_\eta\) 的后验均值都在真值的一个后验标准差左右。\(\mu\) 的后验均值为负而真值为正，但差距在 1.5 个后验标准差以内：日收益的均值本来就很难估。FFBS 一次抽出整条 \(z\)，即使 \(\alpha=0.95\) 这样的高持续性也能很快混合。
• SV 后验均值波动率与真实波动率的相关系数 0.83，高于 GARCH(1,1) 的 0.68。但要注意这不是公平比较：SV 的后验均值用了全样本（相当于平滑），GARCH 的条件波动率只用了过去数据（相当于滤波）。公平的比较应该在每个时点只用当时为止的数据重新估计（练习 8）。同时，数据本来就是按 SV 生成的，GARCH 在这里是错设模型。
• 次日 99% VaR：MCMC 预测分布 4.41%，插入式 4.09%，前者约高 8%。插入式把参数和最后一天的波动率都固定在后验均值，只模拟未来冲击，低估了尾部；这与例 12.7 中"SV 预测高于 GARCH，是因为纳入了参数不确定性"的逻辑一致。100 万元多头的次日 99% VaR 约 4.4 万元（原书习题 12.5 的做法）。

Python 示例二：Markov 转换模型的 Gibbs 抽样、Hamilton 滤波与风控开关

模拟 1500 个日收益，两个状态：平静（\(\mu=0.06\%\)，\(\sigma=0.8\%\)）与动荡（\(\mu=-0.10\%\)，\(\sigma=2.0\%\)），转移概率 \(e_0=0.02\)、\(e_1=0.05\)（期望持续期 50 天与 20 天）。先按 12.10.2 节的步骤做 Gibbs 抽样（各状态内无 GARCH，状态用奇偶分批的单点抽样，转移概率先验与例 12.6 相同，为 Beta(5,95)），再用 statsmodels 的 MarkovRegression 做极大似然（Hamilton 滤波），比较滤波概率与平滑概率，最后做一个"动荡概率高就空仓"的风控开关。

import numpy as np
import statsmodels.api as sm

rng = np.random.default_rng(99)

# ============ 模拟两状态 Markov 转换收益：状态 0 平静、状态 1 动荡 ============
n = 1500
mu_true, sd_true = np.array([0.06, -0.10]), np.array([0.8, 2.0])
e_true = np.array([0.02, 0.05])                        # e0=P(0->1), e1=P(1->0)
s = np.zeros(n, int)
for t in range(1, n):
    s[t] = 1 - s[t - 1] if rng.random() < e_true[s[t - 1]] else s[t - 1]
r = mu_true[s] + sd_true[s] * rng.normal(size=n)

# ============ Gibbs 抽样（原书 12.9 节思路，去掉 GARCH 部分）============
def ms_gibbs(r, n_iter=3000, burn=1000):
    n_ = len(r)
    mu, sig2, e = np.array([0.0, 0.0]), np.array([0.5, 4.0]), np.array([0.1, 0.1])
    st = (np.abs(r) > 2 * r.std()).astype(int)          # 初始状态：大波动日记为动荡
    keep, Pst = [], np.zeros(n_)
    idx_even, idx_odd = np.arange(0, n_, 2), np.arange(1, n_, 2)
    for it in range(n_iter):
        # (1) 状态 s_j：给定邻居 s_{j-1}, s_{j+1}，奇偶位置条件独立 -> 两批向量化更新
        Ptr = np.array([[1 - e[0], e[0]], [e[1], 1 - e[1]]])
        loglik = -0.5 * np.log(sig2)[None, :] - 0.5 * (r[:, None] - mu[None, :]) ** 2 / sig2[None, :]
        for idx in (idx_even, idx_odd):
            lp = loglik[idx].copy()
            has_prev, has_next = idx > 0, idx < n_ - 1
            lp[has_prev] += np.log(Ptr[st[idx[has_prev] - 1]])          # P(s_j | s_{j-1})
            lp[has_next] += np.log(Ptr[:, st[idx[has_next] + 1]].T)      # P(s_{j+1} | s_j)
            p1 = 1 / (1 + np.exp(lp[:, 0] - lp[:, 1]))
            st[idx] = (rng.random(len(idx)) < p1).astype(int)
        # (2) 各状态的均值与方差（共轭：正态 + 逆卡方）
        for i in (0, 1):
            ri = r[st == i]; ni = len(ri)
            prec = ni / sig2[i] + 1 / 1.0
            mu[i] = rng.normal((ri.sum() / sig2[i]) / prec, np.sqrt(1 / prec))
            sig2[i] = (3 * 1.0 + np.sum((ri - mu[i]) ** 2)) / rng.chisquare(3 + ni)
        if sig2[0] > sig2[1]:                           # 识别约束：状态 1 方差更大
            mu, sig2, st = mu[::-1].copy(), sig2[::-1].copy(), 1 - st
        # (3) 转移概率：Beta 先验 + 转换次数 (Result 12.3)
        for i in (0, 1):
            prev = st[:-1] == i
            l_i = np.sum(st[1:][prev] != i); n_i = prev.sum()
            e[i] = rng.beta(5 + l_i, 95 + n_i - l_i)
        if it >= burn:
            keep.append(np.r_[mu, np.sqrt(sig2), e]); Pst += st
    return np.array(keep), Pst / (n_iter - burn)

D, P_smooth_gibbs = ms_gibbs(r)
names = ["mu0", "mu1", "sd0", "sd1", "e0", "e1"]
truth = np.r_[mu_true, sd_true, e_true]
print("参数   真值    Gibbs 后验均值(标准差)")
for k, nm in enumerate(names):
    print(f"{nm:4s} {truth[k]:6.3f}   {D[:, k].mean():7.3f} ({D[:, k].std():.3f})")
dur = 1 / D[:, 4:6]
print("期望持续期(天) 后验均值:", np.round(dur.mean(0), 1), " 真值:", 1 / e_true,
      " 本样本实现值:", np.round([np.sum(s[:-1] == i) / np.sum((s[:-1] == i) & (s[1:] != i)) for i in (0, 1)], 1))

# ============ 极大似然（Hamilton 滤波）与滤波 / 平滑概率 ============
mod = sm.tsa.MarkovRegression(r, k_regimes=2, trend="c", switching_variance=True)
np.random.seed(0)                                      # search_reps 的随机初值用全局随机数
res = mod.fit(search_reps=20)
prm = dict(zip(mod.param_names, res.params))
sd_mle = np.sqrt([prm["sigma2[0]"], prm["sigma2[1]"]])
hi = int(np.argmax(sd_mle))                            # 识别哪个是动荡状态
p_filt = res.filtered_marginal_probabilities[:, hi]
p_smth = res.smoothed_marginal_probabilities[:, hi]
order = np.argsort(sd_mle)                             # 按 平静, 动荡 顺序输出
print("\nMLE 各状态标准差:", np.round(sd_mle[order], 3),
      " 期望持续期:", np.round(np.asarray(res.expected_durations)[order], 1))
acc = lambda p: np.mean((p > 0.5) == (s == 1))
print(f"状态识别准确率: 滤波概率 {acc(p_filt):.3f} | 平滑概率 {acc(p_smth):.3f} | Gibbs 后验 {acc(P_smooth_gibbs):.3f}")

# ============ 风控开关：动荡概率高时降仓 ============
def sharpe(x): return x.mean() / x.std() * np.sqrt(252)
pos_f = np.r_[1.0, (p_filt[:-1] < 0.5).astype(float)]   # t 日仓位只用 t-1 日收盘的滤波概率
pos_s = np.r_[1.0, (p_smth[:-1] < 0.5).astype(float)]   # 用平滑概率 = 偷看未来
print(f"年化Sharpe: 买入持有 {sharpe(r):.2f} | 滤波概率开关 {sharpe(pos_f * r):.2f} | "
      f"平滑概率开关(前视!) {sharpe(pos_s * r):.2f}")
print(f"最大单日亏损: 买入持有 {r.min():.2f}% | 滤波概率开关 {(pos_f * r).min():.2f}%")

关键输出：

参数   真值    Gibbs 后验均值(标准差)
mu0   0.060     0.068 (0.025)
mu1  -0.100    -0.059 (0.109)
sd0   0.800     0.787 (0.019)
sd1   2.000     2.044 (0.087)
e0    0.020     0.028 (0.006)
e1    0.050     0.071 (0.016)
期望持续期(天) 后验均值: [37.9 14.9]  真值: [50. 20.]  本样本实现值: [52.6 18.8]

MLE 各状态标准差: [0.785 2.057]  期望持续期: [37.5 12.2]
状态识别准确率: 滤波概率 0.895 | 平滑概率 0.938 | Gibbs 后验 0.942
年化Sharpe: 买入持有 0.47 | 滤波概率开关 0.63 | 平滑概率开关(前视!) 0.99
最大单日亏损: 买入持有 -5.22% | 滤波概率开关 -4.96%

读输出：

• Gibbs 后验和 MLE 都准确恢复了两个状态的标准差（0.79/2.04 对真值 0.8/2.0）。动荡状态的均值后验标准差达 0.11，远大于其绝对值：区制间波动率的差异很容易识别，均值（风险溢价）的差异很难识别，所以例 12.6 中 \(\beta_2-\beta_1\) 显著是一个有分量的结论。
• 两种方法估计的期望持续期（约 38/15 天与 38/12 天）都短于本样本实现的 53/19 天。原因是平静期里偶尔出现的大波动日被判为短暂的动荡，使转换次数偏多；Beta(5,95) 先验（均值 0.05）也把 \(e_0\) 往上拉。期望持续期是区制模型中最不稳的量，用于择时之前要做稳健性检验。
• 状态识别准确率：滤波概率 0.895，平滑概率 0.938，Gibbs 后验 0.942。后两者用了未来数据。
• 风控开关：用 \(t-1\) 日的滤波概率决定 \(t\) 日是否持仓，年化 Sharpe 从 0.47 提高到 0.63；用平滑概率则"提高"到 0.99，多出来的部分全是前视偏差。最大单日亏损只从 −5.22% 降到 −4.96%：区制切换的第一天总是来不及反应，滤波概率要看到几天大波动后才会确认动荡。

本章小结

SV 模型让对数波动率服从带独立冲击的 AR(1)，似然要对整条波动率路径积分，而 Gibbs 抽样把波动率当作增广参数，就把问题拆成了加权回归、AR(1) 回归和对 \(h_t\) 的抽样。逐点抽 \(h_t\) 的条件后验 (12.23) 可以用 Griddy Gibbs 或 Metropolis 处理，但混合慢；把观测取对数平方后，用 KSC 七正态混合近似 \(\ln\chi^2_1\)，模型就成了条件线性高斯状态空间模型，可以用 FFBS 一次抽出整条路径，杠杆效应则通过把 \(\eta_t\) 分解并线性化转移方程来处理。多元 SV 用 Cholesky 分解，\(q_{21,t}\) 给出时变相关和时变 β。Markov 转换 GARCH-M 模型把每期的状态当作参数逐个抽取，转移概率用 Beta 共轭后验，GE 的例子显示两个区制的风险溢价显著不同。MCMC 预测在每次迭代中模拟未来路径，得到同时包含参数不确定性的预测分布，可以直接算 VaR；实证上，它给出的波动率预测高于把参数当作已知的 GARCH。

练习

基础

1. 由 (12.21) 推导 \(\ln h_t\) 的无条件均值与方差，并用例 12.3 的后验均值计算 \(E(\ln h_t)\) 与对应的典型月度波动率。 提示：\(E\ln h_t=\alpha_0/(1-\alpha_1)=0.454/0.163\approx2.79\)，\(e^{2.79}\approx16.2\)，月度标准差约 4%。
2. 验证 (12.23) 推导中的配方：取 \(A=1\)，\(C=\alpha_1^2\)，\(a=\alpha_0+\alpha_1\ln h_{t-1}\)，\(b=(\ln h_{t+1}-\alpha_0)/\alpha_1\)，算出 \(c\) 等于 \(\mu_t\)。
3. 证明 \(\epsilon_t^*=\ln\epsilon_t^2\) 的均值约为 −1.27、方差为 \(\pi^2/2\)，并用表 12.3 验证混合分布的均值与方差（应很接近）。
4. 对 FFBS，推导 (12.51) 的 \(\mu_t^*\) 与 \(\Sigma_t^*\)。 提示：\((z_t,z_{t+1})|F_t\) 的协方差为 \(\alpha\Sigma_{t|t}\)，套用定理 11.1。
5. 例 12.6 中 \(e_1=0.089\)、\(e_2=0.112\)。写出两状态链的平稳分布，并计算长期处于状态 2 的时间比例。 提示：\(\pi_2=e_1/(e_1+e_2)\approx0.44\)。

进阶

6. 在示例一的代码中加入杠杆效应：用 \(\rho=-0.5\) 模拟数据，先用无杠杆的 FFBS–Gibbs 估计，观察哪些参数有偏；再按 (12.52)–(12.53) 的线性化思路修改前向滤波（提示：此时 FFBS 只是近似）。
7. （原书习题 12.5）对一条日收益序列（可用示例一的模拟数据代替 Cisco 2001–2008）拟合 SV 模型，得到 1 步波动率的预测分布，并计算 100 万美元多头、概率 0.01 的次日 VaR。再把 Gibbs 迭代数减半、加倍，看 VaR 的蒙特卡罗误差有多大。
8. 修改示例一，做一个公平的波动率预测比较：在最后 200 天，每 20 天用截至当时的数据重新跑 SV 的 Gibbs 和 GARCH 的 MLE，只用滤波值（SV 用 \(z_{t|t}\) 的后验，或在 FFBS 的前向步取 \(z_{t|t}\)）预测次日波动，比较 QLIKE 损失。
9. 在示例二中把转移概率先验改为 Beta(1,1)（均匀），比较期望持续期估计的变化；再把各状态内加入 GARCH(1,1)，体会 12.10.2 节所说的 \(O(n^2)\) 计算量。
10. （原书习题 12.4、12.6、12.7）Ford 月度收益的 GARCH 与 SV 比较；Ford 与 S&P 500 的二元 SV，讨论波动关系并计算 Ford 的时变 β；P&G 与价值加权指数的二元 SV 与 BEKK(1,1) 比较。

原书推荐习题：12.5（SV + 预测分布 + VaR，完整串起 12.8 与 12.11，最贴近风控实务）、12.4 与 12.7（GARCH/BEKK 与 SV 的比较）、12.6（二元 SV 求时变 β）。

原书对照

（原书印刷页码约等于 PDF 页码减 20。）