4. 马尔可夫链蒙特卡洛（MCMC）入门：Metropolis-Hastings算法与Gibbs采样，PyMC库实战
说实话，贝叶斯统计在资产定价里最大的拦路虎是什么？
不是模型有多复杂，而是后验分布算不出来。
你想想看，我们搞金融的，手里一堆收益率数据，想估计某个因子的风险暴露。用贝叶斯方法，先验分布好设，似然函数也能写，但一算后验——完了，分母那个积分根本解不出解析解。
这时候MCMC就登场了。它不直接算后验，而是通过采样的方式，从后验分布里抽出一堆样本。样本多了，后验的均值、方差、分位数就都有了。说白了，就是用计算换解析。
核心思想：MCMC = 马尔可夫链（Markov Chain）+ 蒙特卡洛（Monte Carlo）。马尔可夫链负责“怎么走”，蒙特卡洛负责“算平均”。两者结合，就能从复杂分布中高效采样。
4.1 为什么资产定价需要MCMC？
我在做多因子模型的时候遇到过这么个事。当时用Fama-French五因子模型，想估计每个因子的贝塔系数。传统OLS回归当然能做，但有个问题——我想加入一些先验信息，比如某些因子在熊市里应该更显著。
这就得用贝叶斯了。但后验分布是个高维的、非标准的形式。你没法直接积分。怎么办？
MCMC就是干这个的。它能在高维空间里游走，每一步都朝着后验密度高的地方去，最终留下的轨迹就是后验分布的近似样本。
- 高维积分问题：资产定价模型动辄几十个参数，后验积分维度太高
- 非线性关系：很多定价模型是非线性的，比如随机波动率模型
- 先验信息融合：MCMC能自然地把先验分布纳入采样过程
4.2 Metropolis-Hastings算法：最通用的MCMC
Metropolis-Hastings（MH）算法，我习惯叫它“试错采样法”。
它的逻辑特别直观：
- 从当前位置出发，随机走一步（提议一个新位置）
- 看看新位置的后验密度是不是更高
- 如果是，就跳过去；如果不是，以一定概率接受（避免困在局部最优）
这个“一定概率”怎么算？看公式：
接受概率 = min(1, 后验密度(新位置) / 后验密度(当前位置))说白了，新位置比当前位置好，就100%接受。新位置差一些，也有机会接受——但机会大小取决于差多少。
我的经验：提议分布的方差很关键。方差太小，采样效率低，半天走不出一个山头；方差太大，拒绝率太高，原地踏步。我一般先跑一小段，看接受率——理想值在20%-50%之间。
4.3 Gibbs采样：高维问题的利器
Gibbs采样是MH算法的一个特例，但用起来更爽。
它的思路是：高维联合分布不好采，那就一个一个维度来。每次只更新一个参数，其他参数固定住。这样就把高维问题拆成了一堆一维问题。
举个例子。在资产定价中，我们经常要同时估计多个因子的贝塔系数和残差方差。用Gibbs采样：
- 先固定方差，采样贝塔
- 再固定贝塔，采样方差
- 反复迭代，直到收敛
每个步骤都是一维的条件分布，往往有解析形式，采样效率极高。
注意：Gibbs采样要求每个条件分布都能直接采样。如果某个条件分布不是标准形式，那就得用MH了。我遇到过这种情况——条件分布是截断的t分布，没法直接采，只好在Gibbs里嵌套一个MH步骤。
4.4 知识体系与核心逻辑
下面这张图，是我自己梳理的MCMC在资产定价中的应用逻辑。你看一眼就明白了：
4.5 PyMC库实战：三分钟跑通第一个MCMC
理论说完了，咱们直接上手。PyMC是我最常用的贝叶斯建模工具，没有之一。
下面这个例子，模拟一个简单的资产定价场景：用市场收益率解释某只股票的收益率。
import pymc as pm
import numpy as np
import arviz as az
# 生成模拟数据
np.random.seed(42)
n = 200
market_return = np.random.normal(0.01, 0.02, n)
true_alpha = 0.0005
true_beta = 1.2
true_sigma = 0.015
stock_return = true_alpha + true_beta * market_return + np.random.normal(0, true_sigma, n)
# 构建贝叶斯模型
with pm.Model() as capm_model:
# 先验分布
alpha = pm.Normal('alpha', mu=0, sigma=0.01)
beta = pm.Normal('beta', mu=1, sigma=0.5)
sigma = pm.HalfNormal('sigma', sigma=0.02)
# 似然函数
mu = alpha + beta * market_return
likelihood = pm.Normal('likelihood', mu=mu, sigma=sigma, observed=stock_return)
# MCMC采样
trace = pm.sample(1000, tune=500, chains=2, progressbar=False)
# 查看结果
summary = az.summary(trace, hdi_prob=0.95)
print(summary)这段代码干了什么？
- 定义了alpha、beta、sigma三个参数的先验分布
- 用市场收益率和线性模型构建似然函数
- PyMC自动选择采样算法（默认是NUTS，一种改进的MH算法）
- 跑完1000步采样，得到后验样本
避坑指南：我曾经犯过一个低级错误——采样步数设得太少。100步就停了，结果后验分布根本没收敛。后来养成习惯：至少跑1000步，前一半作为“热身期”（tune）丢弃掉。
4.6 诊断：你的MCMC收敛了吗？
MCMC最怕什么？怕你以为收敛了，其实没有。
我常用的诊断指标有三个：
| 指标 | 含义 | 合格标准 |
|---|---|---|
| R-hat | 多条链之间的方差 vs 链内方差 | < 1.01 |
| 有效样本量 | 实际独立的样本数 | > 总样本的10% |
| 迹图 | 采样值随时间的变化轨迹 | 看起来像“白噪声” |
说白了，R-hat接近1说明多条链达成共识了。有效样本量太低，说明采样效率差，得调整提议分布或者跑更长时间。
重要提醒：千万别只看R-hat。我见过R-hat=1.00但迹图明显有趋势的情况——那是链卡在局部区域了。一定要结合迹图一起看。
4.7 我的MCMC实战清单
最后，分享一份我自己的MCMC检查清单。每次跑模型前，我都会过一遍：
- 先验合理吗？ 别设太紧，也别太松。比如beta的先验，我一般用N(1, 0.5)，因为市场贝塔通常在1附近
- 链数够吗？ 至少2条链，4条更好。方便诊断R-hat
- 热身期够长吗？ 热身期至少占总步数的1/3
- 接受率正常吗？ 对于MH算法，20%-50%是理想区间
- 后验预测检查做了吗？ 用后验样本生成预测数据，看和真实数据是否吻合
嗯，这一章的内容就到这。MCMC是个大话题，但核心就两件事：怎么走（算法），怎么判断走好了（诊断）。把这两件事搞明白，贝叶斯资产定价的大门就打开了。