对偶变量法

精度代价高昂

蒙特卡洛估计的标准误按 1 / sqrt(N) 下降：精度提高 10 倍需要 100 倍的路径。当一次模拟代价很高时，比增加路径更好的杠杆是降低每条路径的噪声。这就是方差缩减方法，最简单的是对偶变量。它把同一批随机数用两次：对每个随机数 z，也跑 -z 的路径。

import numpy as np

def terminal_price(s0, rate, sigma, t, z):
    return s0 * np.exp((rate - 0.5 * sigma ** 2) * t + sigma * np.sqrt(t) * z)

z = np.array([-1.0, 0.0, 1.0])
print(np.round(terminal_price(100.0, 0.05, 0.2, 1.0, z), 3))         # the price for z = -1, 0, 1
print(np.round(terminal_price(100.0, 0.05, 0.2, 1.0, -z), 3))        # and for the opposite draws: mirror images

为什么配对有效

如果 z 很大且为正，路径收于高位，看涨期权赔付很多；用 -z 则收于低位，赔付很少。两个收益负相关，所以一对的平均值比单个稳定得多。对看涨期权这类单调收益，一对平均值 (f(z) + f(-z)) / 2 的方差是单个收益方差的 (1 + corr) / 2 倍，其中 corr 为负。价格仍是无偏量的平均，所以不会引入偏差。

import numpy as np

def terminal_price(s0, rate, sigma, t, z):
    return s0 * np.exp((rate - 0.5 * sigma ** 2) * t + sigma * np.sqrt(t) * z)

z = np.random.default_rng(1).standard_normal(20000)
up = np.maximum(terminal_price(100.0, 0.05, 0.2, 1.0, z) - 100.0, 0.0)       # call payoffs for z
down = np.maximum(terminal_price(100.0, 0.05, 0.2, 1.0, -z) - 100.0, 0.0)    # ... and for -z
print(round(float(np.corrcoef(up, down)[0, 1]), 3))                           # negative: the payoffs offset each other

估计量与误差条

把每一对当作一个观测：pair = (payoff(z) + payoff(-z)) / 2。估计值是各对收益的折现均值，标准误是各对收益的折现标准差（ddof=1）除以 sqrt(对数)。计数用的是*对*的个数，因为一对中的两个成员并不独立。要评判改进，就与使用同样数量收益、但来自独立随机数的普通估计比较。

import numpy as np

def terminal_price(s0, rate, sigma, t, z):
    return s0 * np.exp((rate - 0.5 * sigma ** 2) * t + sigma * np.sqrt(t) * z)

def mc_estimate(payoffs, discount):
    return discount * payoffs.mean(), discount * payoffs.std(ddof=1) / np.sqrt(len(payoffs))

def call_payoff(z):
    return np.maximum(terminal_price(100.0, 0.05, 0.2, 1.0, z) - 100.0, 0.0)

rng = np.random.default_rng(3)
z, w = rng.standard_normal(5000), rng.standard_normal(5000)
discount = np.exp(-0.05)

plain = mc_estimate(call_payoff(np.concatenate([z, w])), discount)           # 10,000 independent payoffs
anti = mc_estimate((call_payoff(z) + call_payoff(-z)) / 2, discount)         # 5,000 pairs: the same 10,000 evaluations
print([round(float(v), 4) for v in plain], [round(float(v), 4) for v in anti])
print(round(float(plain[1] / anti[1]), 2), "times smaller error")

# A call profits from pairing; a straddle (symmetric in z) does not
import numpy as np

def terminal_price(z):
    return 100.0 * np.exp((0.05 - 0.5 * 0.2 ** 2) + 0.2 * z)

rng = np.random.default_rng(4)
z = rng.standard_normal(20000)
for name, payoff in (("call", lambda p: np.maximum(p - 100.0, 0.0)), ("straddle", lambda p: np.abs(p - 100.0))):
    pairs = (payoff(terminal_price(z)) + payoff(terminal_price(-z))) / 2
    plain = payoff(terminal_price(np.concatenate([z, rng.standard_normal(20000)])))
    ratio = (plain.std(ddof=1) / np.sqrt(len(plain))) / (pairs.std(ddof=1) / np.sqrt(len(pairs)))
    print(name, round(float(ratio), 2))                       # the call gains (1.41); the straddle gets worse (0.78)

• 每个 z 与 -z 一起使用；把两个收益平均成一个配对收益，对单调收益噪声小得多。
• 标准误按对数计算，用 ddof=1：discount * std(pairs) / sqrt(pairs)。
• 与同成本的普通估计比较；对关于 z 对称的收益无益，并可与控制变量结合。

下一节：控制变量，利用精确值已知的量来抵消噪声。

实验：评分内容

• 用 terminal_price 把随机数变成价格25 分
• 用 mc_estimate 报告估计值及误差25 分
• 用 antithetic_call_price 把每个随机数与其相反数配对25 分
• 用 variance_ratio 度量改进25 分