第 12a 章 MCMC 方法与贝叶斯推断

本章对应 Tsay 原书第 12 章 12.1–12.6 节。前面各章的估计几乎都是极大似然。可是一旦模型里有大量看不见的量——缺失的数据点、可能是异常值的观测、每天的潜在波动率、市场所处的区制——似然就要对这些量积分，常常算不出来。马尔可夫链蒙特卡罗（Markov chain Monte Carlo, MCMC）换了一个思路：把看不见的量也当作"参数"，从它们的联合后验分布中抽样，再用样本做推断。本章讲 MCMC 的基本工具（Gibbs 抽样、共轭先验、Metropolis–Hastings、Griddy Gibbs），以及三个基础应用：带时间序列误差的回归、缺失值插补、加性异常值检测。第 12b 章用这些工具估计随机波动率模型和 Markov 转换模型。

学习目标

1. 说明 MCMC 的基本思想：构造一个平稳分布等于后验分布的马尔可夫链；理解数据增广把潜变量当作参数的做法。
2. 会写 Gibbs 抽样的步骤，知道 burn-in、分块抽样和收敛诊断的作用。
3. 熟记常用共轭先验的后验公式（正态均值、正态精度、Beta–Bernoulli、Gamma–Poisson、逆卡方–方差），理解"后验精度 = 先验精度 + 数据精度"。
4. 会用 Metropolis、Metropolis–Hastings 和 Griddy Gibbs 处理没有闭式条件后验的参数。
5. 能独立推导并编程实现带 AR 误差回归的 Gibbs 抽样，以及 AR 模型中缺失值和加性异常值的条件后验。

读前导读

这一章在解决什么问题。 你熟悉的估计方法是"写出似然、求最大值"。它要求似然能算出来。可是当模型里有一大堆看不见的东西——哪些报价是坏点、停牌那几天的真实价格、每天的潜在波动率——似然就要对所有这些未知量求积分，几百维、几千维，根本算不动。MCMC 换了个办法：不求最大值，而是生成一大批参数的随机样本，让这批样本的分布恰好等于后验分布。有了样本，均值、标准差、分位数、任何函数的分布都只是简单统计。你在 CFA 里见过的蒙特卡罗模拟（用随机路径给期权或 VaR 定价）是"已知分布、抽样算数"；MCMC 是"分布只知道形状、不会直接抽样，就设计一条随机游走去逼近它"。

贝叶斯的部分你可以这样理解：先验是"看数据之前的判断"，数据带来新证据，后验是两者按可信度加权的结果。这和 Black–Litterman 把市场均衡收益和投资者观点按精度加权是同一回事，也像审计里结合内控评估（先验）和实质性测试结果（数据）形成结论。

需要先想起来的数学。

• 条件概率与贝叶斯定理。\(P(A|B)=P(B|A)P(A)/P(B)\)。分布版本 \(f(\theta|X)\propto f(X|\theta)P(\theta)\) 里的 \(\propto\) 读作"正比于"，意思是差一个与 \(\theta\) 无关的常数，最后靠"总概率为 1"补回来。见 第 00 册第 07 章 概率中的分析工具。
• 积分作为"求总和"。\(f(X)=\int f(X|\theta)P(\theta)d\theta\) 是把所有可能的 \(\theta\) 按先验加权求和，离散情形就是 \(\sum_\theta\)。它只是归一化常数，MCMC 的妙处正在于不用算它。见 第 00 册第 03 章 积分。
• 配方。\(a\mu^2-2b\mu+c=a(\mu-b/a)^2+\text{常数}\)。正态密度的指数是二次函数，凡是"指数部分对 \(\mu\) 是二次"的分布，配方后就能读出均值 \(b/a\) 和方差 \(1/a\)。共轭先验的推导几乎都靠这一招。
• 矩阵求逆与精度矩阵。多元情形下"精度"是协方差矩阵的逆 \(\Sigma^{-1}\)。公式 \(\Sigma_*^{-1}=\Sigma_o^{-1}+n\Sigma^{-1}\) 就是标量 \(1/\sigma_*^2=1/\sigma_o^2+n/\sigma^2\) 的矩阵版本。见 第 00 册第 06 章 线性代数速成。
• Gamma 函数与几个分布。\(\Gamma(\alpha)\) 是阶乘的推广，\(\Gamma(n)=(n-1)!\)；它只出现在归一化常数里，读公式时可以忽略。Beta 分布描述 0 到 1 之间的概率（如胜率），Gamma 分布描述正数（如精度、强度），卡方分布 \(\chi^2_v\) 是 \(v\) 个独立标准正态的平方和。

怎么读这一章。 必读：12.1.2（MCMC 的核心想法）、12.2（Gibbs）、12.3 中的 Result 12.1、12.1a、12.3、12.8（其余共轭结果可以当查表用）、12.4.1–12.4.2（Metropolis 和 MH）、12.5（第一个完整的 Gibbs 推导，是后面所有推导的模板）。12.1.3 的 EM 历史、Result 12.4–12.7、Griddy Gibbs 第一次可以只看结论。12.6–12.7 推导较细，建议先读讲解框弄懂"把缺失值当作回归系数"这个技巧，再看公式。最后一定跑一遍或细读量化实战的例一，代码与公式一一对应，是检验理解的最好方式。

12.1 为什么需要 MCMC

12.1.1 马尔可夫过程

先回顾一个概念。随机过程 \(\{X_t\}\) 取值于空间 \(\Omega\)，如果给定现在 \(X_t\)，未来 \(X_h\)（\(h>t\)）与过去 \(X_s\)（\(s<t\)）无关：

就称为马尔可夫过程（Markov process）。离散时间下写作 \(P(X_h|X_t,X_{t-1},\dots)=P(X_h|X_t)\)。转移概率函数 \(P_t(\theta,h,A)=P(X_h\in A|X_t=\theta)\)；如果它只依赖 \(h-t\) 而不依赖 \(t\)，称链有平稳转移分布。离散状态马尔可夫链的基本性质见第 02 册。

12.1.2 马尔可夫链模拟

贝叶斯推断要的是后验分布 \(P(\theta|X)\)。MCMC 的想法是：在参数空间 \(\Theta\) 上构造一个马尔可夫链，让它的平稳分布恰好是 \(P(\theta|X)\)，然后让链跑足够长时间，当前值的分布就足够接近后验。之后链上的每一个值都可以看作（相关的）后验抽样。对同一个目标分布，可以构造出很多满足条件的链，这些方法统称 MCMC。

白话解释："平稳分布"可以用信用评级迁移矩阵来理解。一年期迁移矩阵反复作用很多年后，一家公司处在各评级的概率会稳定下来，而且与它今天的初始评级无关，这个极限就是平稳分布。MCMC 反过来做：先指定想要的极限（后验），再设计转移规则，使得链跑久了各处被访问的频率正好等于后验概率。 两个代价要记住。第一，链的开头还记得初值，不能用（这就是后面的 burn-in）。第二，相邻样本是相关的，2000 个 MCMC 样本的信息量通常少于 2000 个独立样本，"有效样本量"衡量的就是这个折扣。

12.1.3 从 EM 算法到数据增广

MCMC 在统计中的流行，有一条从 EM 算法（Dempster, Laird & Rubin 1977）过来的脉络。EM 处理缺失数据：

• M 步：假装缺失值已知，用完整数据的方法做极大似然估计；
• E 步：给定数据和当前拟合的模型，求缺失值的条件期望并填进去；

从任意初值开始反复迭代直到收敛。Tanner & Wong (1987) 做了两方面推广：

1. 迭代模拟：把"填条件期望"换成"从条件分布中随机抽一个值"；
2. 数据增广（data augmentation）：主动加入一些辅助变量，常常能让条件分布变简单、模拟变快。

后面会反复看到数据增广的威力：异常值模型里每个时点加一个"是不是异常"的 0/1 变量，随机波动率模型里把每天的波动率当作参数，Markov 转换模型里把每天的状态当作参数。加进去之后，每个条件分布都是熟悉的分布。

12.2 Gibbs 抽样

12.2.1 算法

Gibbs 抽样（Gibbs sampling；Geman & Geman 1984；Gelfand & Smith 1990）是最流行的 MCMC 方法。这里"参数"的含义很宽：缺失的数据点可以当作参数；一个交易日内 \(N\) 笔成交背后的 \(N\) 个"真实价格"也可以当作 \(N\) 个参数。

设有三个参数 \(\theta_1,\theta_2,\theta_3\)，数据 \(X\)，模型 \(M\)。似然可能很难写，但三个完全条件分布（full conditional distributions）

都能抽样。不必知道它们的确切形式，只需能从中抽随机数。

给定任意初值 \(\theta_{2,0},\theta_{3,0}\)：

1. 从 \(f_1(\theta_1|\theta_{2,0},\theta_{3,0},X,M)\) 抽 \(\theta_{1,1}\)；
2. 从 \(f_2(\theta_2|\theta_{3,0},\theta_{1,1},X,M)\) 抽 \(\theta_{2,1}\)；
3. 从 \(f_3(\theta_3|\theta_{1,1},\theta_{2,1},X,M)\) 抽 \(\theta_{3,1}\)。

这就完成了一次 Gibbs 迭代。注意每一步都用最新的值作条件。重复 \(m\) 次得到 \((\theta_{1,m},\theta_{2,m},\theta_{3,m})\)。在正则条件下（实质是要求从任意初值出发都能到达整个参数空间，理论见 Tierney 1994），\(m\) 足够大时它近似是联合后验 \(f(\theta_1,\theta_2,\theta_3|X,M)\) 的一次抽样。

为什么平稳分布是后验？ 设 \((\theta_1,\theta_2,\theta_3)\) 已经服从联合后验 \(f\)。第 1 步只替换 \(\theta_1\)，而且是按 \(f(\theta_1|\theta_2,\theta_3)\) 替换，所以替换后 \((\theta_1,\theta_2,\theta_3)\) 的联合分布仍是 \(f(\theta_2,\theta_3)f(\theta_1|\theta_2,\theta_3)=f\)。每一步都保持 \(f\) 不变，整个迭代当然也保持 \(f\) 不变。

推导拆解：把这句话展开成三步。

1. 假设迭代前 \((\theta_1,\theta_2,\theta_3)\sim f\)，那么其中 \((\theta_2,\theta_3)\) 的边际分布是 \(f(\theta_2,\theta_3)\)。
2. 第 1 步丢掉旧的 \(\theta_1\)，按 \(f(\theta_1|\theta_2,\theta_3)\) 重新抽一个。新三元组的密度 \(=\)（\(\theta_2,\theta_3\) 的密度）\(\times\)（给定它们时新 \(\theta_1\) 的密度）\(=f(\theta_2,\theta_3)f(\theta_1|\theta_2,\theta_3)\)。
3. 由条件密度的定义 \(f(\theta_1|\theta_2,\theta_3)=f(\theta_1,\theta_2,\theta_3)/f(\theta_2,\theta_3)\)，乘积正好是 \(f(\theta_1,\theta_2,\theta_3)\)。 这只说明"已经到了后验就不会离开"。"从任意初值出发最终会到后验"需要额外的正则条件，正文括号里提到的就是它。

白话解释：Gibbs 像轮流调整几个相互关联的估计。比如同时估计一只股票的 β 和特质波动率：先假定波动率已知，回归抽一个 β；再用这个 β 算残差，抽一个波动率；如此循环。每一步都只解一个简单的一元或低维问题，合起来却在探索整个联合分布。

12.2.2 实践要点

实践中运行 \(n\) 次迭代，丢弃前 \(m\) 次，得到 Gibbs 样本

点估计与后验方差：

要检验 \(H_0:\theta_1=\theta_2\)，就对每次抽样算差值 \(\theta_{1,j}-\theta_{2,j}\)，再看它的均值与标准差，用 \(t=\hat\theta/\hat\sigma\)。任何参数的函数（比如两个状态的期望持续期之差、某个分位数）都可以这样直接得到后验，这是 MCMC 相对 delta 方法的一大便利。

• burn-in：被丢弃的前 \(m\) 个样本叫 burn-in，用来让链"忘掉"初值。另一种做法是用不同初值跑很多条短链，各取最后一次抽样组成样本。
• 分块：Gibbs 把高维问题拆成一串低维条件抽样，极端情况下是 \(N\) 个一元分布。但如果某些参数高度相关，逐个抽样会收敛得很慢（链在狭长的后验山脊上只能走小碎步），应把它们放在一起联合抽样，如用 \(f(\theta_1,\theta_2|\theta_3)\) 和 \(f_3(\theta_3|\theta_1,\theta_2)\)（Liu, Wong & Kong 1994）。第 12b 章的 FFBS 就是"把整条波动率路径一次抽出来"的分块方法。
• 收敛诊断：理论只保证 \(m\) 足够大时收敛，没有给出具体的 \(m\)。诊断方法很多但没有共识，没有任何方法能 100% 保证收敛（Carlin & Louis 2000；Gelman et al. 2003）。原书的建议是用不同初值重复多次，看结论是否一致。常用的补充工具有：迹图（trace plot）、样本自相关与有效样本量（effective sample size, ESS）、多链的 Gelman–Rubin \(\hat R\) 统计量（接近 1 表示各链混合良好）。

12.3 贝叶斯推断与共轭先验

12.3.1 后验分布

统计推断有两种范式：经典方法（极大似然）和贝叶斯方法（先验 + 数据 → 后验）。本书前面用的都是经典方法，但每个问题都有贝叶斯解法，借助 MCMC 已经可以计算，而且多数情况下两种方法结果相近。有时贝叶斯方法更有优势，例如算 VaR 时可以自然地纳入参数不确定性（第 12b 章例 12.7），代价是计算量大。

设先验 \(P(\theta)\)，似然 \(f(X|\theta)\)，由贝叶斯定理

分母 \(f(X)\) 只是归一化常数。从 (12.5) 看，基于似然的推断等价于取常数先验的贝叶斯推断。完全条件分布在贝叶斯文献中叫条件后验分布（conditional posterior distribution）。

12.3.2 共轭先验

如果先验和后验属于同一个分布族，就称先验是共轭先验（conjugate prior）。在 MCMC 中，这意味着条件后验有闭式，可以直接调用常规随机数生成器（DeGroot 1970 第 9 章）。下面是金融计量中最常用的几个结果。

Result 12.1（正态均值，方差已知）。\(x_1,\dots,x_n\) iid \(N(\mu,\sigma^2)\)，\(\sigma^2\) 已知，先验 \(\mu\sim N(\mu_o,\sigma_o^2)\)，则后验 \(\mu|X\sim N(\mu_*,\sigma_*^2)\)：

用精度（precision，方差的倒数）写更好记。令 \(\eta=1/\sigma^2\)，\(\eta_o=1/\sigma_o^2\)，\(\eta_*=1/\sigma_*^2\)：

后验精度 = 先验精度 + 数据精度（\(\bar x\) 是充分统计量，精度为 \(n\eta\)）；后验均值是先验均值与样本均值按精度的加权平均。\(n\) 增大时先验的影响逐渐消失。

证明思路：\(f(\mu|X)\propto\exp[-\frac{n(\bar x-\mu)^2}{2\sigma^2}-\frac{(\mu-\mu_o)^2}{2\sigma_o^2}]\)，指数是 \(\mu\) 的二次函数，配方即得。

推导拆解：配方的具体过程。

1. 似然 \(\prod_i\exp[-(x_i-\mu)^2/(2\sigma^2)]\) 中，\(\sum(x_i-\mu)^2=\sum(x_i-\bar x)^2+n(\bar x-\mu)^2\)，第一项与 \(\mu\) 无关，并入常数。这就是"\(\bar x\) 是充分统计量"的含义。
2. 用精度写指数：\(-\frac12[n\eta(\mu-\bar x)^2+\eta_o(\mu-\mu_o)^2]\)。展开只看含 \(\mu\) 的项：\(-\frac12[(n\eta+\eta_o)\mu^2-2(n\eta\bar x+\eta_o\mu_o)\mu]\)。
3. 套"\(a\mu^2-2b\mu\) 的中心是 \(b/a\)、精度是 \(a\)"：精度 \(\eta_*=n\eta+\eta_o\)，均值 \(\mu_*=(n\eta\bar x+\eta_o\mu_o)/\eta_*\)。 数值例：先验认为某基金经理的月 α 均值 0、标准差 0.5%（精度 4）；36 个月样本均值 0.6%，月残差标准差 3%，样本均值的标准误 \(3/\sqrt{36}=0.5\%\)（数据精度也是 4）。后验均值 \(=(4\times0+4\times0.6)/8=0.3\%\)，后验标准差 \(1/\sqrt8\approx0.35\%\)。先验和数据一样可信，于是各占一半，结果把 α 收缩了一半。

Result 12.1a（多元版本）。\(x_i\) iid \(N(\mu,\Sigma)\)，\(\Sigma\) 已知，先验 \(\mu\sim N(\mu_o,\Sigma_o)\)，则后验 \(N(\mu_*,\Sigma_*)\)：

把它用在回归上（把 OLS 估计 \(\hat\beta\sim N(\beta,\sigma^2(X'X)^{-1})\) 当作"一个观测"），就得到回归系数的条件后验，这是后面所有回归型 Gibbs 步骤的核心（Box & Tiao 1973）。

Result 12.2（正态精度，均值已知）。Gamma 分布 \(f(\eta|\alpha,\beta)=\frac{\beta^\alpha}{\Gamma(\alpha)}\eta^{\alpha-1}e^{-\beta\eta}\)，均值 \(\alpha/\beta\)，方差 \(\alpha/\beta^2\)。若 \(x_i\) iid \(N(\mu,1/\eta)\)，\(\mu\) 已知，先验 \(\eta\sim\Gamma(\alpha,\beta)\)，则后验

Result 12.3（Bernoulli 概率）。Beta 分布 \(f(\theta|\alpha,\beta)=\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}(1-\theta)^{\beta-1}\)，均值 \(\alpha/(\alpha+\beta)\)，方差 \(\frac{\alpha\beta}{(\alpha+\beta)^2(\alpha+\beta+1)}\)。若 \(x_i\) iid Bernoulli(\(\theta\))，先验 Beta(\(\alpha,\beta\))，则后验

先验参数可以理解为"事先看到了 \(\alpha\) 次成功、\(\beta\) 次失败"。异常值模型中的异常比例 \(\epsilon\)、Markov 转换模型中的转移概率 \(e_i\) 都用这个结果。

Result 12.4（Poisson 强度）。\(x_i\) iid Poisson(\(\lambda\))，先验 \(\lambda\sim\Gamma(\alpha,\beta)\)，则后验 \(\Gamma(\alpha+\sum x_i,\ \beta+n)\)。（比如每分钟成交笔数的强度。）

Result 12.5（指数分布）。\(x_i\) iid 指数分布（率 \(\lambda\)），先验 \(\Gamma(\alpha,\beta)\)，则后验 \(\Gamma(\alpha+n,\ \beta+\sum x_i)\)。（比如交易间隔时间，第 05 章的久期。）

Result 12.6（负二项）。负二项分布 \(p(n|m,\lambda)=\binom{m+n-1}{n}\lambda^m(1-\lambda)^n\)。原书的例子：公司有 \(m\) 个职位，每个 MBA 应聘者独立地以概率 \(\lambda\) 合适，面试总人数为 \(Y\)，则 \(X=Y-m\) 服从负二项。\(m\) 固定、先验 \(\lambda\sim\) Beta(\(\alpha,\beta\)) 时，后验为 Beta(\(\alpha+mn\), \(\beta+\sum x_i\))。

Result 12.7（正态，均值与精度都未知）。正态–Gamma 先验：\(\mu|\eta\sim N(\mu_o,(\tau_o\eta)^{-1})\)，\(\eta\sim\Gamma(\alpha,\beta)\)。则 \(\mu|\eta,X\) 为正态，均值 \(\mu_*=\frac{\tau_o\mu_o+n\bar x}{\tau_o+n}\)，精度 \((\tau_o+n)\eta\)；\(\eta\) 的边际后验为 \(\Gamma(\alpha+n/2,\beta_*)\)，

最后一项惩罚样本均值与先验均值的分歧。

Result 12.8（零均值正态的方差，逆卡方先验）。若 \(1/Y\sim\chi^2_v\)，称 \(Y\) 服从逆卡方分布，密度 \(f(y|v)=\frac{2^{-v/2}}{\Gamma(v/2)}y^{-(v/2+1)}e^{-1/(2y)}\)，\(E(Y)=1/(v-2)\)（\(v>2\)），\(\mathrm{Var}(Y)=\frac{2}{(v-2)^2(v-4)}\)（\(v>4\)）。若 \(a_1,\dots,a_n\) iid \(N(0,\sigma^2)\)，先验 \(v\lambda/\sigma^2\sim\chi^2_v\)，则后验

抽样时就是 \(\sigma^2=(v\lambda+\sum a_i^2)/\chi^2_{v+n}\)。先验可以理解为"事先看过 \(v\) 个观测，其平均平方是 \(\lambda\)"。它和 Result 12.2 是同一件事的两种写法。

12.4 没有闭式时：替代算法

条件后验没有闭式时，需要别的抽样方法。

12.4.1 Metropolis 算法

适用于条件后验只知道到差一个归一化常数的情形（Metropolis & Ulam 1949；Metropolis et al. 1953）。

1. 取初值 \(\theta_0\)，满足 \(f(\theta_0|X)>0\)；
2. 对 \(t=1,2,\dots\)： (a) 从跳跃分布或建议分布（jumping / proposal distribution）\(J_t(\theta^*|\theta_{t-1})\) 中抽候选 \(\theta^*\)。\(J_t\) 必须对称：\(J_t(\theta_i|\theta_j)=J_t(\theta_j|\theta_i)\)； (b) 计算 \(r=f(\theta^*|X)/f(\theta_{t-1}|X)\)； (c) 以概率 \(\min(r,1)\) 令 \(\theta_t=\theta^*\)，否则 \(\theta_t=\theta_{t-1}\)。

正则条件下 \(\{\theta_t\}\) 依分布收敛到 \(f(\theta|X)\)。只用到比值，所以不需要归一化常数。需要做的只有三件事：算 \(r\)、从跳跃分布抽样、抽一个均匀随机数决定是否接受。直观上：往后验密度更高的地方跳一定接受；往更低的地方跳以概率 \(r\) 接受，这样链不会只停在众数，而会按密度比例访问各处。最常用的对称跳跃分布是以当前值为中心的正态或 Student-t，叫随机游走 Metropolis（random-walk Metropolis）。

白话解释：一个小数值例。当前值 \(\theta_{t-1}\) 处后验密度（未归一化）为 10，候选 \(\theta^*\) 处为 4，则 \(r=0.4\)。抽一个 \(U(0,1)\) 随机数，小于 0.4 就跳过去，否则原地不动，并把原值再记录一次。长期下来，链停在第一个点附近的时间约是第二个点的 2.5 倍，正好与密度之比一致。 归一化常数为什么能消掉：\(f(\theta|X)=f(X|\theta)P(\theta)/f(X)\)，分子分母的 \(f(X)\) 相同，比值里直接约掉。所以只要能算"似然 × 先验"就够了。 实务上建议分布的步长要调：太小几乎每次都接受，但每步只挪一点，样本高度相关；太大几乎总被拒绝，链长时间不动。一维问题接受率在 40% 左右、高维问题在 25% 左右通常比较合适。

12.4.2 Metropolis–Hastings 算法

Hastings (1970) 去掉了对称的要求，接受比改为

分子分母各除以"被建议的概率"，抵消建议分布本身的偏向。

为什么正确：设转移核为 \(K(\theta\to\theta^*)=J(\theta^*|\theta)\min\{1,r\}\)。可以验证细致平衡（detailed balance）\(f(\theta)K(\theta\to\theta^*)=f(\theta^*)K(\theta^*\to\theta)\)：两边都等于 \(\min\{f(\theta)J(\theta^*|\theta),\ f(\theta^*)J(\theta|\theta^*)\}\)。满足细致平衡的链以 \(f\) 为平稳分布。提高 MH 效率的方法见 Tierney (1994)。

推导拆解：为什么两边都等于那个 \(\min\)。

1. 记 \(a=f(\theta)J(\theta^*|\theta)\)，\(b=f(\theta^*)J(\theta|\theta^*)\)，则从 \(\theta\) 出发的接受比是 \(r=b/a\)。
2. 左边 \(f(\theta)K(\theta\to\theta^*)=a\min\{1,b/a\}=\min\{a,b\}\)（把 \(a\) 乘进 \(\min\) 的两项）。
3. 反方向的接受比是 \(a/b\)，右边 \(=b\min\{1,a/b\}=\min\{b,a\}\)。两边相等。
4. 细致平衡推出平稳：对 \(\theta\) 积分，\(\int f(\theta)K(\theta\to\theta^*)d\theta=f(\theta^*)\int K(\theta^*\to\theta)d\theta=f(\theta^*)\)，即"一步之后流入 \(\theta^*\) 的总概率 = 原来在 \(\theta^*\) 的概率"。（这里略去了"拒绝后留在原地"那部分概率，把它加回来结论不变。） 直观上，细致平衡说的是任意两点之间"来往的资金流"两两相抵，整体分布自然不再变化。

12.4.3 Griddy Gibbs

金融模型里常有非线性参数，例如 ARMA 的 MA 系数、GARCH 的系数，它们的条件后验没有闭式。Tanner (1996) 的 Griddy Gibbs 适用于一元条件后验，适用面广，但可能低效。对标量参数 \(\theta_i\)（\(\theta_{-i}\) 为其余参数）：

1. 在一个适当的区间内取格点 \(\theta_{i1}\le\theta_{i2}\le\cdots\le\theta_{im}\)，计算 \(w_j=f(\theta_{ij}|X,\theta_{-i})\)（只需到差一个常数）；
2. 用 \(\{w_j\}\) 近似 \(\theta_i\) 的逆累积分布函数；
3. 抽 \(U\sim U(0,1)\)，经近似逆 CDF 变换得到 \(\theta_i\) 的一次抽样。

最简单的近似是把 \(\theta_i\) 当作取值于格点的离散分布 \(p(\theta_{ij})=w_j/\sum_vw_v\)。区间的选择需要检查：如果 Gibbs 样本的直方图在端点处有明显的概率质量，说明区间太窄，要扩大；如果概率集中在区间中间一小块，说明太宽（大多数 \(w_j\approx0\)，浪费计算），要缩小。

Griddy Gibbs 或 MH 可以嵌入 Gibbs 抽样里，只用来抽部分参数，其余参数仍用闭式条件后验，这叫 Metropolis-within-Gibbs。第 12b 章的随机波动率和 Markov 转换 GARCH 都是这种混合结构。

12.5 带时间序列误差的线性回归

12.5.1 模型与思路

原书第 2 章用 SCA 估计过这类模型：

\(x_t=(1,x_{1t},\dots,x_{kt})'\)，可以包含 \(y_t\) 的滞后。参数 \(\theta=(\beta',\phi,\sigma^2)'\)。

思路很自然：已知时间序列模型，回归就好估；已知回归，残差 \(z_t=y_t-x_t'\beta\) 的 AR 模型也好估。Gibbs 抽样正是在两者之间来回迭代，只是每次都"抽"而不是"估"。需要三个条件后验 \(f(\beta|Y,X,\phi,\sigma^2)\)、\(f(\phi|Y,X,\beta,\sigma^2)\)、\(f(\sigma^2|Y,X,\beta,\phi)\)。

独立的共轭先验：

先验中的常数叫超参数（hyperparameters），通常取 \(\beta_o=0\)，\(\phi_o=0\)，\(\Sigma_o\) 为大对角阵，即弱信息先验。

12.5.2 三个条件后验

\(\beta\) 的条件后验。给定 \(\phi\)，做准差分（第 05 册第 16 章的 Cochrane–Orcutt 变换）：\(y_{o,t}=y_t-\phi y_{t-1}\)，\(x_{o,t}=x_t-\phi x_{t-1}\)，得到误差为白噪声的回归

OLS 估计 \(\hat\beta=(\sum x_{o,t}x_{o,t}')^{-1}\sum x_{o,t}y_{o,t}\)，分布为 \(N(\beta,\sigma^2(\sum x_{o,t}x_{o,t}')^{-1})\)。由 Result 12.1a：

推导拆解：这一步把三件熟悉的事串了起来。

1. 准差分：\(y_t=x_t'\beta+z_t\) 减去 \(\phi\) 倍的 \(y_{t-1}=x_{t-1}'\beta+z_{t-1}\)，得 \(y_t-\phi y_{t-1}=(x_t-\phi x_{t-1})'\beta+(z_t-\phi z_{t-1})\)，最后一项正是白噪声 \(a_t\)。这在 CFA 里讲序列相关的修正时出现过。
2. OLS 的分布：误差是白噪声后，\(\hat\beta\sim N(\beta,\sigma^2(X_o'X_o)^{-1})\)，数据精度就是 \(X_o'X_o/\sigma^2\)（这里 \(X_o'X_o=\sum x_{o,t}x_{o,t}'\)）。
3. 套 Result 12.1a：把 \(\hat\beta\) 当作"一个观测"，后验精度 = 先验精度 \(\Sigma_o^{-1}\) + 数据精度；后验均值 = 两个均值按精度矩阵加权。 先验很弱（\(\Sigma_o\) 很大，\(\Sigma_o^{-1}\approx0\)）时，\(\beta_*\approx\hat\beta\)，\(\Sigma_*\approx\sigma^2(X_o'X_o)^{-1}\)，正好回到 GLS 的结果。

\(\phi\) 的条件后验。给定 \(\beta\)，\(z_t=y_t-\beta'x_t\) 已知，\(\hat\phi=\sum z_{t-1}z_t/\sum z_{t-1}^2\)，方差 \(\sigma^2/\sum z_{t-1}^2\)。由 Result 12.1：

\(\sigma^2\) 的条件后验。\(a_t=z_t-\phi z_{t-1}\) 已知，由 Result 12.8：

Gibbs 步骤：(1) 设定超参数；(2) 给任意初值，比如忽略序列相关的 OLS；(3) 从 (12.9) 抽 \(\beta\)；(4) 从 (12.10) 抽 \(\phi\)；(5) 从 (12.11) 抽 \(\sigma^2\)；重复 (3)–(5) 多次，丢掉 burn-in 后用样本均值作点估计。AR(\(p\)) 误差的推广（原书习题 12.2）只需把准差分改为 \(\phi(B)\)，把 (12.10) 换成多元的 Result 12.1a。

12.5.3 例 12.1：3 年期与 1 年期国债利率

数据为美国 1 年期与 3 年期国债恒定期限周利率，1962-01-05 至 2009-04-10（圣路易斯联储），样本量 2466。利率有单位根，所以用周变化 \(c_{3t}=r_{3t}-r_{3,t-1}\)、\(c_{1t}=r_{1t}-r_{1,t-1}\)（%）。原书第 2 章用 MA(1) 误差，这里用 AR(2) 误差。传统方法（R）得到

标准误分别为 0.0075、0.0201、0.0201。除滞后 4、6 处残差 ACF 边际显著外，模型充分。

贝叶斯设定 \(c_{3t}=\beta c_{1t}+z_t\)，\(z_t=\phi_1z_{t-1}+\phi_2z_{t-2}+a_t\)（12.13），先验 \(\beta\sim N(0,4)\)，\(\phi\sim N[0,\mathrm{diag}(0.25,0.16)]\)，\((v\lambda)/\sigma^2=(10\times0.05)/\sigma^2\sim\chi^2_{10}\)。初值取两步 OLS 估计，迭代 2100 次，丢弃前 100 次。后验均值（后验标准差）为：

\(\sigma\approx0.069\)。原书图 12.1 的迹图稳定，图 12.2 是边际后验直方图；换初值结果相似，判断已经收敛。后验均值与 (12.12) 很接近，因为样本大、模型简单，先验几乎不起作用。

12.6 缺失值

12.6.1 缺失值与异常值是同一个问题

加性异常值（additive outlier, AO）的定义是

\(\omega\) 是异常幅度，\(x_t\) 是没有异常的序列。录入错误、坏的报价、测量误差都会造成 AO。即使只有几个，AO 也可能严重扭曲参数估计，甚至导致模型误设。

判断 \(y_h\) 是不是 AO 的思路是：先假装 \(x_h\) 缺失，求它在其余数据下的条件分布；如果观测到的 \(y_h\) 在这个分布下很可能出现，就不是异常；如果概率很小，就判为 AO。所以缺失值处理和 AO 检测基于同一个思想。缺失值可以用 Kalman 滤波处理（第 11a 章 11.7 节、第 11b 章 11.15 节；Jones 1980），也可以用 MCMC（McCulloch & Tsay 1994a）。异常值文献有 Chang, Tiao & Chen (1988)、Tsay (1988)、Tsay, Peña & Pankratz (2000)；异常值按影响方式分四类，这里只讨论 AO。

12.6.2 AR(\(p\)) 中单个缺失值的条件后验

设 \(x_t=\phi_1x_{t-1}+\cdots+\phi_px_{t-p}+a_t\)（12.15），\(x_h\) 缺失（\(1<h<n\)）。把 \(x_h\) 当作未知参数，参数向量 \(\theta=(\phi',x_h,\sigma^2)'\)。先验 \(\phi\sim N(\phi_o,\Sigma_o)\)，\(x_h\sim N(\mu_o,\sigma_o^2)\)，\(v\lambda/\sigma^2\sim\chi^2_v\)。\(\phi\) 与 \(\sigma^2\) 的条件后验同 12.5 节；关键是 \(x_h\) 的条件后验。

给定模型和数据，\(x_h\) 只与它前后各 \(p\) 个邻居 \(\{x_{h-p},\dots,x_{h-1},x_{h+1},\dots,x_{h+p}\}\) 有关，因为 \(x_h\) 只出现在 \(t=h,\dots,h+p\) 这 \(p+1\) 个方程里。把这些方程改写成"关于 \(x_h\) 的回归"：

1. \(t=h\)：令 \(y_h=\phi_1x_{h-1}+\cdots+\phi_px_{h-p}\)，\(b_h=-a_h\)，则 \(y_h=x_h+b_h=\phi_0x_h+b_h\)，\(\phi_0=1\)；
2. \(t=h+j\)（\(j=1,\dots,p\)）：把方程中除 \(x_h\) 以外的项移到左边，令 \(y_{h+j}=x_{h+j}-\phi_1x_{h+j-1}-\cdots-\phi_{j-1}x_{h+1}-\phi_{j+1}x_{h-1}-\cdots-\phi_px_{h+j-p}\)，\(b_{h+j}=a_{h+j}\)，则 \(y_{h+j}=\phi_jx_h+b_{h+j}\)。

于是得到 \(p+1\) 个方程

（由正态对称性，\(-a_h\) 与 \(a_h\) 同分布）这是一个只有 \(p+1\) 个数据点、没有截距的简单线性回归。LS 估计与方差为

\(p=1\) 时

叫 \(x_h\) 的滤波值（filtered value）。前后两个邻居权重相同，因为高斯 AR(1) 是时间可逆的。

推导拆解：\(p=1\) 的情形逐步写出来。

1. \(x_h\) 出现在两个方程里：\(x_h=\phi_1x_{h-1}+a_h\) 和 \(x_{h+1}=\phi_1x_h+a_{h+1}\)。
2. 把它们看成"以 \(x_h\) 为未知系数"的回归：第一式 \(y_h:=\phi_1x_{h-1}=1\cdot x_h-a_h\)；第二式 \(y_{h+1}:=x_{h+1}=\phi_1\cdot x_h+a_{h+1}\)。回归元分别是 \(1\) 和 \(\phi_1\)。
3. 无截距回归的 OLS 是 \(\sum(\text{回归元}\times\text{因变量})/\sum\text{回归元}^2=(\phi_1x_{h-1}+\phi_1x_{h+1})/(1+\phi_1^2)\)。 数值例：\(\phi_1=0.5\)，前一天 \(x_{h-1}=2\)，后一天 \(x_{h+1}=1\)，则 \(\hat x_h=0.4\times3=1.2\)。注意它不是简单的前后平均 1.5，而是向均值 0 收缩：\(\phi_1\) 越小，序列越"健忘"，邻居提供的信息越少，估计越靠近无条件均值。 金融上，这就是停牌一天的股票（用去均值的 AR(1) 近似其超额收益或价差）的合理插补值：同时参考停牌前后的数据，而不是简单沿用前一天。

再由 Result 12.1，\(x_h\) 的条件后验为正态 \(N(\mu_*,\sigma_*^2)\)：

12.6.3 连续缺失

若 \(x_h,x_{h+1}\) 都缺失，它们与 \(\{x_{h-p},\dots,x_{h-1};x_{h+2},\dots,x_{h+p+1}\}\) 相关。有两种做法：

1. 直接推广：构造一个以 \((x_h,x_{h+1})\) 为系数的多元线性回归，结合先验得二元正态后验，在 Gibbs 中联合抽取（原书习题 12.3）；
2. 在一次 Gibbs 迭代中多次使用单个缺失值的公式，逐个抽取。

因为相邻缺失值高度相关，联合抽取更好，缺失段较长时尤其如此；缺失点少时逐个抽取也可以。以上推导假设 \(h-p\ge1\)、\(h+p\le n\)；靠近样本两端时，要相应减少回归中的数据点。

12.7 加性异常值检测

12.7.1 模型与数据增广

在 MCMC 框架下，AO 检测很直接。McCulloch & Tsay (1994a) 的简单 Gibbs 抽样效果良好，除非存在幅度相近的成片 AO（Justel, Peña & Tsay 2001）。模型为

• \(\delta_t\) 是独立的 Bernoulli 变量，\(P(\delta_t=1)=\epsilon\)，表示第 \(t\) 期是否出现 AO；
• \(\beta_t\) 是独立的异常幅度；
• \(x_t=\phi_0+\phi_1x_{t-1}+\cdots+\phi_px_{t-p}+a_t\) 是没有异常的 AR(\(p\))。

每一期都可能出现 AO，概率为 \(\epsilon\)。\(n\) 个数据，却有 \(2n+p+3\) 个参数：\(\phi\)（\(p+1\) 个）、\(\delta\)、\(\beta\)、\(\sigma^2\)、\(\epsilon\)。\(\delta\) 和 \(\beta\) 是数据增广引入的。

共轭先验：\(\phi\sim N(\phi_o,\Sigma_o)\)，\(v\lambda/\sigma^2\sim\chi^2_v\)，\(\epsilon\sim\mathrm{Beta}(\gamma_1,\gamma_2)\)，\(\beta_t\sim N(0,\xi^2)\)。需要的条件后验有五个。

12.7.2 条件后验

\(\phi\) 与 \(\sigma^2\)。给定 \(\delta,\beta\)，去掉异常得 \(x_t=y_t-\delta_t\beta_t\)，就是普通的 AR(\(p\)) 回归：\(\hat\phi=(\sum_{t=p+1}^n\mathbf x_{t-1}\mathbf x_{t-1}')^{-1}\sum\mathbf x_{t-1}x_t\)，\(\mathbf x_{t-1}=(1,x_{t-1},\dots,x_{t-p})'\)，后验与 (12.9) 同形。\(\sigma^2\) 的后验为 \((v\lambda+\sum_{t=p+1}^na_t^2)/\sigma^2\sim\chi^2_{v+(n-p)}\)，\(a_t=x_t-\phi'\mathbf x_{t-1}\)。

\(\epsilon\)。给定 \(\delta\)，由 Result 12.3，\(\epsilon\sim\mathrm{Beta}(\gamma_1+\sum\delta_t,\ \gamma_2+n-\sum\delta_t)\)。（原书没有单列这一步，但它是完整 Gibbs 循环的一部分。）

\(\delta_h\)。它只和 \(j=h-p,\dots,h+p\) 的 \(y_j,\beta_j,\delta_j\)（\(j\ne h\)）以及 \(\phi,\sigma^2\) 有关。令 \(x_j^*=x_j\)（\(j\ne h\)），\(x_h^*=y_h\)（第 \(h\) 点暂且取观测值本身），定义

• 情形 I（\(\delta_h=0\)，先验概率 \(1-\epsilon\)）：\(x_h^*=x_h\)，所以 \(w_j=a_j\sim N(0,\sigma^2)\)。
• 情形 II（\(\delta_h=1\)，先验概率 \(\epsilon\)）：\(x_h^*=x_h+\beta_h\)，于是 \(w_h=a_h+\beta_h\sim N(\beta_h,\sigma^2)\)，而 \(j>h\) 时 \(w_j=a_j-\phi_{j-h}\beta_h\sim N(-\phi_{j-h}\beta_h,\sigma^2)\)。令 \(\psi_0=-1\)，\(\psi_i=\phi_i\)，统一写成 \(w_j\sim N(-\psi_{j-h}\beta_h,\sigma^2)\)。

令 \(m=\min(n,h+p)\)，比较两种情形下的似然并按先验概率加权：

白话解释：(12.19) 是贝叶斯定理在"两个假设"上的直接应用：

\[P(\text{是异常}|\text{数据})=\frac{P(\text{是异常})\times P(\text{数据}|\text{是异常})}{P(\text{是异常})P(\text{数据}|\text{是异常})+P(\text{不是})P(\text{数据}|\text{不是})}.\]

先验概率就是 \(\epsilon\) 与 \(1-\epsilon\)；两个 \(\exp[\cdot]\) 是两种情形下局部残差 \(w_j\) 的正态似然（共同的常数因子约掉了）。如果把 \(y_h\) 当真值时，它和前后几天的残差都很大，而扣掉幅度 \(\beta_h\) 后残差变小，"是异常"的似然就压倒性地大。 金融类比：交易对账时看到一笔离谱的成交价，你会同时看它前后几笔是否也跟着动。若后面的价格没有跟随，更像录入错误；若后面也延续了，更像真实的价格变动。\(w_j\)（\(j>h\)）正是在做"后面有没有跟随"的检查。 代码里写成 \(1/(1+\frac{1-\epsilon}{\epsilon}e^{l_0-l_1})\) 是同一个式子，分子分母同除以分子，避免指数溢出。

\(\beta_h\)。若 \(\delta_h=0\)，数据不含 \(\beta_h\) 的信息，从先验 \(N(0,\xi^2)\) 抽。若 \(\delta_h=1\)，由回归 \(w_j=-\psi_{j-h}\beta_h+a_j\) 得 \(\hat\beta_h=\frac{-\sum_{j=h}^m\psi_{j-h}w_j}{\sum_{j=h}^m\psi_{j-h}^2}\)，方差 \(\sigma^2/\sum\psi_{j-h}^2\)，结合先验得

所有 Gibbs 迭代结束后，\(\delta_h\) 的样本均值就是第 \(h\) 期为 AO 的后验概率。

12.7.3 例 12.2：3 年期国债周变化

数据为美国 3 年期国债周变化 \(c_{3t}\)，1988-03-18 至 1999-09-10，600 个观测（例 12.1 的子序列）。PACF 建议 AR(3)，ML 拟合为

标准误 0.041、0.042、0.041，\(Q(12)=11.4\) 不显著。（笔记中第三项的下标误印为 \(t-2\)，此处已更正为 \(t-3\)。）

Gibbs 联合估计 + AO 检测的设定：先验 \(\phi\sim N(0,0.25I_3)\)，\(v\lambda/\sigma^2=5\times0.00256/\sigma^2\sim\chi^2_5\)，\(\gamma_1=5\)，\(\gamma_2=95\)（预期约 5% 的点是 AO），\(\xi^2=0.1\)。这里 \(0.00256\approx\sigma^2/5\)，\(\xi^2\approx9\sigma^2\)，即先验认为异常幅度的标准差约为 3 倍噪声标准差。初值 \(\epsilon=0.05\)，\(\sigma^2=0.012\)，\(\phi=(0.2,0.02,0.1)\)；迭代 1050 次，丢弃前 50 次。

结果：

后验标准差 0.046、0.045、0.046、0.0008，与 ML 结果相近。原书图 12.3(b) 画出各点为 AO 的后验概率，(c) 为后验异常幅度。两个最突出的点：

• \(t=323\)（1994-05-20）：AO 后验概率 0.83，幅度 −0.304。\(c_{3t}\) 从 0.24 变为 −0.34，两周内利率约下降 0.6 个百分点。
• \(t=201\)（1992-01-17）：AO 后验概率 0.58，幅度 0.176，\(c_{3t}\) 从 −0.02 变为 0.33。

备注：Gibbs 检测的计算量大，但它联合估计参数和异常值。传统方法把估计与检测分开做，速度快，但存在多个异常值时可能产生虚假检测（一个异常值扭曲了参数，使另一个正常点看起来异常）。用 SCA 软件的传统方法也识别出 \(t=323\) 和 \(t=201\) 为最显著的两个 AO，幅度为 −0.39 与 0.36。

量化实战

应用场景

1. 贝叶斯回归与收缩：12.5 节的 Gibbs 框架可以直接推广到因子模型的贝叶斯估计。先验 \(\beta\sim N(\beta_o,\Sigma_o)\) 就是一种收缩（shrinkage）：小样本或高维因子回归中，把系数向先验均值收缩能显著降低估计误差。Black–Litterman 模型（把市场均衡收益作为先验，把投资者观点作为"数据"）用的正是 Result 12.1a 这种精度加权。
2. 行情数据清洗：12.7 节的 AO 模型可以用来识别坏报价、录入错误和闪崩中的离群成交。与"超过 3 倍标准差就剔除"的规则相比，它利用了序列的动态结构，并且同时估计模型参数，不会因为异常值本身而把阈值撑大。
3. 缺失值插补：停牌、节假日不对齐时，12.6 节给出 AR 模型下缺失值的条件分布，可用于多重插补（multiple imputation）——不是只填一个值，而是抽多个值，让下游的风险估计反映插补的不确定性。
4. 参数不确定性：任何 MCMC 输出都是后验样本。把每组参数抽样代入组合优化或风险计算，就得到"参数不确定性下"的风险分布，比单点估计更稳健（第 12b 章的 VaR 例子）。

Python 示例：带 AR(2) 误差回归的 Gibbs 抽样与 AO 检测

例一仿照例 12.1：模拟 2466 周的利率变化，\(c_{3t}=0.782c_{1t}+z_t\)，\(z_t\) 为 AR(2)（0.183, −0.036），\(\sigma_a=0.068\)，先验与原书相同；与 statsmodels 的 MLE 比较，并用一条从极端初值出发的链检查收敛。例二仿照例 12.2：模拟 600 个 AR(3) 观测，在 \(t=201\) 加 \(+0.45\)、在 \(t=323\) 加 \(-0.50\) 两个 AO，用 McCulloch–Tsay 的 Gibbs 抽样联合估计参数与 AO。

import numpy as np
import statsmodels.api as sm

rng = np.random.default_rng(12)

# ============ 例一：带 AR(2) 误差的回归的 Gibbs 抽样（仿例 12.1）============
n = 2466
c1 = rng.normal(0, 0.10, n)                          # "1 年期利率周变化"
z = np.zeros(n); a = rng.normal(0, 0.068, n)
for t in range(2, n):
    z[t] = 0.183 * z[t - 1] - 0.036 * z[t - 2] + a[t]
c3 = 0.782 * c1 + z                                   # "3 年期利率周变化"

# 先验 (同例 12.1)：beta~N(0,4)，phi~N(0,diag(.25,.16))，(v*lam)/sigma2 ~ chi2_v，v=10, lam=0.05
b0, B0 = 0.0, 4.0
phi0, Phi0 = np.zeros(2), np.diag([0.25, 0.16])
v, lam = 10, 0.05
p = 2

def gibbs_reg_ar(y, x, n_iter=2100, burn=100, init=None):
    if init is None:                                  # 初值：两步 OLS
        beta = np.sum(x * y) / np.sum(x * x)
        zz = y - beta * x
        Zl = np.column_stack([zz[1:-1], zz[:-2]])
        phi = np.linalg.lstsq(Zl, zz[2:], rcond=None)[0]
        sig2 = np.var(zz[2:] - Zl @ phi)
    else:
        beta, phi, sig2 = init
    draws = []
    for it in range(n_iter):
        # (1) beta | phi, sigma2：准差分后为普通回归 (12.8)(12.9)
        yo = y[p:] - phi[0] * y[p - 1:-1] - phi[1] * y[:-p]
        xo = x[p:] - phi[0] * x[p - 1:-1] - phi[1] * x[:-p]
        prec = np.sum(xo * xo) / sig2 + 1 / B0
        mean = (np.sum(xo * yo) / sig2 + b0 / B0) / prec
        beta = rng.normal(mean, np.sqrt(1 / prec))
        # (2) phi | beta, sigma2：z_t 的 AR(2) 回归 + 正态先验 (12.10)
        zz = y - beta * x
        Zl = np.column_stack([zz[1:-1], zz[:-2]]); zt = zz[2:]
        Sinv = Zl.T @ Zl / sig2 + np.linalg.inv(Phi0)
        S = np.linalg.inv(Sinv)
        m = S @ (Zl.T @ zt / sig2 + np.linalg.inv(Phi0) @ phi0)
        phi = rng.multivariate_normal(m, S)
        # (3) sigma2 | beta, phi：逆卡方 (12.11)
        res = zt - Zl @ phi
        sig2 = (v * lam + np.sum(res ** 2)) / rng.chisquare(v + len(res))
        draws.append([beta, phi[0], phi[1], sig2])
    return np.array(draws)[burn:]

D = gibbs_reg_ar(c3, c1)
names = ["beta", "phi1", "phi2", "sigma2"]
print("Gibbs 后验均值(后验标准差):")
for k, nm in enumerate(names):
    print(f"  {nm:7s} {D[:, k].mean():8.4f} ({D[:, k].std():.4f})")

def ess(x):                                           # 有效样本量：用自相关和截到首个负值
    x = x - x.mean(); n_ = len(x)
    acf = np.correlate(x, x, "full")[n_ - 1:] / (x @ x)
    s = 0
    for k in range(1, n_):
        if acf[k] < 0: break
        s += acf[k]
    return n_ / (1 + 2 * s)
print("  有效样本量:", [int(ess(D[:, k])) for k in range(4)], "/", len(D))

mle = sm.tsa.SARIMAX(c3, exog=c1, order=(2, 0, 0), trend="n").fit(disp=False)
print("MLE (SARIMAX):", np.round(mle.params, 4), "标准误:", np.round(mle.bse, 4))

# 第二条链从很差的初值出发，检查收敛
D2 = gibbs_reg_ar(c3, c1, burn=0, init=(-2.0, np.array([0.9, -0.5]), 1.0))
print("极端初值链前 3 次抽样的 beta:", np.round(D2[:3, 0], 4))
print("两条链后验均值(OLS 初值 / 极端初值):", np.round(D.mean(0), 4), np.round(D2[100:].mean(0), 4))

# ============ 例二：AR(3) 加性异常值的 Gibbs 检测（仿例 12.2，McCulloch & Tsay 1994a）============
rng = np.random.default_rng(3)
n, p = 600, 3
phi_true = np.array([0.227, 0.006, 0.114])
x = np.zeros(n)
for t in range(p, n):
    x[t] = phi_true @ x[t - p:t][::-1] + rng.normal(0, np.sqrt(0.0128))
y = x.copy()
y[200] += 0.45; y[322] -= 0.50                         # 两个加性异常值 (t=201, 323，1 起算)

# 先验：phi~N(0,0.25 I)，v*lam/sigma2~chi2_5 (v*lam=5*0.00256)，eps~Beta(5,95)，beta_t~N(0,0.1)
Sig0inv = np.eye(p + 1) / 0.25
v, vlam, g1, g2, xi2 = 5, 5 * 0.00256, 5, 95, 0.1
def lagmat(xx):
    return np.column_stack([np.ones(n - p)] + [xx[p - i:n - i] for i in range(1, p + 1)])

delta = np.zeros(n, int); bta = np.zeros(n)
phi = np.array([0.0, 0.2, 0.02, 0.1]); sig2, eps = 0.012, 0.05
n_iter, burn = 1050, 50
P_ao = np.zeros(n); B_ao = np.zeros(n); keep = []
for it in range(n_iter):
    xx = y - delta * bta
    # phi | .
    X = lagmat(xx); xt = xx[p:]
    Sinv = X.T @ X / sig2 + Sig0inv; S = np.linalg.inv(Sinv)
    phi = rng.multivariate_normal(S @ (X.T @ xt / sig2), S)
    # sigma2 | .
    res = xt - X @ phi
    sig2 = (vlam + res @ res) / rng.chisquare(v + n - p)
    # eps | delta
    eps = rng.beta(g1 + delta.sum(), g2 + n - delta.sum())
    # delta_h, beta_h 逐点 (12.19)
    psi = np.r_[-1.0, phi[1:]]                          # psi_0=-1, psi_i=phi_i
    for h in range(p, n):
        xs = y - delta * bta; xs[h] = y[h]              # x*_j：第 h 点取观测值本身
        m = min(n - 1, h + p)
        js = np.arange(h, m + 1)
        w = np.array([xs[j] - phi[0] - phi[1:] @ xs[j - p:j][::-1] for j in js])
        ps = psi[js - h]
        l1 = -np.sum((w + ps * bta[h]) ** 2) / (2 * sig2)
        l0 = -np.sum(w ** 2) / (2 * sig2)
        pr1 = 1 / (1 + (1 - eps) / eps * np.exp(l0 - l1))
        delta[h] = rng.random() < pr1
        if delta[h]:
            den = sig2 + np.sum(ps ** 2) * xi2
            bta[h] = rng.normal(-np.sum(ps * w) * xi2 / den, np.sqrt(sig2 * xi2 / den))
        else:
            bta[h] = rng.normal(0, np.sqrt(xi2))
    if it >= burn:
        P_ao += delta; B_ao += delta * bta
        keep.append(np.r_[phi[1:], sig2, eps])
P_ao /= (n_iter - burn)
B_ao = np.where(P_ao > 0, B_ao / np.maximum(P_ao * (n_iter - burn), 1), 0)
keep = np.array(keep)
print("\nAO 模型后验均值 phi1..3, sigma2, eps:", np.round(keep.mean(0), 4))
top = np.argsort(P_ao)[::-1][:4]
print("  真实 AO 位置 t=201, 323 的后验概率:", np.round(P_ao[[200, 322]], 2), " 其余点平均:", round(np.delete(P_ao, [200, 322]).mean(), 3))
for h in top:
    print(f"  t={h + 1:3d}  AO 后验概率={P_ao[h]:.2f}  后验幅度(条件于是 AO)={B_ao[h]:+.3f}")
ols = sm.OLS(y[p:], lagmat(y)).fit()
print("忽略异常值的 OLS phi1..3:", np.round(ols.params[1:], 4), " sigma2:", round(ols.scale, 4))

关键输出：

Gibbs 后验均值(后验标准差):
  beta      0.7881 (0.0141)
  phi1      0.1803 (0.0206)
  phi2     -0.0186 (0.0205)
  sigma2    0.0048 (0.0001)
  有效样本量: [1799, 2000, 2000, 2000] / 2000
MLE (SARIMAX): [ 0.7875  0.1805 -0.0188  0.0046] 标准误: [0.0131 0.0204 0.0199 0.0001]
极端初值链前 3 次抽样的 beta: [0.906  0.7975 0.7711]
两条链后验均值(OLS 初值 / 极端初值): [ 0.7881  0.1803 -0.0186  0.0048] [ 0.7879  0.1799 -0.0176  0.0048]

AO 模型后验均值 phi1..3, sigma2, eps: [0.2109 0.0177 0.1061 0.012  0.0294]
  真实 AO 位置 t=201, 323 的后验概率: [1.   0.98]  其余点平均: 0.023
  t=201  AO 后验概率=1.00  后验幅度(条件于是 AO)=+0.478
  t=323  AO 后验概率=0.97  后验幅度(条件于是 AO)=-0.401
  t=341  AO 后验概率=0.64  后验幅度(条件于是 AO)=-0.337
  t=567  AO 后验概率=0.48  后验幅度(条件于是 AO)=-0.322
忽略异常值的 OLS phi1..3: [2.060e-01 1.000e-04 1.048e-01]  sigma2: 0.0136

读输出：

• 例一：Gibbs 后验均值和后验标准差与 MLE 的点估计和标准误几乎一样，与原书"样本大、模型简单时贝叶斯与经典结果相近"的结论一致。有效样本量接近名义样本量（2000），说明这条链的自相关很弱、混合很好——这得益于条件后验都是闭式的，且 \(\beta\) 与 \(\phi\) 的相关性不强。从 \(\beta=-2\)、\(\sigma^2=1\) 这样离谱的初值出发，第 2 次迭代就回到了 0.8 附近，两条链的后验均值一致。
• 例二：两个真实 AO 的后验概率分别为 1.00 和 0.97，幅度估计 +0.48、−0.40 接近真值 +0.45、−0.50。另有 \(t=341\)、\(t=567\) 两点后验概率在 0.5 上下，查看模拟数据可知它们是无异常序列本身的极端值（约 −3.2、−3.6 倍标准差）：在正态假设下，这样的点确实"像"异常值。\(\epsilon\) 的后验均值约 3%，低于先验均值 5%。忽略异常值的 OLS 把 \(\sigma^2\) 估成 0.0136，比 AO 模型的 0.0120 大，这就是两个异常值把噪声方差"撑大"的效果。
• 实务提醒：真实收益率有厚尾，正态 AO 模型会把一部分正常的尾部观测也标为 AO。用在数据清洗时，应把后验概率作为排查优先级而不是自动删除的依据；也可以把 \(a_t\) 换成 Student-t 分布（用正态尺度混合的数据增广实现）。另外，如果某个 AO 恰好把观测推向其条件均值附近，它既检测不出来，也几乎不影响估计——这种情况换个随机种子就能观察到。

本章小结

MCMC 构造一个平稳分布等于后验的马尔可夫链，用链上的样本做推断。Gibbs 抽样依次从各参数的完全条件分布中抽样，每一步都保持后验不变；实践中要丢弃 burn-in，对高相关参数分块联合抽样，并用多初值、迹图、有效样本量检查收敛。共轭先验让条件后验有闭式，核心规律是"后验精度 = 先验精度 + 数据精度"；没有闭式时用 Metropolis（对称建议）、Metropolis–Hastings（一般建议，接受比含建议密度修正）或 Griddy Gibbs（一元格点逆 CDF），并可嵌入 Gibbs 之中。数据增广把缺失值、异常指示和异常幅度都当作参数：带 AR 误差的回归在准差分回归与 AR 回归之间交替抽样；AR 模型中的缺失值等价于一个 \(p+1\) 个点的无截距回归；AO 检测按先验概率加权比较"是/不是异常"两种情形下的局部似然。

练习

基础

1. （原书习题 12.1）\(x\sim N(\mu,4)\)，先验 \(\mu\sim N(0,25)\)，观测到一个 \(x\)，求 \(\mu\) 的后验。 提示：由 Result 12.1，后验均值 \(25x/29\)，方差 \(100/29\)。
2. 某策略过去 40 个交易日有 26 天盈利。取先验 Beta(10,10)，求胜率的后验分布、后验均值和 95% 可信区间；与均匀先验 Beta(1,1) 比较。 提示：后验 Beta(36,24)，均值 0.6。
3. 证明 Metropolis–Hastings 的转移核满足细致平衡。
4. 对 AR(1) 中的单个缺失值，用 (12.16) 推导滤波值 \(\hat x_h=\frac{\phi_1}{1+\phi_1^2}(x_{h-1}+x_{h+1})\)，并解释为什么前后权重相同。
5. 在例 12.2 的先验设定下，解释 \(\gamma_1=5,\gamma_2=95\) 与 \(\xi^2=0.1\) 的含义；如果把 \(\xi^2\) 改成 0.01，你预期检测结果怎样变化？ 提示：\(\xi^2\) 太小时，先验认为异常幅度很小，大幅异常会被低估、甚至部分归入噪声。

进阶

6. （原书习题 12.2）回归误差为 AR(\(p\)) 时，在共轭先验 \(\beta\sim N(\beta_o,\Sigma_o)\)，\(\phi\sim N(\phi_o,A_o)\)，\(v\lambda/\sigma^2\sim\chi^2_v\) 下推导三个条件后验。
7. （原书习题 12.3）AR(\(p\)) 中 \(x_h,x_{h+1}\) 两个相邻缺失值，给定二元正态先验，推导其联合条件后验。
8. 修改例一代码，把 \(\beta\) 与 \(\phi\) 改为用随机游走 Metropolis 抽样（建议分布标准差分别取 0.01、0.02），比较接受率与有效样本量，体会闭式 Gibbs 步骤的效率优势。
9. 把例二中的 \(a_t\) 改为自由度 4 的 Student-t 噪声（不加任何 AO），运行 AO 检测，统计后验概率 > 0.5 的点有多少个，讨论正态 AO 模型用于厚尾金融数据的局限。
10. （原书习题 12.8）30 年期抵押贷款利率与 3 月期国库券利率（1971-04 至 2009-09）：先用传统方法建立带时间序列误差的回归，再用 MCMC 重新估计并比较。没有数据时，可用例一的模拟框架设定一个 ARMA 误差版本。

原书推荐习题：12.2（AR(\(p\)) 误差回归的条件后验推导，掌握 Gibbs 推导套路）、12.3（连续缺失值的联合后验）、12.1（共轭先验基础计算）、12.8。

原书对照

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