第 10 章 随机模拟与方差缩减

本章对应 Ross 原书第 10 章「模拟」。蒙特卡洛模拟是量化工程师每天都在用的工具：期权定价、风险情景、策略压力测试、置换检验、bootstrap 都离不开它。本章讲清三件事：怎样从均匀随机数生成任意分布（逆变换法、拒绝法、正态专用算法、离散分布）；模拟估计的误差由什么决定；怎样在不改变期望的前提下降低方差（对偶变量、条件化、控制变量，以及习题中的重要性抽样）。

学习目标

1. 理解伪随机数与线性同余生成器的原理，会用 Fisher–Yates 算法生成随机排列。
2. 掌握逆变换法，能对指数、Weibull、几何等分布写出生成公式。
3. 掌握拒绝法的步骤、正确性证明与效率（平均迭代 \(c\) 次），会用指数拒绝法、Box–Muller 和 Marsaglia 极坐标法生成正态。
4. 会生成二项、Poisson、Gamma、卡方等常见分布。
5. 理解蒙特卡洛估计的均方误差 \(\mathrm{Var}(Y)/k\)，掌握对偶变量、条件化、控制变量三种方差缩减技术及其适用条件，能用于期权定价。

读前导读

这一章在解决什么问题。 很多期望没法用公式算出来：路径依赖期权的价格、复杂组合在压力情景下的损失分布、一个交易规则的回测统计量在零假设下的分布。蒙特卡洛的办法很朴素：按模型随机生成大量情景，在每个情景里算出结果，再取平均。第 08 章的大数定律保证平均值收敛到真实期望，CLT 告诉你误差有多大。

你在 CFA 里可能见过「蒙特卡洛 VaR」或「蒙特卡洛退休规划」，但没学过两件事，本章补上。第一，计算机只会生成 0 到 1 之间的均匀随机数，怎么把它变成指数、正态、Poisson 等分布的样本？第二，误差只按 \(1/\sqrt k\) 下降，精度提高 10 倍要多跑 100 倍，所以要靠「方差缩减」技巧，在样本量不变的情况下把误差压下来。控制变量法和你熟悉的最小方差对冲是同一个公式：用一个已知价格的资产，对冲掉估计中的噪声。

需要先想起来的数学。

• 反函数。 若 \(y=F(x)\)，反函数 \(x=F^{-1}(y)\) 就是「反过来解出 \(x\)」。例：\(F(x)=1-e^{-2x}\)，令 \(u=1-e^{-2x}\)，解得 \(x=-\frac12\ln(1-u)\)。见 第 00 册第 01 章 函数极限与连续。
• 对数与指数的运算。 \(\log(ab)=\log a+\log b\)，\(\log\prod U_i=\sum\log U_i\)；\(0<U<1\) 时 \(\log U<0\)，所以 \(-\log U>0\)。见 第 00 册第 04 章 级数与收敛。
• 积分作为期望、蒙特卡洛积分。 \(\int_0^1g(x)\,dx=E[g(U)]\)，\(U\sim U(0,1)\)。例：\(\int_0^1x^2dx=1/3\)，随机抽很多个 \(U\)，算 \(U^2\) 的平均，也会接近 1/3。见 第 00 册第 03 章 积分。
• 极坐标与雅可比。 平面点 \((x,y)\) 可写成 \((r\cos\theta,r\sin\theta)\)；变量替换时面积元 \(dx\,dy=r\,dr\,d\theta\)。Box–Muller 方法的依据就在这里。见 第 00 册第 05 章 多元微积分与优化。
• 取整记号。 \(\lfloor x\rfloor\) 是不超过 \(x\) 的最大整数（向下取整），\(\lceil x\rceil\) 是不小于 \(x\) 的最小整数（向上取整）。例：\(\lfloor2.7\rfloor=2\)，\(\lceil2.1\rceil=3\)。见 第 00 册第 08 章 读懂数学证明与符号。

怎么读这一章。 核心必读：10.2.1 逆变换法、10.2.2 拒绝法的步骤与效率（证明可后看）、10.4 全节，以及量化实战第 2 部分（期权定价）。10.1 的线性同余法和 Fisher–Yates 看懂思路即可。正态生成的三种方法（例 2c、2d 和 Marsaglia 极坐标法）了解原理就够，实际工作中直接调用库函数。10.3 离散分布可以快速浏览。10.4.4 后半段关于自测题 10.5 的讨论值得一读，它说明了「期望不存在时模拟会怎样失败」。

10.1 引言

动机 一副牌按某固定策略玩接龙（solitaire），赢的概率是多少？\(52!\) 种排列无法枚举，只能「做实验」：玩 \(n\) 局，令 \(X_i\) 为第 \(i\) 局获胜的示性变量，由强大数定律（第 08 章），以概率 1

随机数 模拟的起点是生成 \(U(0,1)\) 随机变量。计算机生成的是伪随机数（pseudorandom numbers）：由确定的递推产生、但统计上看起来像独立均匀样本的序列。经典方法是线性同余法：给定种子（seed）\(X_0\) 和正整数 \(a,c,m\)，

现代量化系统中，线性同余法已被 Mersenne Twister、PCG、Philox 等生成器取代（numpy 默认 PCG64），原理相同。对回测和研究而言更重要的是种子管理：固定种子保证结果可复现，并行计算时每个进程要用独立的随机流（numpy 的 SeedSequence.spawn），否则不同进程可能生成相关甚至相同的随机数。

例 1a（生成随机排列：Fisher–Yates / Knuth 洗牌）

1. 取任一排列 \(X(1),\dots,X(n)\)，如 \(X(i)=i\)；
2. 令 \(I=n\)；
3. 生成随机数 \(U\)，令 \(N=\lfloor IU\rfloor+1\)（在 \(1,\dots,I\) 上均匀）；
4. 交换 \(X(N)\) 与 \(X(I)\)；
5. \(I\leftarrow I-1\)，若 \(I>1\) 回到第 3 步；
6. 输出 \(X(1),\dots,X(n)\)。

例：\(n=4\)，依次得 \(N=3\) → (1,2,4,3)；\(N=1\) → (4,2,1,3)；\(N=1\) → (2,4,1,3)。复杂度 \(O(n)\)，所有 \(n!\) 个排列等可能。应用：随机化实验设计——把受试者随机分组，以消除分组偏差。在量化中，随机排列是置换检验（permutation test）的核心：打乱因子与收益的对应关系，看观测到的统计量在「无关联」的零假设下有多罕见。

10.2 连续随机变量的一般生成方法

10.2.1 逆变换法

命题 2.1 \(U\sim U(0,1)\)，\(F\) 为连续分布函数，令 \(Y=F^{-1}(U)\)，则 \(Y\) 的分布函数为 \(F\)。

证明：\(F\) 单调，所以 \(F^{-1}(U)\le a\iff U\le F(a)\)，\(P\{Y\le a\}=P\{U\le F(a)\}=F(a)\)。（第 05 章 5.7 节已作为「概率积分变换」介绍过这一结论，并在其量化实战中用它抽样违约时间；本节把它系统化为通用的生成方法。）

白话解释：把 \(U\) 想成一个随机的「分位数」。\(U=0.05\) 就取分布的 5% 分位点，\(U=0.5\) 就取中位数。因为 \(U\) 在 0 到 1 上均匀，落在「5% 分位点以下」的概率恰好是 5%，生成的样本自然就服从 \(F\)。你在 Excel 里写 NORM.INV(RAND(), μ, σ) 生成正态收益，用的就是这个方法，只不过正态的 \(F^{-1}\) 由软件用数值方法近似计算。 证明里的「\(\iff\)」：\(F\) 递增，所以对不等式两边同时作用 \(F\) 或 \(F^{-1}\)，方向不变。 例 2a 的数值检查：\(U=0.5\) 时 \(-\log0.5=0.693\)，正是 Exp(1) 的中位数 \(\ln2\)。

例 2a（指数分布） \(F(x)=1-e^{-x}\)，\(F^{-1}(u)=-\log(1-u)\)。由于 \(1-U\) 也是均匀的，\(-\log U\sim\) Exp(1)；\(-\frac1\lambda\log U\sim\) Exp\((\lambda)\)。

例 2b（Gamma\((n,\lambda)\)，\(n\) 为整数） \(n\) 个独立 Exp\((\lambda)\) 之和：

更多例子（原书习题）：Weibull 分布 \(F(t)=1-e^{-at^\beta}\)（习题 10.5），解 \(F(t)=u\) 得 \(t=\big(-\frac1a\log(1-u)\big)^{1/\beta}\)；给定失效率函数 \(\lambda(t)\) 时 \(F(t)=1-e^{-\int_0^t\lambda(s)ds}\)（习题 10.6），信用模型中按违约强度生成违约时间用的就是它；\(F(x)=x^n\) 可取 \(U^{1/n}\) 或 \(n\) 个均匀数的最大值（习题 10.7）；混合分布 \(pF_1+(1-p)F_2\) 可先以概率 \(p\) 选成分，再从该成分抽样（组合法 composition，习题 10.9）。

逆变换法的局限是需要 \(F^{-1}\) 的显式表达式；正态分布就没有，需要别的方法。

10.2.2 拒绝法

设已经能从密度 \(g\) 抽样，想从密度 \(f\) 抽样，且存在常数 \(c\) 使

拒绝法（rejection method，acceptance–rejection）：

1. 从 \(g\) 抽 \(Y\)，并生成随机数 \(U\)；
2. 若 \(U\le\dfrac{f(Y)}{cg(Y)}\)，令 \(X=Y\) 并停止；否则回到第 1 步。

命题 2.2 拒绝法得到的 \(X\) 密度为 \(f\)。

证明：令 \(K=P\{U\le f(Y)/cg(Y)\}\) 为一次迭代被接受的概率。

推导拆解： （1）最终输出的 \(X\) 是「某一轮被接受的 \(Y\)」，各轮相互独立、规则相同，所以 \(X\) 的分布等于「在一轮被接受的条件下 \(Y\) 的分布」。 （2）条件概率 = 联合概率 / 条件事件概率。联合概率 \(P\{Y\le x,\ U\le\frac{f(Y)}{cg(Y)}\}\)：先对 \(Y=y\) 条件化，\(U\) 均匀，所以 \(P\{U\le h\}=h\)；再对 \(y\) 按密度 \(g\) 积分，得 \(\int_{-\infty}^x\frac{f(y)}{cg(y)}g(y)\,dy\)。\(g\) 约掉，只剩 \(f/c\)。 （3）\(x\to\infty\) 时左边是 1，右边 \(\frac1{cK}\int f=\frac1{cK}\)，所以 \(K=1/c\)。 白话：在 \(g\) 的密度曲线下撒点，再把高出 \(f\) 曲线的点扔掉，剩下的点就恰好分布在 \(f\) 曲线下面。\(c\) 是把 \(g\) 放大到能完全盖住 \(f\) 所需的倍数，被扔掉的比例是 \(1-1/c\)。

效率：每次迭代独立地以概率 \(1/c\) 被接受，迭代次数服从几何分布，均值为 \(c\)。\(g\) 的形状越接近 \(f\)，\(c\) 越接近 1，效率越高。

例 2c（用指数分布拒绝生成正态） \(|Z|\) 的密度 \(f(x)=\frac{2}{\sqrt{2\pi}}e^{-x^2/2}\)（\(x>0\)），取 \(g(x)=e^{-x}\)：

进一步改进：\(U\le e^{-(Y-1)^2/2}\iff-\log U\ge(Y-1)^2/2\)，而 \(-\log U\) 本身就是 Exp(1)。所以等价于生成两个 Exp(1) 变量 \(Y_1,Y_2\)，若 \(Y_2\ge(Y_1-1)^2/2\) 则接受 \(Y_1\)。由无记忆性，接受时的超出量 \(Y_2-(Y_1-1)^2/2\) 又是一个独立的 Exp(1)，可以留给下一轮当 \(Y_1\) 用。这样平均每个正态只需约 1.64 个指数变量。

例 2d（Box–Muller 方法） 两个独立标准正态 \((X,Y)\) 的极坐标满足：\(R^2\sim\) Exp（均值 2），\(\Theta\sim U(0,2\pi)\)，两者独立（第 06b 章例 7b）。所以

推导拆解：把 \((X,Y)\) 的生成拆成两半。半径：\(R^2\) 是均值 2 的指数分布，由例 2a，\(R^2=-2\log U_1\)，所以 \(R=\sqrt{-2\log U_1}\)。角度：\(\Theta=2\pi U_2\) 在 \((0,2\pi)\) 上均匀。再用 \(X=R\cos\Theta\)、\(Y=R\sin\Theta\) 换回直角坐标。为什么 \(R^2\) 是指数分布？二维标准正态密度 \(\frac1{2\pi}e^{-(x^2+y^2)/2}\) 只依赖 \(r^2=x^2+y^2\)，换成极坐标后（面积元 \(r\,dr\,d\theta\)）可以分解成「只含 \(\theta\) 的常数」乘以「只含 \(r\) 的部分」，前者说明角度均匀且与半径独立；令 \(s=r^2\)，\(ds=2r\,dr\)，后者变成 \(\frac12e^{-s/2}\)，正是均值 2 的指数密度。

Marsaglia 极坐标法（polar method）避开三角函数：令 \(V_i=2U_i-1\)，\((V_1,V_2)\) 在正方形 \([-1,1]^2\) 上均匀；不断重抽直到 \(S=V_1^2+V_2^2\le1\)，此时 \((V_1,V_2)\) 在单位圆盘上均匀，其极坐标满足 \(R^2=S\sim U(0,1)\)、\(\Theta\sim U(0,2\pi)\)，相互独立（习题 10.13）。而 \(\cos\Theta=V_1/R\)、\(\sin\Theta=V_2/R\)，且 \(S\) 本身可以充当 Box–Muller 中的 \(U_1\)：

例 2e（卡方分布） \(Z_1^2+Z_2^2\sim\) Exp（速率 1/2），所以 \(\chi^2_{2k}\sim\) Gamma\((k,1/2)\)，可取 \(-2\log\big(\prod_{i=1}^kU_i\big)\)；奇数自由度 \(\chi^2_{2k+1}=Z^2-2\log\big(\prod_{i=1}^kU_i\big)\)。

10.3 离散分布的模拟

离散逆变换 \(P\{X=x_j\}=P_j\)。生成 \(U\)，若 \(\sum_{i<j}P_i<U\le\sum_{i\le j}P_i\)，令 \(X=x_j\)；其概率恰为 \(P_j\)。按概率从大到小排列取值再比较，可以减少平均比较次数（自测题 10.3）。

例 3a（几何分布） \(P\{X=i\}=(1-p)^{i-1}p\)，\(\sum_{i<j}P\{X=i\}=1-(1-p)^{j-1}\)。\(X=j\) 当且仅当 \((1-p)^j\le1-U<(1-p)^{j-1}\)；用 \(U\) 代替 \(1-U\)，注意 \(\log(1-p)<0\) 会使不等号反向：

例 3b（二项分布） \(X=\sum_{i=1}^nI(U_i<p)\sim\text{Bin}(n,p)\)。

例 3c（Poisson 分布） 不断生成随机数，令 \(N=\min\{n:\prod_{i=1}^nU_i<e^{-\lambda}\}\)，则 \(X=N-1\sim\) Poisson\((\lambda)\)。理由：

10.4 方差缩减技术

10.4.1 蒙特卡洛估计的误差

要估计 \(\theta=E[g(X_1,\dots,X_n)]\)。独立生成 \(k\) 组 \(\mathbf X^{(j)}\)，令 \(Y_j=g(\mathbf X^{(j)})\)，

金融直觉：数值感受：一个价格约 5.6 元的期权，单条路径收益的标准差约 10 元。跑 10 万条路径，标准误 \(10/\sqrt{100000}\approx0.03\) 元；想把标准误降到 0.003 元，要跑 1000 万条。方差缩减的价值可以直接折算成算力：方差降到 1/4，相当于免费多跑 4 倍路径。下面三种技术的共同原则是：期望不变（无偏），只换一个噪声更小的估计量。

10.4.2 对偶变量

两个同分布估计的平均：

\(g\) 不单调时（如跨式期权、对称的收益函数）对偶变量可能无效甚至有害。

金融直觉：对偶变量相当于「每抽一个上涨情景，就配一个对称的下跌情景」。对看涨期权这种单调收益，一个情景收益高，它的镜像情景收益就低，两者平均后噪声互相抵消一部分。跨式期权收益 \(|S_T-K|\) 关于平值附近大致对称，\(Z\) 和 \(-Z\) 两种情景的收益差不多大，正相关而不是负相关，平均后噪声反而没减少（练习 5）。 公式中 \(\mathrm{Cov}/2\) 一项是关键：\(\mathrm{Cov}=0\) 时就退化为两次独立模拟；\(\mathrm{Cov}<0\) 才有收益。量化实战里，看涨期权的方差只降到 1/1.42，因为收益函数在 \(S_T<K\) 一侧是平的，镜像情景经常两边都是 0，负相关不强。

10.4.3 条件化方差缩减

由全方差公式（第 07b 章）

例 4a（估计 \(\pi\)） \((V_1,V_2)\) 在 \([-1,1]^2\) 上均匀，\(I\) 为落在单位圆内的示性变量，\(E[I]=\pi/4\)，用 \(4\times\) 落入比例估计 \(\pi\)。

改进 1（条件化）：\(E[I\mid V_1]=P\{V_2^2\le1-V_1^2\mid V_1\}=\sqrt{1-V_1^2}\)；又 \(E\big[\sqrt{1-V_1^2}\big]=E\big[\sqrt{1-U^2}\big]\)，所以用 \(\sqrt{1-U^2}\) 的样本均值估计 \(\pi/4\)。

推导拆解：给定 \(V_1=v\)，点落在圆内要求 \(V_2^2\le1-v^2\)，即 \(V_2\in[-\sqrt{1-v^2},\sqrt{1-v^2}]\)。\(V_2\) 在 \([-1,1]\)（长度 2）上均匀，这个区间长度为 \(2\sqrt{1-v^2}\)，概率就是 \(\sqrt{1-v^2}\)。第二个等号：\(V_1\) 在 \([-1,1]\) 上均匀，\(|V_1|\) 在 \([0,1]\) 上均匀，而 \(\sqrt{1-V_1^2}\) 只依赖 \(|V_1|\)。 白话：原方法每次只记录「中 / 不中」（0 或 1），信息很粗。条件化之后，第二个坐标的随机性被精确积分掉，换成了一个光滑的数，噪声自然小。期权定价里用「给定波动率路径后的 Black–Scholes 公式」代替继续模拟股价，道理一样：能解析算的部分就别再随机模拟。

改进 2（再加对偶）：\(\sqrt{1-u^2}\) 在 \([0,1]\) 上单调递减，用 \(\frac12\big[\sqrt{1-U^2}+\sqrt{1-(1-U)^2}\big]\)。

原书 \(n=10000\) 的结果：落点比例 3.1612；\(\sqrt{1-U^2}\) 均值 3.128448；加对偶 3.139578（\(n=64000\) 时 3.143288）。下面的实战代码会量化这三种方法的方差差别。

在期权定价中，「条件蒙特卡洛」是同一思想：例如随机波动率模型中，先模拟波动率路径，再对给定波动率路径用 Black–Scholes 公式算出条件期望价格，而不是继续模拟股价。

10.4.4 控制变量

若已知某个函数的期望 \(E[f(\mathbf X)]=\mu\)，则对任意常数 \(a\)，

\(a^*\) 就是 \(g\) 对 \(f\) 回归斜率的相反数，与第 07b 章的最优线性预测、与最小方差对冲比率完全同构：控制变量就是「用一个已知期望的资产对冲掉估计中的噪声」。

推导拆解：\(\mathrm{Var}(W)\) 来自和的方差公式（第 07a 章式 4.1），\(\mu\) 是常数，不影响方差。它是 \(a\) 的二次函数，求导 \(2a\mathrm{Var}[f]+2\mathrm{Cov}[g,f]=0\) 得 \(a^*\)。代回：\(\mathrm{Var}[g]+\frac{\mathrm{Cov}^2}{\mathrm{Var}[f]}-2\frac{\mathrm{Cov}^2}{\mathrm{Var}[f]}=\mathrm{Var}[g]-\frac{\mathrm{Cov}^2}{\mathrm{Var}[f]}=\mathrm{Var}[g](1-\rho^2)\)。 金融直觉：对照 CFA 里的最小方差对冲比率 \(h^*=\rho\frac{\sigma_S}{\sigma_F}=\frac{\mathrm{Cov}(S,F)}{\mathrm{Var}(F)}\)：持有现货 \(g\)、卖出 \(h^*\) 份期货 \(f\)，对冲后的方差是 \(\sigma_S^2(1-\rho^2)\)。控制变量做的是同一件事，只是对冲的对象是「模拟误差」。\(f\) 的期望 \(\mu\) 已知，所以加上 \(a(f-\mu)\) 不改变期望，就像对冲工具按公允价格入场，不改变组合的期望价值。相关性越高，对冲越干净：亚式期权例子里 \(\rho^2=0.999\)，方差降了一千多倍。

原书自测题 10.5 及其两处问题 题目：\(X,Y\) 为独立的均值 1 的指数随机变量，用模拟估计 \(E[e^{XY}]\)，并用控制变量改进。原书解答给出两种控制变量：\(XY\)（期望 1），以及 Taylor 展开的前两项 \(XY+X^2Y^2/2\)。

• 常数印刷错误：原书写作 \(c\big[X_iY_i+X_i^2Y_i^2/2-1/2\big]\)，但
  \[E\Big[XY+\frac{X^2Y^2}{2}\Big]=E[X]E[Y]+\frac12E[X^2]E[Y^2]=1+\frac12\cdot2\cdot2=3,\]
  所以应减去 3 而不是 1/2，否则估计量有偏。（同理，用 \(XY\) 作控制变量时应写成 \(e^{XY}+c(XY-1)\)。）
• 题设本身的问题：在这个题设下 \(E[e^{XY}]\) 其实是无穷大。因为 \(E[e^{XY}\mid Y=y]=E[e^{yX}]=\frac{1}{1-y}\) 只在 \(y<1\) 时有限，\(y\ge1\) 时为无穷，而 \(P\{Y\ge1\}=e^{-1}>0\)。因此样本均值不会收敛，控制变量也无从改进。这个例子提醒我们：蒙特卡洛的前提是被估计的期望存在，最好方差也有限；对重尾的被积函数，模拟结果会表现为「样本量越大，估计越大、越不稳定」。若把题设改为 \(X,Y\) 独立 \(U(0,1)\)，\(E[e^{XY}]\) 有限，控制变量 \(XY+X^2Y^2/2\) 的期望变为 \(\frac14+\frac1{18}=\frac{11}{36}\)，方法完全适用。

10.4.5 重要性抽样（习题 10.16）

要估计 \(\theta=\int_0^1g(x)\,dx\)，可以从任一密度 \(f\) 抽 \(X\)，用 \(\frac{g(X)}{f(X)}\) 的均值估计，因为 \(E_f\big[\frac{g(X)}{f(X)}\big]=\int g\)（\(E_f\) 表示 \(X\) 按密度 \(f\) 抽样时的期望；按定义 \(E_f[\frac{g}{f}]=\int\frac{g(x)}{f(x)}f(x)\,dx\)，\(f\) 约掉即得）。当 \(f\) 的形状与 \(g\) 接近时，比值 \(g/f\) 接近常数，方差很小。这就是重要性抽样（importance sampling）：把样本集中到「对积分贡献大」的区域。估计深度虚值期权价格、极端损失概率（VaR 尾部）时，大部分普通样本对结果毫无贡献，重要性抽样（如把正态的均值平移到行权价附近，再按似然比加权）可以把方差降低几个数量级。

金融直觉：估计 \(P\{Z>4\}\approx3.2\times10^{-5}\) 这类尾部概率时，直接模拟 10 万次平均只有约 3 次落进尾部，估计几乎全是噪声。重要性抽样相当于「故意多抽危机情景，再按真实概率给它们降权」：从 \(N(4,1)\) 抽样，大约一半样本落在尾部，每个样本乘上似然比 \(\frac{\varphi(x)}{\varphi(x-4)}=e^{-4x+8}\)（真实密度 / 抽样密度）纠正偏差（练习 9）。压力测试里对极端情景加大采样、再按发生概率加权汇总，是同一个思路。 风险：如果抽样密度 \(f\) 在 \(g\) 不为 0 的某些区域太小，权重 \(g/f\) 会偶尔极大，方差反而爆炸。选 \(f\) 时要保证它的尾部不比 \(g\) 更薄。

量化实战

1. 生成器的正确性检查。 下面先实现本章的各种生成算法并核对理论值，再重现原书估计 \(\pi\) 的三种方法。

import numpy as np
from scipy import stats
rng = np.random.default_rng(10)
n = 200_000

# ---------- 1. 逆变换：指数、Weibull ----------
expo = -np.log(rng.random(n)) / 2.0                     # Exp(λ=2)
weib = (-np.log(rng.random(n)) / 0.5) ** (1 / 1.5)       # F(t)=1-exp(-a t^β), a=0.5, β=1.5
wb = stats.weibull_min(1.5, scale=0.5 ** (-1 / 1.5))
print(f"Exp(2):  均值 {expo.mean():.4f} (理论 0.5000), 中位数 {np.median(expo):.4f} (理论 {np.log(2)/2:.4f})")
print(f"Weibull: 均值 {weib.mean():.4f} (理论 {wb.mean():.4f}), 中位数 {np.median(weib):.4f} (理论 {wb.median():.4f})")

def describe(z):
    return f"均值 {z.mean():+.4f}, 标准差 {z.std():.4f}, P(|Z|>1.96) {np.mean(np.abs(z) > 1.96):.4f}"

# ---------- 2. 拒绝法：用 Exp(1) 生成 |Z|，再随机取号 ----------
def normal_by_rejection(m):
    out, tries = [], 0
    while len(out) < m:
        y1 = rng.exponential(1.0, m); y2 = rng.exponential(1.0, m)
        tries += m
        out.extend(y1[y2 >= (y1 - 1) ** 2 / 2])
    z = np.array(out[:m]); sign = np.where(rng.random(m) < 0.5, -1, 1)
    return z * sign, tries / len(out)
z_rej, iters = normal_by_rejection(n)
print(f"拒绝法正态: 每个接受样本平均尝试 {iters:.3f} 次 (理论 c = {np.sqrt(2*np.e/np.pi):.3f})\n  " + describe(z_rej))

# ---------- 3. Marsaglia 极坐标法 ----------
V = 2 * rng.random((n, 2)) - 1
S = (V ** 2).sum(1)
keep = (S <= 1) & (S > 0)
fac = np.sqrt(-2 * np.log(S[keep]) / S[keep])
z_polar = np.concatenate([fac * V[keep, 0], fac * V[keep, 1]])
print(f"极坐标法: 接受率 {keep.mean():.4f} (理论 π/4 = {np.pi/4:.4f})\n  " + describe(z_polar))

# ---------- 4. Poisson：均匀数连乘直到 < e^{-λ} ----------
def poisson_by_product(lam):
    k, prod, thresh = 0, rng.random(), np.exp(-lam)
    while prod >= thresh:
        prod *= rng.random(); k += 1
    return k
pois = np.array([poisson_by_product(4.0) for _ in range(50_000)])
print(f"Poisson(4): 均值 {pois.mean():.3f}, 方差 {pois.var():.3f}")

# ---------- 5. 估计 π：原始 / 条件化 / 条件化 + 对偶（原书 10.4 节） ----------
m = 10_000
V = 2 * rng.random((m, 2)) - 1
est_hit = 4 * ((V ** 2).sum(1) <= 1)                     # 每个估计用 2 个随机数
U = rng.random(m)
est_cond = 4 * np.sqrt(1 - U ** 2)                       # 每个用 1 个随机数
U2 = rng.random(m // 2)
est_anti = 2 * (np.sqrt(1 - U2 ** 2) + np.sqrt(1 - (1 - U2) ** 2))   # 每个用 1 个随机数
for name, e, per in [("落点比例", est_hit, 2), ("条件化 √(1-U²)", est_cond, 1),
                     ("条件化+对偶", est_anti, 1)]:
    print(f"{name:<14s} π估计 {e.mean():.5f}  每个随机数的方差 {e.var() * per:.4f}")

关键输出：

Exp(2):  均值 0.5002 (理论 0.5000), 中位数 0.3480 (理论 0.3466)
Weibull: 均值 1.4276 (理论 1.4330), 中位数 1.2376 (理论 1.2433)
拒绝法正态: 每个接受样本平均尝试 1.316 次 (理论 c = 1.315)
  均值 -0.0027, 标准差 1.0001, P(|Z|>1.96) 0.0498
极坐标法: 接受率 0.7857 (理论 π/4 = 0.7854)
  均值 -0.0016, 标准差 1.0003, P(|Z|>1.96) 0.0503
Poisson(4): 均值 3.997, 方差 4.000
落点比例           π估计 3.13320  每个随机数的方差 5.4317
条件化 √(1-U²)    π估计 3.13086  每个随机数的方差 0.8098
条件化+对偶         π估计 3.14145  每个随机数的方差 0.1105

「每个随机数的方差」把随机数消耗统一起来比较：同样的随机数预算下，条件化把方差降为原来的约 1/7，再加对偶又降为约 1/7，合计约 50 倍——相当于把样本量放大 50 倍。

2. 期权定价中的方差缩减。 欧式看涨期权有 Black–Scholes 解析价，可以用来检验各种方法。控制变量选折现后的标的价格 \(e^{-rT}S_T\)，在风险中性测度下其期望就是 \(S_0\)（第 07c 章的对数正态均值公式）。算术平均亚式期权没有解析解，但几何平均亚式期权有（几何平均的对数是正态的），两者高度相关，是教科书级的控制变量。最后核对自测题 10.5 的常数。

import numpy as np
from scipy import stats
rng = np.random.default_rng(100)
S0, K, r, sig, T = 100.0, 105.0, 0.03, 0.25, 0.5
n = 100_000
disc = np.exp(-r * T)

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

def report(name, samples, n_normals):
    est, se = samples.mean(), samples.std(ddof=1) / np.sqrt(samples.size)
    print(f"{name:<18s} 估计 {est:.4f}  标准误 {se:.4f}  (用 {n_normals:,} 个正态数)")
    return se

# ---------- 1. 欧式看涨：普通 / 对偶 / 控制变量 ----------
print(f"Black-Scholes 解析价 {bs_call(S0, K, r, sig, T):.4f}")
Z = rng.standard_normal(n)
ST = lambda z: S0 * np.exp((r - sig**2 / 2) * T + sig * np.sqrt(T) * z)
pay = disc * np.maximum(ST(Z) - K, 0)
se0 = report("普通蒙特卡洛", pay, n)

Zh = Z[: n // 2]                                            # 同样用 n 个「路径」的成本
pay_anti = 0.5 * disc * (np.maximum(ST(Zh) - K, 0) + np.maximum(ST(-Zh) - K, 0))
se1 = report("对偶变量 Z,-Z", pay_anti, n // 2)

ctrl = disc * ST(Z)                                         # 控制变量：折现标的价，E = S0
c = -np.cov(pay, ctrl)[0, 1] / ctrl.var(ddof=1)             # a* = -Cov/Var
pay_cv = pay + c * (ctrl - S0)
se2 = report("控制变量 S_T", pay_cv, n)
print(f"  估计的 a* = {c:.4f}, Corr(收益, 控制)^2 = {np.corrcoef(pay, ctrl)[0,1]**2:.3f}")
print(f"  方差缩减倍数: 对偶 {(se0/se1)**2:.2f}x, 控制变量 {(se0/se2)**2:.2f}x")

# ---------- 2. 算术平均亚式期权：用几何平均亚式（有解析解）作控制变量 ----------
m_steps, n_paths = 50, 50_000
dt = T / m_steps
Zp = rng.standard_normal((n_paths, m_steps))
logS = np.log(S0) + np.cumsum((r - sig**2 / 2) * dt + sig * np.sqrt(dt) * Zp, axis=1)
arith = disc * np.maximum(np.exp(logS).mean(1) - K, 0)
geo = disc * np.maximum(np.exp(logS.mean(1)) - K, 0)
# 几何平均的对数服从正态：均值与方差可显式计算
t_i = dt * np.arange(1, m_steps + 1)
mu_g = np.log(S0) + (r - sig**2 / 2) * t_i.mean()
var_g = sig**2 * np.sum(np.minimum.outer(t_i, t_i)) / m_steps**2
d1 = (mu_g - np.log(K) + var_g) / np.sqrt(var_g)
geo_exact = disc * (np.exp(mu_g + var_g / 2) * stats.norm.cdf(d1) - K * stats.norm.cdf(d1 - np.sqrt(var_g)))
b = -np.cov(arith, geo)[0, 1] / geo.var(ddof=1)
arith_cv = arith + b * (geo - geo_exact)
print(f"\n几何亚式解析价 {geo_exact:.4f}, 模拟 {geo.mean():.4f}")
s_a = report("算术亚式 普通", arith, n_paths * m_steps)
s_b = report("算术亚式 控制变量", arith_cv, n_paths * m_steps)
print(f"  Corr^2 = {np.corrcoef(arith, geo)[0,1]**2:.4f}, 方差缩减 {(s_a/s_b)**2:.0f}x")

# ---------- 3. 原书自测题 10.5 的控制变量常数 ----------
X, Y = rng.exponential(1.0, (2, 2_000_000))
f = X * Y + (X * Y) ** 2 / 2
print(f"\nE[XY + X²Y²/2] 模拟 = {f.mean():.4f}（应为 1 + E[X²]E[Y²]/2 = 3）")
g = np.exp(np.minimum(X * Y, 700))
print("E[e^{XY}] 的样本均值随样本量: " +
      ", ".join(f"n={k:.0e}: {g[:k].mean():.3g}" for k in [10**4, 10**5, 10**6, 2 * 10**6]))

关键输出：

Black-Scholes 解析价 5.5760
普通蒙特卡洛             估计 5.6042  标准误 0.0327  (用 100,000 个正态数)
对偶变量 Z,-Z          估计 5.5699  标准误 0.0275  (用 50,000 个正态数)
控制变量 S_T           估计 5.5812  标准误 0.0170  (用 100,000 个正态数)
  估计的 a* = -0.4956, Corr(收益, 控制)^2 = 0.729
  方差缩减倍数: 对偶 1.42x, 控制变量 3.70x

几何亚式解析价 2.3348, 模拟 2.3443
算术亚式 普通            估计 2.4597  标准误 0.0228  (用 2,500,000 个正态数)
算术亚式 控制变量          估计 2.4498  标准误 0.0007  (用 2,500,000 个正态数)
  Corr^2 = 0.9991, 方差缩减 1172x

E[XY + X²Y²/2] 模拟 = 3.0021（应为 1 + E[X²]E[Y²]/2 = 3）
E[e^{XY}] 的样本均值随样本量: n=1e+04: 1.88e+06, n=1e+05: 3.16e+15, n=1e+06: 3.78e+22, n=2e+06: 1.89e+22

解读：

• 对偶变量在同样的收益函数计算次数下把方差降到 1/1.42，同时只用了一半的随机数；控制变量的缩减倍数 3.70 与理论值 \(1/(1-\rho^2)=1/(1-0.729)\approx3.69\) 吻合。看涨期权收益只在 \(S_T>K\) 时与 \(S_T\) 线性相关，所以相关系数不算很高；越是实值的期权，控制变量越有效。
• 算术亚式期权与几何亚式期权的相关系数平方高达 0.999，控制变量把方差降了一千多倍，相当于把模拟路径数扩大一千倍。几何亚式的模拟值 2.3443 与解析价 2.3348 相差 0.4 个标准误左右，说明解析式正确。
• 自测题 10.5：\(E[XY+X^2Y^2/2]\) 的模拟值为 3.0021，确认常数应为 3。\(E[e^{XY}]\) 的样本均值随样本量暴涨到 \(10^{22}\) 量级（代码里把指数截断在 \(e^{700}\) 以防溢出），说明这个期望是无穷大，题设下的模拟没有意义。

3. 随机排列与置换检验。 用 Fisher–Yates 算法打乱收益，构造「因子与收益无关」的零分布，检验截面 IC 的显著性。

import numpy as np
rng = np.random.default_rng(11)

def fisher_yates(x):
    """原书例 1a：从最后一个位置开始，与 1..I 中均匀选出的位置交换"""
    x = x.copy()
    for I in range(len(x) - 1, 0, -1):
        N = int((I + 1) * rng.random())        # 0..I 上均匀
        x[N], x[I] = x[I], x[N]
    return x

# 均匀性检查：4 个元素的 24 种排列出现频率
counts = {}
for _ in range(240_000):
    k = tuple(fisher_yates(np.arange(4)))
    counts[k] = counts.get(k, 0) + 1
freq = np.array(list(counts.values())) / 240_000
print(f"排列种数 {len(counts)}, 频率范围 [{freq.min():.4f}, {freq.max():.4f}] (理论 {1/24:.4f})")

# 置换检验：一个截面因子的 IC 是否显著
n = 500
factor = rng.standard_normal(n)
ret = 0.08 * factor + rng.standard_normal(n)
ic = np.corrcoef(factor, ret)[0, 1]
null = np.array([np.corrcoef(factor, fisher_yates(ret))[0, 1] for _ in range(5000)])
print(f"观测 IC = {ic:.4f}, 置换分布标准差 = {null.std():.4f} (≈1/√n = {1/np.sqrt(n):.4f}), "
      f"双侧 p 值 = {(np.abs(null) >= abs(ic)).mean():.4f}")

关键输出：

排列种数 24, 频率范围 [0.0410, 0.0428] (理论 0.0417)
观测 IC = 0.0637, 置换分布标准差 = 0.0445 (≈1/√n = 0.0447), 双侧 p 值 = 0.1512

解读：24 种排列的频率都在 1/24 附近；500 只股票的单期截面上，IC 的零分布标准差约为 \(1/\sqrt{500}\approx0.045\)，所以 0.064 的单期 IC 远不显著。真实 IC 约 0.08 的因子要靠多期累积才能确认（第 08 章的样本量逻辑）。

本章小结

模拟以伪随机数为起点，强大数定律保证模拟均值收敛，CLT 给出误差 \(\sigma/\sqrt k\)。逆变换法 \(F^{-1}(U)\) 适用于分布函数可显式求逆的情形；拒绝法借助易抽样的 \(g\)，以 \(f/(cg)\) 的概率接受，平均迭代 \(c\) 次；正态分布可用指数拒绝法、Box–Muller 或 Marsaglia 极坐标法；Gamma、卡方、Poisson、几何、二项都有基于均匀数的简单构造。由于误差只按 \(1/\sqrt k\) 下降，方差缩减往往比增加样本更划算：对偶变量要求被积函数单调；条件化利用全方差公式；控制变量的最优系数是 \(-\mathrm{Cov}/\mathrm{Var}\)，缩减比例为 \(\rho^2\)；重要性抽样把样本集中到重要区域。所有这些都以被估计的期望存在为前提。

练习

基础

1. 写出用逆变换法从密度 \(f(x)=\frac{e^x}{e-1}\)（\(0<x<1\)）抽样的公式。 答案：\(F(x)=\frac{e^x-1}{e-1}\)，\(X=\log\big(U(e-1)+1\big)\)（自测题 10.1）。
2. 用拒绝法（\(g\) 为 \(U(0,1)\)）从密度 \(f(x)=30(x^2-2x^3+x^4)\) 抽样，求 \(c\) 与平均迭代次数。 答案：\(f\) 在 \(x=1/2\) 处最大，\(c=15/8\)，平均迭代 1.875 次；接受条件 \(U_2\le16(U_1^2-2U_1^3+U_1^4)\)（自测题 10.2）。
3. 离散分布 \(P\{X=1\}=0.15\)，\(P\{X=2\}=0.2\)，\(P\{X=3\}=0.35\)，\(P\{X=4\}=0.3\)。写出平均比较次数最少的逆变换算法。 答案：依次判断 \(U\le0.35\to3\)，\(\le0.65\to4\)，\(\le0.85\to2\)，否则 1（自测题 10.3）。
4. 用逆变换法生成违约强度为常数 \(h\) 的违约时间；若强度为 \(h(t)=ct\) 呢？ 答案：\(-\frac1h\log U\)；\(\sqrt{-2\log U/c}\)（习题 10.6）。
5. 说明为什么对跨式期权（收益 \(|S_T-K|\)）使用 \(Z,-Z\) 对偶变量效果可能很差。

进阶

6. 证明拒绝法中所需迭代次数服从几何分布，均值为 \(c\)；并说明 \(g\) 选为 Exp\((\lambda)\) 生成 \(|Z|\) 时 \(\lambda=1\) 最优（习题 10.10）。
7. 证明：若 \((V_1,V_2)\) 在单位圆盘上均匀，则 \(R^2=V_1^2+V_2^2\sim U(0,1)\)，且与角度 \(\Theta\) 独立（习题 10.13）。
8. 推导控制变量的最优系数 \(a^*\) 及最小方差（习题 10.15），并说明它与最小方差对冲比率的关系。
9. 用重要性抽样估计 \(P\{Z>4\}\)（\(Z\) 标准正态）：从 \(N(4,1)\) 抽样，写出估计量并说明为什么方差比直接模拟小得多。 提示：估计量为 \(I(X>4)\,e^{-4X+8}\)，\(X\sim N(4,1)\)；直接模拟时每 3 万个样本才约有 1 个落入区域。
10. 修正自测题 10.5：把 \(X,Y\) 改为独立 \(U(0,1)\)，写出以 \(XY+X^2Y^2/2\) 为控制变量的无偏估计量，并编程比较方差缩减倍数。 答案要点：\(e^{XY}+a\big(XY+\frac{X^2Y^2}{2}-\frac{11}{36}\big)\)，\(a\) 用样本 \(-\mathrm{Cov}/\mathrm{Var}\) 估计。

原书推荐习题：Problems 10.1（另一种随机排列算法的正确性）、10.5–10.6（Weibull 与失效率模拟）、10.9（组合法）、10.10（拒绝法参数优化）、10.12（蒙特卡洛积分）、10.13、10.15（控制变量最优系数）、10.16（重要性抽样，强烈推荐）；Self-Test 10.1–10.5（10.5 注意本章指出的两处问题）。进一步阅读：Ross《Simulation》第 5 版。

原书对照