第 21b 章 数值方法之二：蒙特卡洛模拟与方差缩减

树和有限差分从到期日向后倒推，适合美式期权，但一旦收益依赖整条价格路径（亚式期权、回望期权），或者依赖三个以上的标的变量（一篮子期权、多资产结构化产品），它们就力不从心了。**蒙特卡洛模拟（Monte Carlo simulation）**从今天出发向前随机生成大量价格路径，在每条路径上算出收益，平均后贴现。它的计算量大致随标的变量个数线性增长，而树和有限差分是指数增长，因此在高维问题上几乎是唯一可行的方法。

蒙特卡洛的代价是收敛慢：误差与试验次数的平方根成反比，要把精度提高 10 倍，需要 100 倍的计算量。所以本章的第二部分——方差缩减——与第一部分同样重要。

对量化交易者来说，本章的技术用途远不止期权定价。模拟 GBM 路径、用 Cholesky 分解生成相关的随机冲击、用标准误判断模拟是否充分，这些是风险模型（第 22 章的蒙特卡洛 VaR）、策略压力测试和回测稳健性分析的通用工具。

学习目标

1. 写出风险中性定价下的蒙特卡洛估值步骤，会用 \(\ln S\) 的精确离散化 (21.16)(21.17) 生成路径。
2. 会计算蒙特卡洛估计的标准误和置信区间，理解 \(1/\sqrt M\) 收敛速度对计算成本的含义。
3. 会用 Cholesky 分解生成相关的正态随机数，理解相关矩阵非正定时为什么无解。
4. 知道如何在二叉树上抽样，以及如何用公共随机数估计希腊字母。
5. 掌握六种方差缩减技术：对偶变量、控制变量、重要性抽样、分层抽样、矩匹配、准随机序列，并能判断各自适用的场合。

读前导读

这一章在解决什么问题。 期权的价值是"风险中性期望收益的贴现值"，这个结论你在 CFA 的二叉树里已经用过。当期望算不出封闭解时，最朴素的办法就是：模拟很多次，取平均。这就是蒙特卡洛。它的统计基础你也熟悉：样本均值是总体均值的无偏估计，样本均值的标准误是 \(s/\sqrt n\)，95% 置信区间是均值加减 1.96 个标准误。本章把这套 CFA 一级就学过的抽样理论用在期权定价上。

难点不在"取平均"，而在两件事。第一，怎样生成"正确"的随机路径：股价要服从风险中性世界的对数正态分布，多个资产之间还要有正确的相关性。这要用到 Itô 引理的结论和 Cholesky 分解。第二，怎样用更少的模拟次数达到同样精度。由于误差只按 \(1/\sqrt M\) 下降，暴力加次数很贵，于是有了六种方差缩减技术。它们的共同思想是：要么让每次模拟本身波动更小（对偶、控制变量、重要性抽样），要么让样本在分布上铺得更均匀（分层、矩匹配、准随机）。

这些技术在期权定价之外同样常用：蒙特卡洛 VaR、策略压力测试、组合情景分析，都要生成相关的随机冲击并判断模拟是否足够。

需要先想起来的数学。

• 期望、方差与协方差的运算规则。 \(\mathrm{Var}(aX+bY)=a^2\mathrm{Var}(X)+b^2\mathrm{Var}(Y)+2ab\,\mathrm{Cov}(X,Y)\)。这正是 CFA 两资产组合方差的公式。对偶变量和控制变量的原理都只是这个公式的应用。见 第 00 册第 07 章 概率中的分析工具。
• 中心极限定理与标准误。 \(M\) 个独立同分布样本的均值近似服从正态，标准差为 \(\omega/\sqrt M\)。例：单次模拟标准差 8，模拟 10,000 次，标准误 \(8/100=0.08\)。
• 逆变换法。 若 \(U\) 是 \([0,1]\) 均匀随机数，\(N^{-1}(U)\) 就是标准正态随机数（\(N^{-1}\) 是标准正态累积分布函数的反函数，Excel 的 NORMSINV）。直观：\(U=0.975\) 映到 \(1.96\)，\(U=0.5\) 映到 0。分层抽样、重要性抽样和准随机序列都靠它把均匀数变成正态数。
• 矩阵乘法与 Cholesky 分解。 若 \(\mathbf x\) 的分量独立、方差为 1，则 \(L\mathbf x\) 的协方差矩阵是 \(LL^{\mathsf T}\)。Cholesky 分解就是把给定的相关矩阵 \(R\) 写成"下三角矩阵乘以它的转置"。正定矩阵指对任意非零向量 \(\mathbf w\) 都有 \(\mathbf w^{\mathsf T}R\mathbf w>0\)，金融含义是"任何组合的方差都为正"。见 第 00 册第 06 章 线性代数速成。
• 对数与指数。 \(\ln S\) 的增量正态 ⇔ \(S\) 对数正态。\(E[e^{X}]=e^{\mu+\sigma^2/2}\)（\(X\sim N(\mu,\sigma^2)\)），这是 \(-\sigma^2/2\) 修正项的来源。见 第 00 册第 04 章 级数与收敛。

怎么读这一章。 21.6 节的"路径的生成""生成相关的正态随机数""试验次数与标准误"三小节是核心必读，其中 \(-\sigma^2/2\) 的来源和 Cholesky 的构造要真正弄懂。"在树上抽样"和"计算希腊字母"可以快速读过，但"公共随机数"的道理要记住。21.7 节建议先读对偶变量和控制变量（两者原理最清楚、最常用），再读重要性抽样；分层、矩匹配、准随机可以先看方法选择小结表，需要时再回来细读。最后的代码输出是很好的"看得见的证据"，建议对照"读输出"逐条核对。

21.6 蒙特卡洛模拟

从计算 \(\pi\) 说起

原书 Business Snapshot 21.1 用一个简单例子说明蒙特卡洛的思想。在边长为 1 的正方形内画一个半径 0.5 的内切圆，圆面积为 \(\pi/4\)。向正方形随机投掷飞镖，落入圆内的比例乘以 4 就是 \(\pi\) 的估计。

在 Excel 中：A1、B1 = RAND() 生成两个 \([0,1]\) 均匀随机数作为坐标；C1 = IF((A1−0.5)^2+(B1−0.5)^2<0.5^2, 4, 0)；把这一行复制 100 行，C102 = AVERAGE(C1:C100) 得 3.04，C103 = STDEV(C1:C100) 得 1.69（原书表 21.1）。增加试验次数可以提高精度，但收敛到 3.14159 非常慢。

这个例子包含了蒙特卡洛的全部要素：一个可以写成期望值的量（\(\pi=4\times P(\text{落在圆内})\)）、大量独立随机抽样、样本均值作为估计、样本标准差衡量精度。

用于衍生品估值

蒙特卡洛估值基于风险中性定价：在风险中性世界中抽样路径，求期望收益，再按无风险利率贴现。对于只依赖单一市场变量 \(S\)、在 \(T\) 时刻支付、利率恒定的衍生品：

1. 在风险中性世界中抽样 \(S\) 的一条随机路径；
2. 计算衍生品在这条路径上的收益；
3. 重复第 1–2 步，得到大量收益样本；
4. 求样本均值，作为风险中性期望收益的估计；
5. 按无风险利率贴现，得到衍生品价值的估计。

路径的生成

风险中性世界中，\(S\) 服从

\(\hat\mu\) 为风险中性期望收益率：不付息股票 \(\hat\mu=r\)，支付收益率 \(q\) 的资产 \(\hat\mu=r-q\)，汇率 \(\hat\mu=r-r_f\)，等等。风险中性世界与真实世界的波动率相同（第 13 章 13.7 节）。

朴素离散化。把期限分成 \(N\) 段，每段 \(\Delta t\)：

\(\epsilon\) 是标准正态随机数。一次试验需要 \(N\) 个独立的正态样本来构造完整路径。这个离散化只有在 \(\Delta t\to0\) 时才精确。

白话解释：(21.13) 中的 \(dz\) 是维纳过程（布朗运动）的增量，可以理解为"一段极短时间里的标准化随机冲击"，它在长度 \(\Delta t\) 的区间上的变化是 \(\epsilon\sqrt{\Delta t}\)。所以 (21.14) 读作：每一小段时间里，股价的百分比变化 = 漂移 \(\hat\mu\Delta t\) + 随机部分 \(\sigma\epsilon\sqrt{\Delta t}\)。这正是你熟悉的"收益率 = 期望收益 + 波动率 × 标准正态 × \(\sqrt{时间}\)"。 为什么"只有 \(\Delta t\to0\) 时才精确"：(21.14) 让每步的变化服从正态，所以 \(S\) 有一定概率跳成负数，而且多步复合后分布并不是精确的对数正态。步长越小，这种偏差越小。

对数离散化。实务中模拟 \(\ln S\) 更精确。由 Itô 引理

\(\ln S\) 的漂移和波动率都是常数，所以它的增量精确服从正态分布：

当 \(\hat\mu,\sigma\) 为常数时，(21.16) 对任意 \(\Delta t\) 都精确成立。特别地，若收益只依赖终值，可以一步走到 \(T\)：

(21.17) 可以为在 \(T\) 时刻支付任意非标准收益的衍生品估值，也可以用来检验 BSM 公式。

一个常见错误：在 (21.16) 中忘记 \(-\sigma^2/2\) 项。那样模拟出的 \(S(T)\) 的期望会是 \(S_0e^{(\hat\mu+\sigma^2/2)T}\)，而不是 \(S_0e^{\hat\mu T}\)，价格系统性偏高。检验方法很简单：模拟的 \(e^{-rT}S(T)\) 的均值应等于 \(S_0\)（不付息时）。

推导拆解：\(-\sigma^2/2\) 从哪里来，为什么少了它会偏高。 设 \(X=\left(\hat\mu-\frac{\sigma^2}{2}\right)T+\sigma\sqrt T\epsilon\)，它是正态，均值 \(m=(\hat\mu-\sigma^2/2)T\)，方差 \(v=\sigma^2T\)。 对数正态的期望公式：\(E[e^X]=e^{m+v/2}\)。代入得 \(E[S(T)]=S_0e^{(\hat\mu-\sigma^2/2)T+\sigma^2T/2}=S_0e^{\hat\mu T}\)，正好是风险中性要求的增长。 如果去掉 \(-\sigma^2/2\)，\(m=\hat\mu T\)，期望就变成 \(S_0e^{\hat\mu T+\sigma^2T/2}\)，多出 \(e^{\sigma^2T/2}\)。以 \(\sigma=0.3,T=0.5\) 为例，多出约 2.3%。 直觉：这就是算术平均收益与几何平均收益之差。CFA 里讲过"几何平均 ≈ 算术平均 − 方差/2"，波动越大，二者差距越大。\(\hat\mu\) 是算术意义上的期望增长率，\(\hat\mu-\sigma^2/2\) 是对数（几何）意义上的增长率。

优点与缺点

优点：

• 收益依赖路径（如 0 到 \(T\) 的平均价格）或只依赖终值都可以处理；
• 收益可以在期内多次发生，不必只在到期日；
• 可以容纳几乎任何随机过程（跳跃、随机波动率、均值回复）；
• 可以推广到多个标的变量；
• 自动给出估计的标准误。

缺点：

• 计算非常耗时；
• 难以处理提前行权：模拟从前往后走，而行权决策需要知道未来的继续持有价值（第 27b 章介绍若干扩展方法，如最小二乘蒙特卡洛）。

Business Snapshot 21.2：用 Excel 检验 BSM

参数：\(S_0=K=50\)，\(r=0.05\)，\(\sigma=0.3\)，\(T=0.5\)。\(d_1=0.2239\)，\(d_2=0.0118\)，BSM 看涨价为 4.817（原书表 21.2）。

• A1 = $C$2*EXP(($E$2−$F$2*$F$2/2)*$G$2+$F$2*NORMSINV(RAND())*SQRT($G$2))：这就是 (21.17)，NORMSINV(RAND()) 用逆变换法生成标准正态样本；
• B1 = EXP(−$E$2*$G$2)*MAX(A1−$D$2,0)：收益的现值；
• 复制 1,000 行，B 列均值 4.98，标准差 7.68。

依赖多个市场变量的衍生品

收益依赖 \(n\) 个变量 \(\theta_i\)（\(i=1,\dots,n\)）。设 \(s_i\) 为 \(\theta_i\) 的波动率，\(\hat m_i\) 为风险中性期望增长率，\(\rho_{ik}\) 为驱动 \(\theta_i\) 和 \(\theta_k\) 的维纳过程之间的相关系数（这些量都可以依赖于 \(\theta_i\)）。离散过程为

每次试验从多元标准正态分布中抽取 \(N\) 组 \((\epsilon_1,\dots,\epsilon_n)\)，生成每个变量的路径，再计算衍生品的一个样本值。关键是如何生成相关的正态样本。

生成相关的正态随机数

两个变量。要得到相关系数为 \(\rho\) 的两个标准正态样本，先取两个独立标准正态 \(x_1,x_2\)，令

验证：\(\text{Var}(\epsilon_2)=\rho^2+(1-\rho^2)=1\)，\(\text{Cov}(\epsilon_1,\epsilon_2)=\rho\)。

\(n\) 个变量。设所需相关系数为 \(\rho_{ij}\)，令

逐个确定系数，使方差为 1、相关系数正确：

用矩阵语言，若 \(\boldsymbol\epsilon=L\mathbf x\)，\(L\) 为下三角矩阵，则 \(\text{Cov}(\boldsymbol\epsilon)=LL^{\mathsf T}\)。上面的方程组就是要求 \(LL^{\mathsf T}=R\)（相关矩阵）。这正是 Cholesky 分解（Cholesky decomposition）（矩阵理论见第 01 册第 07a 章）。

推导拆解：为什么 \(\text{Cov}(\boldsymbol\epsilon)=LL^{\mathsf T}\)，以及方程组怎么逐个解。 第一步，\(\boldsymbol\epsilon=L\mathbf x\) 的第 \(i\) 个分量是 \(\epsilon_i=\sum_k\alpha_{ik}x_k\)。由于 \(x_k\) 相互独立、方差为 1，\(\text{Cov}(\epsilon_i,\epsilon_j)=\sum_k\alpha_{ik}\alpha_{jk}\)（不同 \(k\) 的交叉项协方差为零）。这个求和正是矩阵 \(LL^{\mathsf T}\) 第 \(i\) 行第 \(j\) 列的元素。 第二步，"下三角"的好处是可以按顺序解：\(\epsilon_1\) 只用 \(x_1\)；\(\epsilon_2\) 用 \(x_1,x_2\)，先由与 \(\epsilon_1\) 的相关确定 \(\alpha_{21}\)，再由方差为 1 确定 \(\alpha_{22}=\sqrt{1-\alpha_{21}^2}\)；以此类推。每一行的最后一个系数都要开平方：\(\alpha_{ii}=\sqrt{1-\sum_{k<i}\alpha_{ik}^2}\)。 金融直觉：可以把 \(x_1,x_2,\dots\) 看作相互独立的"风险因子"。\(\epsilon_1\) 完全暴露于因子 1；\(\epsilon_2\) 一部分暴露于因子 1（这部分带来与 \(\epsilon_1\) 的相关），剩下的由自己特有的因子 2 补足方差。这与 CFA 里"系统性风险 + 特有风险"的分解很像。

如果方程组没有实数解（某一步要对负数开平方），说明假设的相关结构内部不一致，即相关矩阵不是正定的（第 23 章 23.7 节）。这在实务中很常见：分别估计的两两相关系数、或者人为设定的压力情景相关矩阵，组合起来可能不满足正定性。

金融直觉：相关矩阵"正定"的意思是：用这些资产构成的任何非零组合，算出来的方差都大于零。如果相关矩阵不正定，就存在某个组合，按这个矩阵算出的方差是负数，这在现实中不可能。练习 4 的矩阵就是例子：资产 1 与 2、与 3 都高度正相关，按传递性 2 与 3 也应正相关，可矩阵却说它们强负相关。Cholesky 在某一步要对负数开方，就是这个矛盾在计算上的表现。

试验次数与标准误

设 \(M\) 次试验得到的折现收益样本均值为 \(\mu\)、标准差为 \(\omega\)。由中心极限定理，估计的**标准误（standard error）**为 \(\omega/\sqrt M\)，衍生品价值 \(f\) 的 95% 置信区间为

估值的不确定性与试验次数的平方根成反比：精度翻倍需要 4 倍试验次数，精度提高 10 倍需要 100 倍试验次数。

例 21.7：\(\pi\) 的估计，\(\omega=1.69\)，\(M=100\)，标准误 0.169，95% 置信区间 2.71–3.37，包含真值。

例 21.8：表 21.2 中 \(\omega=7.68\)，\(M=1000\)，标准误 \(7.68/\sqrt{1000}=0.24\)，95% 置信区间 4.51–5.45，包含 BSM 价格 4.817。

在树上抽样

不一定要从连续随机过程抽样，也可以在 \(N\) 步二叉树的 \(2^N\) 条路径中随机抽样。例如上升概率 0.6 时，在每个节点抽一个 0–1 均匀随机数，小于 0.4 走下支、大于 0.4 走上支。得到完整路径后计算收益，重复多次，平均收益按无风险利率贴现（Mintz 1997 讨论了高效实现）。

例 21.9：用图 21.3 的树（\(S_0=50\)，\(u=1.1224\)，\(d=0.8909\)，\(p=0.5073\)，5 步）为收益为 \(\max(S_{ave}-50,0)\) 的期权估值，\(S_{ave}\) 是 5 个月的平均股价（包括首尾共 6 个价格）。这是一种亚式期权（Asian option）。10 次试验（原书表 21.3）：

平均收益 7.08，价值 \(7.08e^{-0.1\times5/12}=6.79\)。这只是示意，实际应用需要多得多的步数和试验次数。注意：重合树对路径依赖期权无法直接倒推（节点丢失了路径信息），但在树上抽样路径没有这个问题。

计算希腊字母

对变量 \(x\)（标的价格或某个参数）的偏导数：先用蒙特卡洛得到 \(\hat f\)，再把 \(x\) 增加 \(\Delta x\) 重新模拟得到 \(\hat f^*\)，估计为 \((\hat f^*-\hat f)/\Delta x\)。

为了减小这个差分估计的标准误，两次计算应使用相同的时间区间数 \(N\)、相同的随机数样本和相同的试验次数 \(M\)——这称为公共随机数（common random numbers）。道理是：\(\hat f\) 和 \(\hat f^*\) 各自都有模拟误差，如果用独立的随机数，差分的误差是两个误差之差，再除以一个小的 \(\Delta x\)，会被极度放大；用相同的随机数，两个误差高度正相关，在相减时大部分抵消。本章"量化实战"会看到，二者的标准误相差近 20 倍。

推导拆解：用方差公式看清放大效应。 差分估计是 \((\hat f^*-\hat f)/\Delta x\)，其方差为 \(\frac{1}{\Delta x^2}[\text{Var}(\hat f^*)+\text{Var}(\hat f)-2\text{Cov}(\hat f^*,\hat f)]\)。 独立随机数时协方差为零，方差约为 \(2\text{Var}(\hat f)/\Delta x^2\)。\(\Delta x\) 越小，除以 \(\Delta x^2\) 后放大得越厉害。 公共随机数时，两次模拟只差一个很小的参数扰动，同一条路径上的两个收益几乎相同，相关系数接近 1，中括号里的三项几乎相互抵消。 这和配对交易的逻辑一样：两个高度相关的头寸单独看都波动很大，但价差的波动很小。

应用范围

当有三个或更多随机变量时，蒙特卡洛通常比其他方法更有效：它的计算时间大致随变量数线性增长，而树和有限差分随变量数指数增长。此外，它能给出估计的标准误，能容纳复杂的收益和随机过程，收益可以是整条路径的任意函数。

21.7 方差缩减方法

按 (21.13)–(21.18) 直接模拟，通常需要极多次试验才能达到合理的精度。由于标准误为 \(\omega/\sqrt M\)，降低误差有两条路：增加 \(M\)，或者减小每次试验的方差 \(\omega^2\)。后者就是方差缩减（variance reduction），往往能以零成本换来数倍乃至数十倍的效率提升。

对偶变量技术

对偶变量技术（antithetic variable technique）：每次试验计算两个值。\(f_1\) 按常规计算；\(f_2\) 把所有标准正态样本取相反号（\(\epsilon\to-\epsilon\)）后计算。试验的样本值取二者平均：

当一个值高于真值时，另一个往往低于真值，平均后误差减小。最终估计为 \(\bar f\) 的均值，标准误为 \(\omega/\sqrt M\)，其中 \(\omega\) 是 \(\bar f\) 的标准差、\(M\) 是配对数。这个标准误通常远小于用 \(2M\) 次独立试验得到的标准误。

为什么有效：\(\text{Var}(\bar f)=\frac14[\text{Var}(f_1)+\text{Var}(f_2)+2\text{Cov}(f_1,f_2)]\)。若收益是 \(\epsilon\) 的单调函数（如普通看涨期权），\(f_1\) 与 \(f_2\) 负相关，协方差项为负，方差下降。若收益是 \(\epsilon\) 的对称函数（如跨式组合），对偶变量就几乎没有帮助。

控制变量技术

控制变量技术（control variate technique）：待估衍生品 A 没有解析解，另一个相似的衍生品 B 有解析解 \(f_B\)。用相同的随机数流和相同的 \(\Delta t\) 并行模拟两者，得到 \(f_A^*\) 和 \(f_B^*\)，改进的估计为

这和 21a 章树方法中的控制变量是同一思想：模拟误差 \(f_B^*-f_B\) 可以精确算出，假设 A 的模拟误差与之相同并扣除。

Hull & White（1987）在研究随机波动率对欧式看涨价格的影响时使用了这一方法：A 是随机波动率下的期权，B 是常数波动率下的期权（BSM 有解析解）。

推广：(21.20) 相当于假设 A 与 B 的误差系数为 1。更一般地可以用 \(f_A^*-b(f_B^*-f_B)\)，其中最优的 \(b=\text{Cov}(f_A^*,f_B^*)/\text{Var}(f_B^*)\)，可以从样本中估计。这时方差缩减的比例为 \(1-\rho_{AB}^2\)，其中 \(\rho_{AB}\) 是两个样本的相关系数。控制变量与目标越相关，效果越好。

金融直觉：最优 \(b\) 的公式你一定眼熟：\(\text{Cov}/\text{Var}\) 就是回归斜率，也就是 CAPM 里 beta 的定义，或者最小方差套期保值比率。控制变量本质上就是"用 B 给 A 做最小方差对冲"：把 A 的模拟值对 B 的模拟值回归，只保留 B 解释不了的残差部分。方差缩减比例 \(1-\rho^2\) 对应回归的"未解释部分" \(1-R^2\)。例如 \(\rho=0.9\) 时方差剩 19%，相当于模拟次数放大约 5 倍；\(\rho=0.99\) 时剩约 2%，相当于放大约 50 倍。

常见的选择：用几何平均亚式期权（有解析解）作为算术平均亚式期权的控制变量；用标的资产自身 \(e^{-rT}S_T\)（期望精确为 \(S_0e^{-qT}\)）作为任何期权的控制变量。

重要性抽样

为深度虚值欧式看涨定价时，常规抽样的大部分路径收益为零，这些计算都浪费了。**重要性抽样（importance sampling）**把抽样集中在"重要"的区域。

设 \(F\) 为 \(T\) 时刻股价的无条件分布，\(q\) 为 \(S_T>K\) 的概率，且 \(q\) 可以解析算出。则 \(G=F/q\) 是在 \(S_T>K\) 条件下的股价分布。直接从 \(G\) 中抽样，每条路径都有正收益；期权价值的估计为

道理是 \(E[\text{收益}]=E[\text{收益}\mid S_T>K]\cdot P(S_T>K)\)，而 \(S_T\le K\) 时收益为零。

对 GBM 而言，\(S_T>K\) 等价于 \(\epsilon>z^*\)，其中 \(z^*=[\ln(K/S_0)-(\hat\mu-\sigma^2/2)T]/(\sigma\sqrt T)\)，\(q=1-N(z^*)\)。从截尾正态 \(\epsilon\mid\epsilon>z^*\) 中抽样，可以用逆变换法：\(\epsilon=N^{-1}(1-uq)\)，\(u\) 为 \([0,1]\) 均匀随机数。

推导拆解：两个式子各一步。 \(z^*\) 的来源：\(S_T>K\) ⇔ \(\ln S_0+(\hat\mu-\sigma^2/2)T+\sigma\sqrt T\epsilon>\ln K\)，两边移项、除以 \(\sigma\sqrt T\) 就得到 \(\epsilon>z^*\)。它和 BSM 里的 \(-d_2\) 是同一个量（当 \(\hat\mu=r\) 时），所以 \(q=1-N(z^*)=N(d_2)\)，即风险中性下期权到期实值的概率。 逆变换公式：\(u\) 在 \([0,1]\) 均匀，\(1-uq\) 就在 \([1-q,1]\) 上均匀，而 \(1-q=N(z^*)\)。\(N^{-1}\) 把 \([N(z^*),1]\) 映回 \([z^*,\infty)\)，所以抽出的 \(\epsilon\) 总是大于 \(z^*\)，并且在这个区间内的分布正好是截尾正态。 这里的 \(q\) 是概率，不是前文的股息率，注意区分。

分层抽样

分层抽样（stratified sampling）：不随机抽取，而是抽取代表值。要取 1,000 个样本，就把分布分成 1,000 个等概率区间，每个区间取一个代表值（通常是区间的均值或中位数）。

对标准正态分布分成 \(n\) 个区间时，第 \(i\) 个代表值可取

例如 \(n=4\) 时为 \(N^{-1}(0.125)\)、\(N^{-1}(0.375)\)、\(N^{-1}(0.625)\)、\(N^{-1}(0.875)\)。\(N^{-1}\) 在 Excel 中是 NORMSINV。

分层抽样消除了"样本恰好在某个区域扎堆"的随机性，对单维问题（如只依赖 \(S_T\) 的欧式期权）效果极好。缺点是必须事先确定样本数，且在高维下分层数会指数增长。

矩匹配

矩匹配（moment matching）：调整标准正态样本，使样本的一阶、二阶（甚至更高阶）矩与理论值精确匹配。设样本 \(\epsilon_i\)（\(i=1,\dots,n\)）的均值为 \(m\)、标准差为 \(s\)，调整后的样本为

其均值恰为 0、标准差恰为 1，在所有计算中用 \(\epsilon_i^*\) 代替 \(\epsilon_i\)。

矩匹配能节省计算时间，但可能带来内存问题：所有样本须保存到模拟结束才能算出 \(m\) 和 \(s\)。它又称二次重抽样（quadratic resampling），常与对偶变量结合使用：对偶变量自动匹配所有奇数阶矩（样本关于 0 对称），矩匹配只需再匹配二阶、可能还有四阶矩。

准随机序列

准随机序列（quasi-random sequences），又称低差异序列（low-discrepancy sequences），是概率分布的代表性样本序列。"准随机"这个名字有误导性——序列是完全确定性的。

用准随机序列可以让标准误与 \(1/M\) 而非 \(1/\sqrt M\) 成正比（在维度不太高、被积函数足够光滑时）。这是数量级的改进：\(M=10^4\) 时，前者误差约 \(10^{-4}\)，后者约 \(10^{-2}\)。

白话解释：伪随机数是"独立地乱撒"，难免有的地方扎堆、有的地方空着，误差来自这种不均匀，按 \(1/\sqrt M\) 缓慢消失。准随机序列是"有计划地铺开"，每个新点都放在现有点最稀疏的地方，所以覆盖得均匀很多。严格的理论结果是误差约为 \((\ln M)^d/M\)（\(d\) 为维度），维度低时 \((\ln M)^d\) 增长很慢，近似看作 \(1/M\)；维度高时这个因子会变大，这就是练习 10 要观察的现象。正文说"误差约 \(10^{-4}\)"是忽略常数和对数因子后的量级示意。

准随机序列与分层抽样相似，但更灵活。分层抽样要事先知道样本数；准随机序列的每个新样本都"填补"已有样本之间的空隙，所以在任何阶段，样本点都大致均匀地分布在概率空间中。原书图 21.14 展示了二维 Sobol' 序列的前 1,024 个点，点与点之间分布非常均匀，没有伪随机数常见的空洞和聚集。实现可参考 Brotherton-Ratcliffe（1994）和 Press 等的《Numerical Recipes》；Python 中 scipy.stats.qmc 提供了 Sobol' 序列。

准随机序列是确定性的，无法直接给出标准误。实务中常用随机化（scrambled）准随机序列：用不同随机种子生成若干组序列，用各组估计的离散程度衡量误差。

方法选择小结

量化实战

本章内容在量化交易中的用途

1. 奇异期权与结构化产品定价。亚式、回望、障碍、一篮子期权、雪球等自动赎回结构，几乎都靠蒙特卡洛定价。生产系统中控制变量、准随机序列和对偶变量是标配。
2. 风险模型。第 22 章的蒙特卡洛 VaR 需要从多元正态（或其他分布）中抽取相关的风险因子变动，Cholesky 分解是核心步骤；相关矩阵非正定时要先修正（如把负特征值截断为零后重新标准化）。
3. 策略研究与回测。用 (21.16) 模拟 GBM 路径，可以检验一个期权策略或仓位规则在不同市场环境下的收益分布；用分层或准随机抽样可以用更少的路径得到稳定的结论。
4. 判断模拟是否充分。标准误 \(\omega/\sqrt M\) 告诉你结果的精度。量化研究中同样的逻辑适用于回测：一个策略的平均日收益 \(\bar r\)、标准差 \(s\)、天数 \(T\)，其标准误是 \(s/\sqrt T\)，用它判断 Sharpe 比率是否显著（第 03 册第 10a、10b 章）。
5. 敏感度计算。用公共随机数做差分是蒙特卡洛希腊字母的最简方法；进阶方法有路径导数（pathwise）和似然比（likelihood ratio）方法，见 Glasserman 的专著（本册不展开）。

Python 示例：六种方差缩减的对比

下面以 Business Snapshot 21.2 的参数（\(S_0=K=50\)，\(r=5\%\)，\(\sigma=30\%\)，\(T=0.5\)，BSM 价格 4.817）为例，对比普通蒙特卡洛与各种方差缩减方法，展示 \(1/\sqrt M\) 的收敛速度，用重要性抽样为深度虚值期权定价，用 Cholesky 分解生成相关随机数，用公共随机数估计 delta，最后复算例 21.9。

import numpy as np
from scipy.stats import norm, qmc

rng = np.random.default_rng(2024)
S0, K, r, sig, T = 50.0, 50.0, 0.05, 0.30, 0.5          # Business Snapshot 21.2 参数

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

def summary(name, x):
    m, se = x.mean(), x.std(ddof=1)/np.sqrt(len(x))
    print(f"{name:<22s} 估计 {m:.4f}  标准误 {se:.4f}  95%CI [{m-1.96*se:.3f}, {m+1.96*se:.3f}]")
    return se

ST = lambda z: S0*np.exp((r - 0.5*sig**2)*T + sig*np.sqrt(T)*z)   # 式 (21.17)
disc = np.exp(-r*T)
M = 10_000
print("BSM 解析价 %.4f" % exact)

z = rng.standard_normal(M)
se0 = summary("普通蒙特卡洛", disc*np.maximum(ST(z)-K, 0))

# 对偶变量：M/2 对，总共仍用 M 个正态数
zh = z[:M//2]
pair = 0.5*disc*(np.maximum(ST(zh)-K, 0) + np.maximum(ST(-zh)-K, 0))
summary("对偶变量", pair)

# 控制变量：以 e^{-rT} S_T（真值 S0）为控制，b 取最优回归系数
fA = disc*np.maximum(ST(z)-K, 0); fB = disc*ST(z)
b = np.cov(fA, fB)[0, 1]/fB.var(ddof=1)
summary("控制变量(S_T)", fA - b*(fB - S0))

# 分层抽样：M 个等概率区间取中点分位数
zs = norm.ppf((np.arange(M) + 0.5)/M)
print(f"{'分层抽样':<22s} 估计 {np.mean(disc*np.maximum(ST(zs)-K,0)):.4f}")

# 矩匹配
zm = (z - z.mean())/z.std()
print(f"{'矩匹配':<22s} 估计 {np.mean(disc*np.maximum(ST(zm)-K,0)):.4f}")

# 准随机 Sobol' 序列（随机化以便估计误差）
ests = []
for s in range(20):
    u = qmc.Sobol(d=1, scramble=True, seed=s).random(2**13).ravel()
    ests.append(np.mean(disc*np.maximum(ST(norm.ppf(u))-K, 0)))
print(f"{'Sobol(8192点×20次)':<22s} 估计 {np.mean(ests):.4f}  各次标准差 {np.std(ests, ddof=1):.4f}")

# 收敛速度：误差 ∝ 1/sqrt(M)
print("\n普通MC 标准误随 M:")
for m in (1_000, 4_000, 16_000, 64_000):
    x = disc*np.maximum(ST(rng.standard_normal(m))-K, 0)
    print(f"  M={m:6d}  标准误 {x.std(ddof=1)/np.sqrt(m):.4f}")

# 重要性抽样：深度虚值看涨 K=80
K2 = 80.0
exact2 = bs_call(S0, K2, T, r, sig)
x_naive = disc*np.maximum(ST(rng.standard_normal(M))-K2, 0)
zstar = (np.log(K2/S0) - (r-0.5*sig**2)*T)/(sig*np.sqrt(T))   # S_T>K2 <=> z>z*
qprob = norm.sf(zstar)
u = rng.uniform(size=M)
zt = norm.isf(u*qprob)                                         # 从 z | z>z* 抽样
x_is = qprob*disc*(ST(zt) - K2)
print("\n深度虚值看涨 K=80  解析 %.5f" % exact2)
summary("  普通MC", x_naive); summary("  重要性抽样", x_is)

# Cholesky：三个相关正态
C = np.array([[1, 0.6, 0.3], [0.6, 1, 0.5], [0.3, 0.5, 1]])
L = np.linalg.cholesky(C)
eps = rng.standard_normal((100_000, 3)) @ L.T
print("\nCholesky L=\n", np.round(L, 4), "\n样本相关=\n", np.round(np.corrcoef(eps.T), 3))

# 公共随机数估计 delta
h = 0.5
zc = rng.standard_normal(M)
f = lambda s0, zz: disc*np.maximum(s0*np.exp((r-0.5*sig**2)*T+sig*np.sqrt(T)*zz)-K, 0)
d_crn = (f(S0+h, zc) - f(S0-h, zc))/(2*h)
d_ind = (f(S0+h, rng.standard_normal(M)) - f(S0-h, rng.standard_normal(M)))/(2*h)
d1 = (np.log(S0/K)+(r+0.5*sig**2)*T)/(sig*np.sqrt(T))
print("\ndelta 解析 %.4f | 公共随机数 %.4f (标准误 %.4f) | 独立随机数 %.4f (标准误 %.4f)" %
      (norm.cdf(d1), d_crn.mean(), d_crn.std(ddof=1)/np.sqrt(M), d_ind.mean(), d_ind.std(ddof=1)/np.sqrt(M)))

# 例 21.9：在图 21.3 的树上抽样，亚式期权
u_, d_ = 1.1224, 0.8909
paths = ["UUUUD","UUUDD","DDDUU","UUUUU","UUDDU","UDUUD","DDUDD","UUDDU","UUUDU","DDUUD"]
pays = []
for pth in paths:
    s = [50.0]
    for c in pth: s.append(s[-1]*(u_ if c == "U" else d_))
    pays.append(max(np.mean(s) - 50, 0))
print("\n例21.9 平均收益 %.2f  价值 %.2f" % (np.mean(pays), np.mean(pays)*np.exp(-0.1*5/12)))

关键输出：

BSM 解析价 4.8174
普通蒙特卡洛                 估计 4.8043  标准误 0.0737  95%CI [4.660, 4.949]
对偶变量                   估计 4.7504  标准误 0.0569  95%CI [4.639, 4.862]
控制变量(S_T)              估计 4.7650  标准误 0.0308  95%CI [4.705, 4.825]
分层抽样                   估计 4.8172
矩匹配                    估计 4.8024
Sobol(8192点×20次)       估计 4.8175  各次标准差 0.0009

普通MC 标准误随 M:
  M=  1000  标准误 0.2336
  M=  4000  标准误 0.1149
  M= 16000  标准误 0.0575
  M= 64000  标准误 0.0294

深度虚值看涨 K=80  解析 0.08585
  普通MC                 估计 0.0944  标准误 0.0105  95%CI [0.074, 0.115]
  重要性抽样                估计 0.0866  标准误 0.0008  95%CI [0.085, 0.088]

Cholesky L=
 [[1.    0.    0.   ]
 [0.6   0.8   0.   ]
 [0.3   0.4   0.866]] 
样本相关=
 [[1.    0.601 0.303]
 [0.601 1.    0.503]
 [0.303 0.503 1.   ]]

delta 解析 0.5886 | 公共随机数 0.5828 (标准误 0.0059) | 独立随机数 0.5779 (标准误 0.1047)

例21.9 平均收益 7.07  价值 6.79

读输出：

• 收敛速度：\(M\) 每乘 4，标准误约减半（0.234 → 0.115 → 0.058 → 0.029），印证了 \(1/\sqrt M\) 规律。
• 对偶变量：用同样多的正态随机数（5,000 对），标准误从 0.074 降到 0.057，方差约降 40%。平值看涨的收益对 \(\epsilon\) 单调但高度非线性（一半路径收益为 0），所以改进有限。
• 控制变量：以 \(e^{-rT}S_T\) 为控制变量，标准误降到 0.031，相当于方差降低约 83%，即把试验次数放大约 6 倍。
• 分层抽样与 Sobol' 序列：一维问题上几乎是精确的。20 组随机化 Sobol' 估计的标准差只有 0.0009，比同样点数的普通蒙特卡洛（约 0.08）小两个数量级，体现了接近 \(1/M\) 的收敛。
• 重要性抽样：执行价 80 的深度虚值看涨，普通蒙特卡洛 1 万条路径中只有极少数有正收益，标准误 0.0105（相对误差 12%）；重要性抽样标准误 0.0008，降低了一个数量级以上。
• Cholesky：\(L\) 的第二行 \((0.6,0.8)\) 正是两变量公式 \((\rho,\sqrt{1-\rho^2})\)；样本相关矩阵与目标一致。
• 公共随机数：同样 1 万条路径，独立随机数的 delta 标准误是公共随机数的近 18 倍，估计几乎没有用处。
• 例 21.9：用原书四位小数的 \(u,d\) 复算得平均收益 7.07，与原书 7.08 的差异来自舍入；贴现后价值同为 6.79。

本章小结

蒙特卡洛模拟在风险中性世界中从今天向前随机生成大量路径，计算每条路径的收益，平均后按无风险利率贴现得到衍生品价值。模拟 \(\ln S\) 的离散化 (21.16) 在参数为常数时对任意步长精确，只依赖终值的收益可以用 (21.17) 一步到位。多个标的变量时，用 Cholesky 分解从独立正态样本构造相关样本；无解意味着相关矩阵非正定。估计的标准误为 \(\omega/\sqrt M\)，精度提高 10 倍需要 100 倍的试验。也可以在二叉树上抽样路径，这让路径依赖期权可以用树的参数估值。希腊字母用扰动参数后的差分估计，必须使用公共随机数。蒙特卡洛擅长路径依赖和高维问题，不擅长提前行权。方差缩减技术——对偶变量、控制变量、重要性抽样、分层抽样、矩匹配、准随机序列——通过降低单次试验的方差或使样本更均匀，可以把效率提高数倍到数百倍。

练习

基础

1. 用 Business Snapshot 21.1 的方法估计 \(\pi\)，\(M=100\) 时标准误约 0.169。要使标准误降到 0.01，需要多少次试验？ 提示：\(1.69/\sqrt M=0.01\)，\(M\approx28{,}561\)。

2. 用分层抽样改进 \(\pi\) 的估计：把正方形分成 \(10\times10\) 个小格，每格取中心点。结果是多少？与 100 次随机投掷相比如何？ 提示：只有一个确定的估计，没有随机误差；误差来自"格子被圆周切开"的部分（原书习题 21.8）。

3. 写出三个相关标准正态样本 \(\epsilon_1,\epsilon_2,\epsilon_3\) 由独立样本 \(x_1,x_2,x_3\) 生成的公式（相关系数 \(\rho_{12},\rho_{13},\rho_{23}\)）。 提示：\(\epsilon_1=x_1\)；\(\epsilon_2=\rho_{12}x_1+\sqrt{1-\rho_{12}^2}x_2\)；\(\epsilon_3=\alpha_{31}x_1+\alpha_{32}x_2+\alpha_{33}x_3\)，\(\alpha_{31}=\rho_{13}\)，\(\alpha_{32}=(\rho_{23}-\rho_{12}\rho_{13})/\sqrt{1-\rho_{12}^2}\)，\(\alpha_{33}=\sqrt{1-\alpha_{31}^2-\alpha_{32}^2}\)（原书习题 21.24）。

4. 若相关矩阵为 \(\begin{pmatrix}1&0.9&0.9\\0.9&1&-0.9\\0.9&-0.9&1\end{pmatrix}\)，用上题的公式会发生什么？解释原因。 提示：\(\alpha_{32}=(-0.9-0.81)/0.436=-3.92\)，\(\alpha_{33}^2<0\)，无实数解。变量 1 与 2、3 都强正相关，2 与 3 却强负相关，内部不一致，矩阵非正定。

5. 为什么蒙特卡洛难以处理美式期权？ 提示：模拟从前往后，在某个时点决定是否行权需要知道"继续持有"的价值，而这依赖于未来所有路径的条件期望，单条路径上无法得到（原书习题 21.9）。

进阶

6. 在表 21.2 的基础上使用对偶变量。具体说明需要改哪几列、新的估计和标准误如何计算。 提示：增加一列用 \(-\epsilon\) 生成 \(S_T\)（需要先把 NORMSINV(RAND()) 存在单独一列），再计算两者收益现值的平均，对平均列求均值与标准差（原书习题 21.22）。

7. 用随机波动率模型为欧式看涨定价时，若同时使用对偶变量和控制变量（控制变量为常数波动率下的期权），解释为什么每次试验需要计算六个期权价值。 提示：随机波动率模型有两组随机数：驱动股价的 \(\epsilon_1\) 和驱动波动率的 \(\epsilon_2\)。对偶变量要对两者分别取反号，组合 \((\pm\epsilon_1,\pm\epsilon_2)\) 给出 4 个随机波动率期权（A）的值；常数波动率期权（B）只依赖 \(\epsilon_1\)，只有 \(\pm\epsilon_1\) 两个值。合计 6 个（原书习题 21.16）。

8. 证明控制变量 \(f_A^*-b(f_B^*-f_B)\) 的方差在 \(b^*=\text{Cov}(f_A^*,f_B^*)/\text{Var}(f_B^*)\) 时最小，最小方差为 \((1-\rho^2)\text{Var}(f_A^*)\)。原书 (21.20) 对应 \(b\) 取何值？什么时候这不是最优选择？ 提示：对 \(b\) 求导；(21.20) 对应 \(b=1\)，当 A 与 B 的波动幅度不同（如 A 的价值远大于 B）时 \(b=1\) 不是最优。

9. 用蒙特卡洛为算术平均亚式看涨期权定价（\(S_0=K=50\)，\(r=5\%\)，\(\sigma=30\%\)，\(T=1\)，按月平均 12 个点），分别用：(a) 普通方法；(b) 几何平均亚式期权作为控制变量。比较标准误。 提示：离散几何平均 \(G=(\prod S_{t_i})^{1/n}\) 的对数服从正态，可以写出其看涨期权的解析解；算术与几何平均高度相关（通常 \(\rho>0.99\)），方差缩减非常显著。

10. 用本章代码的思路，比较在 \(d=10\) 维（10 个独立资产等权一篮子）问题上，普通蒙特卡洛与 Sobol' 序列的误差。准随机序列的优势在高维下还明显吗？ 提示：可以把一篮子看涨期权作为测试；随着维度上升，低差异序列的优势会减弱，但对"有效维度"低的问题（收益主要由少数方向决定）仍然显著。

原书推荐习题：21.22、21.8（对偶变量与分层抽样的方差缩减效果）；21.24（Cholesky 分解）；21.9（蒙特卡洛为何难以处理美式期权）；21.16（随机波动率下对偶变量 + 控制变量）。

原书对照

树方法见第 21a 章，有限差分法与三种方法的比较见第 21c 章。延伸阅读：Boyle (1977) 蒙特卡洛期权定价的开创性论文；Boyle, Broadie & Glasserman (1997) 证券定价中的蒙特卡洛方法；Broadie, Glasserman & Jain (1997) 用随机树为美式期权定价；Brotherton-Ratcliffe (1994) 准随机序列；Press 等《Numerical Recipes》。