第 24a 章蒙特卡洛积分与重要性抽样

本章与第 24b 章共同对应 Wasserman 原书第 24 章《Simulation Methods》。本章讲前半部分：基本蒙特卡洛积分和重要性抽样；第 24b 章讲马尔可夫链蒙特卡洛（MCMC）。原书以贝叶斯推断中的积分为主线，但这些技术是通用的。对量化交易而言，蒙特卡洛是衍生品定价（路径依赖期权、多资产期权）、风险计量（VaR、ES、信用组合损失）和情景分析的主力工具；重要性抽样则是稀有事件（深度虚值期权、尾部损失、违约相关）估计的标准方法。本章是本册的核心章节之一。

学习目标

理解蒙特卡洛积分的原理（大数定律）和误差（中心极限定理），会计算标准误、构造置信区间，并据此确定所需模拟次数。
理解蒙特卡洛误差 $O(N^{-1/2})$ 与维数无关的含义，知道它何时优于确定性数值积分。
会用后验抽样计算后验均值、后验区间及参数复杂函数的后验分布。
掌握重要性抽样的推导、方差公式和"$g$ 要与 $|h|f$ 相似且尾部更厚"的原则，会证明最优重要性密度。
会用自归一化重要性抽样处理只知道比例常数的目标分布，了解接受–拒绝抽样。
能用蒙特卡洛为期权定价，并用均值平移的重要性抽样估计深度虚值期权和尾部风险。

读前导读

这一章在解决什么问题

结论：本章讲"算不出来的积分，就用随机抽样去平均"，以及当关心的是罕见事件时，怎样聪明地抽样。

你在 CFA 里见过蒙特卡洛模拟：给定收益分布，生成上万条情景，统计组合损失的分布，读出 VaR。本章说明这套做法背后的数学：任何期望（期权价格、VaR 对应的概率、ES、贝叶斯后验均值）都是一个积分，而积分可以写成"某个随机变量的平均值"。按分布抽样、求样本平均，大数定律保证它收敛到真值，中心极限定理告诉你误差有多大、需要多少条路径。

第二个主题是重要性抽样。估计深度虚值期权或 1‰ 的尾部损失时，绝大多数模拟路径落在"什么都没发生"的区域，白白浪费。重要性抽样故意多从尾部抽样，再用权重把结果纠正回来。这个"换一个分布抽样、再用比值加权"的操作，和你熟悉的风险中性定价是同一个数学结构：真实概率下的期望，可以在另一个概率测度下计算，只要乘上两个测度之间的密度比。

需要先想起来的数学

1. 期望就是积分。 连续型随机变量 $\mathbb E_fh(X)=\int h(x)f(x)dx$：把每个可能的 $x$ 处的值 $h(x)$ 按密度 $f(x)$ 加权求和。例：$X\sim\text{Uniform}(0,1)$，$\mathbb EX^3=\int_0^1x^3dx=\frac14$。参见第 00 册第 03 章积分（密度与期望作为积分）。

2. 大数定律与中心极限定理。 样本均值收敛到期望；误差近似正态，标准差为 $\sigma/\sqrt N$。例：单条路径收益标准差 15，1 万条路径的标准误是 $15/100=0.15$。要把它降到 0.015，需要 100 万条路径。

3. Jensen 不等式与 $\mathbb E W^2\ge(\mathbb E|W|)^2$。 后者是方差非负 $\mathbb V|W|=\mathbb EW^2-(\mathbb E|W|)^2\ge0$ 的直接改写。参见第 00 册第 07 章概率中的分析工具的常用不等式。

4. 反常积分的收敛与发散。 积分区间延伸到无穷时，被积函数在尾部衰减得不够快，积分就是无穷。例：$\int_1^\infty x^{-2}dx=1$ 有限，$\int_1^\infty x^{-1}dx=\infty$；$\int e^{-x^2/2}dx$ 有限，$\int e^{+x^2/2}dx$ 无穷。重要性抽样的"方差无穷"就是这种情况。参见第 00 册第 03 章积分的反常积分部分。

怎么读这一章

核心必读：24a.2.1（原理、标准误、维数无关）、24a.3.1–24a.3.2（重要性抽样与方差无穷的陷阱）、24a.4 实战三部分及解读。24a.2.2 的后验抽样对理解第 24b 章很重要，例 24.3 值得细读，例 24.4 可以只看思路和勘误。定理 24.5 的证明、24a.3.3 自归一化、24a.3.4 接受–拒绝第一次可以只看结论。建议顺序：24a.2.1 → 实战 (1) → 24a.3.1–24a.3.2 → 实战 (2)(3) → 24a.2.2 → 其余各节。

24a.1 动机：贝叶斯推断中的积分

给定先验 $f(\theta)$ 与数据 $x^n=(x_1,\dots,x_n)$，后验密度为

f(\theta\mid x^n)=\frac{\mathcal L(\theta)f(\theta)}c,\qquad c=\int\mathcal L(\theta)f(\theta)d\theta,

后验均值为

\bar\theta=\int\theta f(\theta\mid x^n)d\theta=\frac{\int\theta\mathcal L(\theta)f(\theta)d\theta}{\int\mathcal L(\theta)f(\theta)d\theta}.

$\theta=(\theta_1,\dots,\theta_k)$ 是多维时，某个分量的边际后验

f(\theta_1\mid x^n)=\int\cdots\int f(\theta_1,\dots,\theta_k\mid x^n)d\theta_2\cdots d\theta_k

是高维积分，通常没有解析解。贝叶斯推断的计算困难几乎全部集中在积分上，模拟方法就是为此而生的。金融中同样的问题无处不在：期权价格是风险中性测度下的期望（积分），VaR 是损失分布的分位数，ES 是尾部条件期望。

24a.2 基本蒙特卡洛积分

24a.2.1 原理

要计算 $I=\int_a^bh(x)dx$。把它改写为期望：

I=\int_a^bh(x)dx=\int_a^b\underbrace{h(x)(b-a)}_{w(x)}\cdot\underbrace{\frac1{b-a}}_{f(x)}dx=\mathbb E_f\big(w(X)\big),\quad X\sim\text{Uniform}(a,b).\tag{24.1}

抽 $X_1,\dots,X_N\sim\text{Uniform}(a,b)$，由大数定律

\hat I=\frac1N\sum_{i=1}^Nw(X_i)\xrightarrow{P}\mathbb E(w(X))=I.\tag{24.2}

由中心极限定理，标准误和 $1-\alpha$ 置信区间为

\widehat{\text{se}}=\frac s{\sqrt N},\quad s^2=\frac1{N-1}\sum_i\big(Y_i-\hat I\big)^2,\ Y_i=w(X_i);\qquad\hat I\pm z_{\alpha/2}\widehat{\text{se}}.

与统计推断不同的是，这里的 $N$ 由我们决定，可以任意大，所以区间可以任意窄。

推导拆解：(24.1) 的技巧是"乘一个再除一个"：$h(x)=h(x)(b-a)\cdot\frac{1}{b-a}$，后一个因子正好是 $[a,b]$ 上均匀分布的密度，于是积分变成了"$w(X)$ 在均匀分布下的期望"。

用例 24.1 验证标准误：$Y=U^3$，$\mathbb EY=\frac14$，$\mathbb EY^2=\mathbb EU^6=\int_0^1u^6du=\frac17$，$\mathbb VY=\frac17-\frac1{16}\approx0.0804$，标准差 $\approx0.284$。$N=10000$ 时标准误 $0.284/100=0.0028$，与原文一致。

所需样本量的公式 $N\approx(z_{\alpha/2}s/\text{目标半宽})^2$，就是把"半宽 $=z_{\alpha/2}s/\sqrt N$"反解出 $N$。

更一般的形式。 若 $f$ 是任何能直接抽样的概率密度，

I=\int h(x)f(x)dx=\mathbb E_f\,h(X)\approx\frac1N\sum_{i=1}^Nh(X_i),\quad X_i\sim f.\tag{24.3}

误差与维数无关。 蒙特卡洛误差 $\sigma_h/\sqrt N$ 中的 $\sigma_h^2=\mathbb V_f\,h(X)$ 与积分的维数无关。确定性数值积分（梯形法、Simpson 法、Gauss 求积）在一维时精度远高于蒙特卡洛，但在 $d$ 维中需要 $m^d$ 个格点，误差随维数迅速恶化。一个 252 个交易日的路径依赖期权就是一个 252 维积分，此时蒙特卡洛几乎是唯一可行的方法。代价是收敛慢：精度提高 10 倍，计算量要增加 100 倍。

金融直觉：格点法的维数困境可以直接算出来。对 252 维积分，每一维只取 10 个格点，就需要 $10^{252}$ 个点，远超宇宙中原子的数量（约 $10^{80}$）。蒙特卡洛的误差 $\sigma_h/\sqrt N$ 中没有 $d$，10 万条路径在 1 维和 252 维时误差的阶相同（$\sigma_h$ 本身可能随问题变化，但不会随维数爆炸）。这就是为什么美式期权之外的路径依赖产品几乎都用蒙特卡洛定价。

例 24.1。 $h(x)=x^3$，$I=\int_0^1x^3dx=\frac14$。$N=10{,}000$ 个均匀随机数得 $\hat I=0.248$，标准误 0.0028。

例 24.2（正态分布函数）。 $\Phi(x)=\int I(s<x)\phi(s)ds$，于是 $\hat\Phi(x)=\#\{X_i\le x\}/N$，$X_i\sim N(0,1)$。$x=2$ 时真值 0.9772，$N=10{,}000$ 得 0.9751，$N=100{,}000$ 得 0.9771。

24a.2.2 后验抽样：直接得到复杂函数的后验分布

例 24.3（两个比例之差）。 $X\sim\text{Binomial}(n,p_1)$，$Y\sim\text{Binomial}(m,p_2)$，关心 $\delta=p_2-p_1$。频率派方法：MLE $\hat\delta=Y/m-X/n$，Delta 方法给出 $\widehat{\text{se}}=\sqrt{\hat p_1(1-\hat p_1)/n+\hat p_2(1-\hat p_2)/m}$。

贝叶斯方法取平坦先验 $f(p_1,p_2)=1$，后验

f(p_1,p_2\mid X,Y)\propto p_1^X(1-p_1)^{n-X}\,p_2^Y(1-p_2)^{m-Y}.

直接对 $\delta$ 的后验分布积分很繁琐：要先求 $\{p_2-p_1\le c\}$ 上的二重积分再求导。模拟则非常简单：后验可分解为 $f(p_1\mid X)f(p_2\mid Y)$，两者独立，

p_1\mid X\sim\text{Beta}(X+1,n-X+1),\qquad p_2\mid Y\sim\text{Beta}(Y+1,m-Y+1).

抽 $N$ 对 $(p_1^{(i)},p_2^{(i)})$，令 $\delta^{(i)}=p_2^{(i)}-p_1^{(i)}$。后验均值约为 $\frac1N\sum\delta^{(i)}$，后验 95% 区间是 $\delta^{(i)}$ 的 2.5% 与 97.5% 分位数，后验密度可以用 $\delta^{(i)}$ 的直方图或核密度估计（第 20 章）得到。原书数据：$n=m=10$，$X=8$，$Y=6$，95% 后验区间为 $(-0.52,0.20)$。

这是模拟方法最强大的地方：一旦有了参数的后验样本，任何函数 $g(\theta)$ 的后验分布只需把 $g$ 作用到每个样本上。 在量化中，$\theta$ 可以是收益均值和协方差，$g(\theta)$ 可以是最优组合权重、Sharpe 比、VaR——后验样本直接给出这些量的不确定性。

金融直觉：比较两种做法。Delta 方法要先对 $g$ 求导、再做线性近似，$g$ 一复杂（例如"最优权重 $\Sigma^{-1}\mu$ 中第 3 只股票的权重"或"VaR 的分位数"）就很难写出；后验抽样只需"对每一组参数样本算一次 $g$"。例如抽 1 万组 $(\mu,\Sigma)$，每组算一次最优组合，就得到 1 万个权重向量，直接看它们的分布——你会发现均值–方差最优权重的不确定性大得惊人，这正是 Black–Litterman 和收缩估计出现的原因。

例 24.4（剂量–反应与 LD50）。 10 个递增剂量 $x_1<\cdots<x_{10}$，每个剂量 15 只小鼠，$Y_i\sim\text{Binomial}(15,p_i)$ 独立，$p_i$ 是剂量 $x_i$ 下的死亡概率。生物学知识要求 $p_1\le p_2\le\cdots\le p_{10}$。目标是半数致死剂量 LD50：$\delta=x_j$，$j=\min\{i:p_i\ge0.5\}$。它是 $(p_1,\dots,p_{10})$ 的复杂函数，后验均值要在受约束区域 $A=\{p_1\le\cdots\le p_{10}\}$ 上做 10 维积分。

取截断在 $A$ 上的平坦先验，模拟步骤为：(1) 抽 $P_i\sim\text{Beta}(Y_i+1,15-Y_i+1)$，$i=1,\dots,10$；(2) 若 $P_1\le\cdots\le P_{10}$ 则保留，否则丢弃重抽（这是拒绝抽样：从无约束后验抽样、只保留落在约束区域内的样本，等价于从截断后验抽样）；(3) 令 $j=\min\{i:P_i>0.5\}$，$\delta=x_j$。重复 $N$ 次即得 $\delta$ 的后验分布 $\mathbb P(\delta=x_j\mid Y)\approx\frac1N\sum_iI(\delta^{(i)}=x_j)$。

数据为 $Y=(0,0,2,2,8,10,12,14,15,14)$。关于原书的一处问题：原书把这组计数称为"存活数"，但计数随剂量递增，而模型要求 $p_i$（死亡概率）随剂量递增且 $Y_i\sim\text{Binomial}(15,p_i)$，二者只有把 $Y_i$ 理解为死亡数才自洽，本教材按死亡数处理。原书报告的结果为 $\bar\delta=4.04$，95% 区间 $(3,5)$。按上述数据与算法，下面的复现得到的死亡率首次超过 0.5 的剂量序号后验均值约为 5.6，95% 区间为第 5 到第 7 个剂量；已对照原书（PDF p.410）确认剂量就是 1–10、数据与此相同，原书的 4.04 无法复现，可能是原书计算或排印有误（见章末勘误说明）。方法本身不受影响。

代码：复现例 24.1、24.3、24.4、24.6。

import numpy as np
from scipy import stats
rng = np.random.default_rng(24)

# 例 24.1：∫_0^1 x^3 dx = 1/4
u = rng.random(10000); w = u ** 3
print("例24.1  I_hat=%.4f  se=%.4f" % (w.mean(), w.std(ddof=1) / np.sqrt(w.size)))

# 例 24.3：两个二项比例之差的后验（平坦先验 → 两个独立 Beta 后验）
n = m = 10; X, Y = 8, 6
p1 = rng.beta(X + 1, n - X + 1, 100000); p2 = rng.beta(Y + 1, m - Y + 1, 100000)
d = p2 - p1
print("例24.3  后验均值=%.3f  95%%后验区间=(%.2f, %.2f)  P(δ>0|data)=%.3f"
      % (d.mean(), *np.quantile(d, [0.025, 0.975]), (d > 0).mean()))

# 例 24.4：LD50，单调约束下的拒绝抽样（计数按“死亡数”理解）
Yd = np.array([0, 0, 2, 2, 8, 10, 12, 14, 15, 14]); k = 15
keep, tried = [], 0
while sum(len(a) for a in keep) < 5000:
    P = rng.beta(Yd + 1, k - Yd + 1, size=(200000, 10)); tried += 200000
    keep.append(P[np.all(np.diff(P, axis=1) >= 0, axis=1)])
P = np.vstack(keep)
j = np.argmax(P > 0.5, axis=1) + 1                      # j = min{i : P_i > 0.5}
vals, cnts = np.unique(j, return_counts=True)
print("例24.4  接受率=%.4f  E[j|data]=%.2f  95%%区间=(%d, %d)" % (len(P) / tried, j.mean(), *np.quantile(j, [0.025, 0.975])))
print("        后验 pmf:", dict(zip(vals.tolist(), (cnts / len(j)).round(3).tolist())))

# 例 24.6：P(Z>3)，基本 MC 与重要性抽样（g = N(4,1)），N=100，重复 20000 次
N, R = 100, 20000
I1 = (rng.standard_normal((R, N)) > 3).mean(1)
x = rng.normal(4, 1, (R, N)); I2 = (stats.norm.pdf(x) / stats.norm.pdf(x, 4, 1) * (x > 3)).mean(1)
print("例24.6  真值=%.5f  基本MC: 均值=%.5f 标准差=%.5f   重要性抽样: 均值=%.5f 标准差=%.5f"
      % (stats.norm.sf(3), I1.mean(), I1.std(), I2.mean(), I2.std()))

关键输出：

例24.1  I_hat=0.2536  se=0.0028
例24.3  后验均值=-0.166  95%后验区间=(-0.51, 0.20)  P(δ>0|data)=0.182
例24.4  接受率=0.0072  E[j|data]=5.58  95%区间=(5, 7)
        后验 pmf: {4: 0.002, 5: 0.468, 6: 0.479, 7: 0.051, 8: 0.0}
例24.6  真值=0.00135  基本MC: 均值=0.00142 标准差=0.00374   重要性抽样: 均值=0.00135 标准差=0.00031

例 24.1 的标准误 0.0028 与原书一致；例 24.3 的后验区间 $(-0.51,0.20)$ 与原书的 $(-0.52,0.20)$ 在模拟误差内一致，并且顺手得到 $\mathbb P(p_2>p_1\mid\text{data})=0.18$。例 24.4 的接受率只有 0.7%：约束区域在 10 维空间中只占后验质量的很小一部分，拒绝抽样虽然正确但浪费，这正是第 24b 章 Gibbs 抽样要解决的问题。

24a.3 重要性抽样

24a.3.1 原理

很多时候我们不会从 $f$ 抽样（例如后验 $\propto\mathcal L(\theta)f(\theta)$ 往往不是标准分布），或者从 $f$ 抽样效率很低（例如要估计的事件极少发生）。重要性抽样（importance sampling）改从一个我们会抽样的密度 $g$ 抽样，再用权重修正：

I=\int h(x)f(x)dx=\int\frac{h(x)f(x)}{g(x)}g(x)dx=\mathbb E_g(Y),\quad Y=\frac{h(X)f(X)}{g(X)}.\tag{24.4}

抽 $X_1,\dots,X_N\sim g$， $\hat I=\frac1N\sum_{i=1}^N\frac{h(X_i)f(X_i)}{g(X_i)}.\tag{24.5}$ 大数定律保证 $\hat I\xrightarrow PI$（只要在 $hf\ne0$ 的地方 $g>0$）。$f(X_i)/g(X_i)$ 称为重要性权重：$g$ 过多抽到的区域权重小于 1，抽少了的区域权重大于 1。

金融直觉：这和风险中性定价是同一个操作。真实测度 $\mathbb P$ 下的期望可以改在风险中性测度 $\mathbb Q$ 下计算，只要乘上两者的密度比（Radon–Nikodym 导数，也就是定价核）：$\mathbb E_{\mathbb P}[X]=\mathbb E_{\mathbb Q}\big[X\cdot\frac{d\mathbb P}{d\mathbb Q}\big]$。重要性抽样中，$f$ 相当于 $\mathbb P$，$g$ 相当于 $\mathbb Q$，权重 $f/g$ 就是密度比。区别在于定价时换测度是为了让贴现价格成为鞅，而重要性抽样换测度是为了让模拟更有效率。

推导拆解：实战中的均值平移 $g=N(\mu,1)$、$f=N(0,1)$，权重有简洁的形式：$\frac{f(x)}{g(x)}=\frac{e^{-x^2/2}}{e^{-(x-\mu)^2/2}}=e^{-\mu x+\mu^2/2}$（展开 $(x-\mu)^2=x^2-2\mu x+\mu^2$，$x^2$ 项相消）。这种指数形式的权重叫"指数倾斜"，在连续时间中对应 Girsanov 定理的漂移变换。$\mu$ 取负值（平移到左尾）时，左尾样本的权重 $e^{-\mu x+\mu^2/2}$ 远小于 1，正好抵消被过量抽到的部分。

24a.3.2 陷阱：方差可能无穷

$\hat I$ 的方差为 $\frac1N\mathbb V_g(Y)$，其中

\mathbb E_g(Y^2)=\int\Big(\frac{h(x)f(x)}{g(x)}\Big)^2g(x)dx=\int\frac{h^2(x)f^2(x)}{g(x)}dx.\tag{24.6}

若 $g$ 的尾部比 $f$ 薄，在尾部 $f^2/g$ 可能比 $f$ 大得多，积分可能发散，$\hat I$ 的标准误为无穷。这时 $\hat I$ 仍然相合，但收敛极其缓慢且不稳定：大多数时候估计偏低（尾部没抽到），偶尔抽到一个尾部点，权重巨大，估计突然跳高。更隐蔽的是，样本标准误在有限样本中总是有限的，看不出问题。

推导拆解：一个可以手算的例子。$f=N(0,1)$，$g=N(0,\sigma^2)$，$h\equiv1$（只是检验权重本身）。$\frac{f^2}{g}\propto\exp\big(-x^2+\frac{x^2}{2\sigma^2}\big)=\exp\big(x^2(\frac{1}{2\sigma^2}-1)\big)$。若 $\sigma^2<\frac12$，指数中 $x^2$ 的系数为正，被积函数在尾部像 $e^{+cx^2}$ 一样爆炸，(24.6) 的积分为无穷。也就是说，提议分布的标准差只要比目标小 30% 以上（$\sigma<0.707$），重要性权重的方差就是无穷。反过来，$\sigma>1$（$g$ 更胖）时系数为负，积分总是有限的。

为什么"样本标准误总是有限"：样本方差只用了你抽到的点，而导致无穷方差的是那些概率极小、权重极大、你多半还没抽到的点。因此不能靠样本标准误发现问题，要从 $f$、$g$ 的尾部形状事先判断。

原则：从尾部比 $f$ 更厚的 $g$ 抽样。 另一方面，若 $g$ 在 $f$ 大的地方很小，$f/g$ 也会很大。所以：好的 $g$ 形状与 $f$（更准确地说与 $|h|f$）相似，但尾部更厚。

定理 24.5（最优重要性密度）。 使 $\hat I$ 方差最小的 $g$ 为

g^*(x)=\frac{|h(x)|f(x)}{\int|h(s)|f(s)ds}.

证明。 $W=fh/g$ 的方差为 $\int\frac{h^2f^2}g-\big(\int hf\big)^2$，第二项与 $g$ 无关。由 Jensen 不等式，$\mathbb E_g(W^2)\ge\big(\mathbb E_g|W|\big)^2=\big(\int|h|f\big)^2$，这是任何 $g$ 都无法突破的下界。代入 $g^*$ 直接计算：$\mathbb E_{g^*}(W^2)=\int\frac{h^2f^2}{|h|f}\int|h|f=\big(\int|h|f\big)^2$，恰好达到下界。$\square$

推导拆解：证明中三步的依据。(1) 方差 = 二阶矩 − 均值²，而 $\mathbb E_gW=\int\frac{hf}{g}g=\int hf=I$ 与 $g$ 无关，所以只需最小化二阶矩。(2) 下界：对任何随机变量，$\mathbb E W^2\ge(\mathbb E|W|)^2$（即 $|W|$ 的方差非负），且 $\mathbb E_g|W|=\int\frac{|h|f}{g}g=\int|h|f$，同样与 $g$ 无关。(3) 代入 $g^*=|h|f/c$（$c=\int|h|f$）：$W=\frac{hf}{|h|f/c}=c\cdot\text{sign}(h)$，$W^2=c^2$ 恒定，二阶矩恰为 $c^2$，达到下界。

白话解释：$g^*$ 的意思是"按贡献大小分配抽样预算"：哪里 $|h|f$ 大（对积分贡献大），就在哪里多抽。深虚值看跌期权的 $|h|f$ 集中在行权边界 $z^*$ 左侧附近，所以把 $g$ 平移到 $z^*$ 附近，就是在近似 $g^*$。

若 $h\ge0$，$g^*$ 下的方差为零——每个样本都恰好等于 $I$。当然 $g^*$ 的归一化常数正是我们要求的 $I$，所以它只有理论意义。实践中的指导是：找一个易于抽样、形状与 $|h|f$ 相似、尾部更厚的 $g$。

例 24.6（尾概率）。 $I=\mathbb P(Z>3)=0.0013$，$h(x)=I(x>3)$，$f$ 为标准正态。基本蒙特卡洛 $N=100$，重复多次得 $\mathbb E\hat I=0.0015$，波动约 0.0039——绝大部分样本浪费在远离右尾的地方，100 个样本中平均只有 0.13 个落在 3 以右。改用 $g=N(4,1)$：样本集中在关心的区域，$\hat I=N^{-1}\sum f(X_i)h(X_i)/g(X_i)$ 的均值 0.0011，波动 0.0002。关于原书的一处问题：原书把"0.0039"和"0.0002"标为方差 $\mathbb V(\hat I)$，按数值它们实为标准差（基本 MC 的理论标准差为 $\sqrt{0.00135\times0.99865/100}=0.0037$）。上面的复现给出基本 MC 标准差 0.0037、重要性抽样 0.0003，标准差约降低 12 倍（方差约降低 150 倍）；原书报告约 20 倍，差异来自原书单次模拟的随机性。

24a.3.3 自归一化重要性抽样

贝叶斯后验通常只知道到比例常数：$f(\theta\mid x)\propto\mathcal L(\theta)f(\theta)$。这时把分子分母都用重要性抽样估计：

\bar\theta=\frac{\int\theta\mathcal L(\theta)f(\theta)d\theta}{\int\mathcal L(\theta)f(\theta)d\theta}\approx\frac{\sum_j\theta_j\mathcal L(\theta_j)f(\theta_j)/g(\theta_j)}{\sum_j\mathcal L(\theta_j)f(\theta_j)/g(\theta_j)},\quad\theta_j\sim g.

归一化常数在分子分母中抵消。这个比值估计量有 $O(1/N)$ 的偏差，但相合。

白话解释：把 $w_j=\mathcal L(\theta_j)f(\theta_j)/g(\theta_j)$ 归一化成 $\tilde w_j=w_j/\sum_kw_k$，估计量就是 $\sum_j\tilde w_j\theta_j$——一个加权平均，权重和为 1，就像按市值加权的指数。未知常数 $c$ 会让所有 $w_j$ 同乘一个数，归一化后消失。有偏的原因是分母本身也是随机的：比值的期望一般不等于期望的比值；但分子分母都收敛，比值也收敛（相合）。

实务中常用"有效样本量" $\text{ESS}=1/\sum_j\tilde w_j^2$ 检查权重是否过于集中：权重均匀时 ESS $=N$；若一个样本拿走了几乎全部权重，ESS 接近 1，说明 $g$ 选得不好。

例 24.7（含离群值的测量模型）。 测量 $X_i=\theta+\epsilon_i$。正态误差的尾部太薄，不适合偶有"野值"的测量；改用自由度 $\nu=3$ 的 $t$ 分布误差，似然 $\mathcal L(\theta)=\prod_it_3(X_i-\theta)$，平坦先验。后验均值没有闭式解，用自归一化重要性抽样计算。原书在一个很小的数据集上比较：数值积分给出的后验均值为 $-0.54$；用正态分布作 $g$ 得 $-0.74$（尾部太薄，失败）；用 Cauchy（$t_1$）作 $g$ 得 $-0.53$（尾部厚，准确）。这再次印证了"$g$ 的尾部要比目标厚"。

24a.3.4 接受–拒绝抽样

原书习题 3 介绍了另一种基础抽样方法。若存在易抽样的 $g$ 和常数 $M$ 使 $f(x)\le Mg(x)$ 处处成立：

抽 $X\sim g$，$U\sim\text{Uniform}(0,1)$；
若 $U\le\frac{f(X)}{Mg(X)}$，接受 $Y=X$；否则回到第 1 步。

正确性。 $\mathbb P(X\le y,\text{接受})=\int_{-\infty}^yg(x)\frac{f(x)}{Mg(x)}dx=\frac{F(y)}M$，令 $y\to\infty$ 得接受概率 $1/M$，相除得 $\mathbb P(Y\le y)=F(y)$。$\square$

白话解释：几何图像最直观。画出曲线 $Mg(x)$，它处处盖住 $f(x)$。在 $Mg$ 下方的区域里均匀撒点：先按 $g$ 选横坐标 $X$，再在 $[0,Mg(X)]$ 上均匀选纵坐标 $UMg(X)$。落在 $f$ 曲线下方的点保留（即 $U\le f/(Mg)$），其余丢弃。保留下来的点在 $f$ 下方均匀分布，其横坐标的密度正好是 $f$。$Mg$ 下方的面积是 $M$，$f$ 下方的面积是 1，所以接受率为 $1/M$。

每次接受平均需要 $M$ 次尝试，所以 $g$ 越接近 $f$（$M$ 越接近 1）越高效。例 24.4 的约束抽样就是 $M$ 很大（约 140）的情形。

补充：方差缩减。 原书未展开、但在定价实务中与重要性抽样同样常用的还有：对偶变量（同时使用 $Z$ 与 $-Z$，利用负相关抵消误差）、控制变量（用一个期望已知且与目标高度相关的量，如几何平均亚式期权之于算术平均亚式期权，回归掉部分噪声）、分层抽样和准蒙特卡洛（低差异序列，误差接近 $O(N^{-1})$）。详见第 08 册。

24a.4 量化实战：期权定价与尾部风险

场景。 几何布朗运动 $S_T=S_0\exp\{(r-\sigma^2/2)T+\sigma\sqrt TZ\}$，$S_0=K=100$，$r=3\%$，$\sigma=20\%$，$T=1$。(1) 用基本蒙特卡洛为欧式看涨期权定价，与 Black–Scholes 闭式解比较，观察误差随 $N$ 的变化；再为没有闭式解的算术平均亚式期权定价，计算达到指定精度所需的路径数。(2) 深度虚值看跌期权（$K=60$，行权概率约 0.46%）：基本 MC 与均值平移的重要性抽样。(3) 厚尾损失 $L\sim t_3$ 的超额损失 $\mathbb E(L-4)_+$（ES 计算的核心量）：比较尾部比目标薄的正态 $g$ 与厚尾的 $t$ 型 $g$。

import numpy as np
from scipy import stats
rng = np.random.default_rng(241)
S0, K, r, sig, T = 100.0, 100.0, 0.03, 0.2, 1.0

def bs_call(S, K, r, sig, T):
    d1 = (np.log(S / K) + (r + sig ** 2 / 2) * T) / (sig * np.sqrt(T)); d2 = d1 - sig * np.sqrt(T)
    return S * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2)

# ---------- (1) 基本 MC：欧式看涨（有闭式解可对照）与算术平均亚式看涨（无闭式解） ----------
for N in [1_000, 10_000, 100_000]:
    Z = rng.standard_normal(N)
    ST = S0 * np.exp((r - sig ** 2 / 2) * T + sig * np.sqrt(T) * Z)
    pay = np.exp(-r * T) * np.maximum(ST - K, 0)
    print("欧式看涨 N=%6d  MC=%.3f ± %.3f (1.96se)   BS闭式=%.3f" % (N, pay.mean(), 1.96 * pay.std(ddof=1) / np.sqrt(N), bs_call(S0, K, r, sig, T)))
N, steps = 50_000, 252
dt = T / steps
Z = rng.standard_normal((N, steps))
paths = S0 * np.exp(np.cumsum((r - sig ** 2 / 2) * dt + sig * np.sqrt(dt) * Z, axis=1))
pay = np.exp(-r * T) * np.maximum(paths.mean(1) - K, 0)
se = pay.std(ddof=1) / np.sqrt(N)
print("亚式看涨 MC=%.3f  se=%.3f  → 要把 se 降到 0.005 需 N≈%d" % (pay.mean(), se, int(N * (se / 0.005) ** 2)))

# ---------- (2) 尾部事件：深度虚值看跌（K=60）的价格，均值平移的重要性抽样 ----------
Kp = 60.0
z_star = (np.log(Kp / S0) - (r - sig ** 2 / 2) * T) / (sig * np.sqrt(T))   # 行权边界对应的 z
put_bs = bs_call(S0, Kp, r, sig, T) - S0 + Kp * np.exp(-r * T)
def put_pay(z): return np.exp(-r * T) * np.maximum(Kp - S0 * np.exp((r - sig ** 2 / 2) * T + sig * np.sqrt(T) * z), 0)
N = 10_000
z = rng.standard_normal(N); y0 = put_pay(z)
print("深虚值看跌 BS=%.5f  z*=%.2f  行权概率=%.5f" % (put_bs, z_star, stats.norm.cdf(z_star)))
print("  基本MC        : %.5f ± %.5f   (有正收益的样本 %d 个)" % (y0.mean(), 1.96 * y0.std(ddof=1) / np.sqrt(N), (y0 > 0).sum()))
for name, mu, s in [("IS g=N(z*,1)", z_star, 1.0), ("IS g=z*+t3", None, None)]:
    if mu is None:
        x = z_star + rng.standard_t(3, N); g = stats.t.pdf(x - z_star, 3)
    else:
        x = rng.normal(mu, s, N); g = stats.norm.pdf(x, mu, s)
    w = stats.norm.pdf(x) / g; y = put_pay(x) * w
    print("  %-14s: %.5f ± %.5f   方差缩减倍数≈%.0f" % (name, y.mean(), 1.96 * y.std(ddof=1) / np.sqrt(N), y0.var() / y.var()))

# ---------- (3) 尾部比 f 薄的 g 为什么危险：厚尾损失 L~t3 的超额损失 E[(L-k)+] ----------
from scipy import integrate
k = 4.0
true = integrate.quad(lambda v: (v - k) * stats.t.pdf(v, 3), k, np.inf)[0]
est = {"g=N(5,1.5²) 薄尾": [], "g=5+1.5·t3 厚尾": []}
for _ in range(500):
    x = rng.normal(k + 1, 1.5, 2000); g = stats.norm.pdf(x, k + 1, 1.5)
    est["g=N(5,1.5²) 薄尾"].append((np.maximum(x - k, 0) * stats.t.pdf(x, 3) / g).mean())
    x = k + 1 + 1.5 * rng.standard_t(3, 2000); g = stats.t.pdf((x - k - 1) / 1.5, 3) / 1.5
    est["g=5+1.5·t3 厚尾"].append((np.maximum(x - k, 0) * stats.t.pdf(x, 3) / g).mean())
print("E[(L-4)+], L~t3 真值=%.4f" % true)
for name, v in est.items():
    v = np.array(v)
    print("  %-16s 500次估计: 中位数=%.4f 均值=%.4f 标准差=%.4f 最大/中位数=%.1f" % (name, np.median(v), v.mean(), v.std(), v.max() / np.median(v)))