控制变量法

利用已知信息

模拟看涨期权时，你同时也模拟了到期股价 S_T，它的真实均值 S0 * exp(r * T) 是精确已知的。样本里 S_T 的平均偏离了一点，比如偏高。由于看涨期权收益与 S_T 同向变动，看涨期权的样本均值很可能也偏高。控制变量正是利用这一点：从想要的量的估计中，减去已知量观测误差的若干倍。

import numpy as np

rng = np.random.default_rng(1)
z = rng.standard_normal(2000)
s_t = 100.0 * np.exp((0.05 - 0.5 * 0.2 ** 2) + 0.2 * z)             # terminal prices
call = np.maximum(s_t - 100.0, 0.0)                                  # the payoff we want

known = 100.0 * np.exp(0.05)                                         # E[S_T], exactly
print(round(float(s_t.mean()), 3), round(known, 3))                  # the sample mean of S_T vs its true mean
print(round(float(s_t.mean() - known), 3))                           # the observed error of the control

最优倍数

应该减去多少控制变量的误差？使 payoff - b * control 方差最小的倍数是 b = cov(payoff, control) / var(control)，恰好是收益对控制变量回归的斜率（与 CAPM 一节中贝塔用的是相同的各项和）。修正后的估计是 mean(payoff) - b * (mean(control) - known_mean)。它是无偏的，因为被减去的项期望为零。

import numpy as np

def control_coefficient(payoff, control):
    dev = control - control.mean()
    return float(np.sum(dev * (payoff - payoff.mean())) / np.sum(dev ** 2))

rng = np.random.default_rng(1)
z = rng.standard_normal(2000)
s_t = 100.0 * np.exp((0.05 - 0.5 * 0.2 ** 2) + 0.2 * z)
call = np.maximum(s_t - 100.0, 0.0)
known = 100.0 * np.exp(0.05)

b = control_coefficient(call, s_t)
controlled = call.mean() - b * (s_t.mean() - known)
print(round(b, 4), round(float(call.mean()), 4), round(float(controlled), 4))       # the plain mean and the corrected one

# A payoff that is exactly proportional to the control is estimated without any error
import numpy as np

x = np.random.default_rng(2).standard_normal(50) + 5.0          # the control, known mean 5.0
y = 2.0 * x                                                      # the payoff
dev = x - x.mean()
b = np.sum(dev * (y - y.mean())) / np.sum(dev ** 2)
print(round(float(b), 6), round(float(y.mean() - b * (x.mean() - 5.0)), 6))      # 2.0 and exactly 10.0

节省多少误差

携带噪声的是调整后的样本 payoff - b * control，所以修正后估计的标准误就是其标准差（ddof=1）除以 sqrt(n)。设收益与控制变量的相关系数为 rho，方差会降为原来的 1 - rho^2：相关 0.9 去掉 81%，0.99 去掉 98%。控制变量不需要额外的模拟，只多几行代码。

import numpy as np

def control_variate_error(payoff, control):
    dev = control - control.mean()
    b = np.sum(dev * (payoff - payoff.mean())) / np.sum(dev ** 2)
    return float((payoff - b * control).std(ddof=1) / np.sqrt(len(payoff)))

rng = np.random.default_rng(1)
z = rng.standard_normal(20000)
s_t = 100.0 * np.exp((0.05 - 0.5 * 0.2 ** 2) + 0.2 * z)
for strike in (80.0, 100.0, 120.0):
    call = np.maximum(s_t - strike, 0.0)
    plain = float(call.std(ddof=1) / np.sqrt(len(call)))
    print(strike, round(plain, 4), round(control_variate_error(call, s_t), 4), round(plain / control_variate_error(call, s_t), 2))

• 控制变量是均值已知的模拟量；从估计中减去其观测误差的若干倍：mean(Y) - b * (mean(X) - E[X])。
• 最优倍数是回归斜率 b = cov(Y, X) / var(X)；误差为 Y - b * X 的标准差除以 sqrt(n)。
• 方差降为 1 - rho^2 倍；好的控制变量与收益高度相关、均值精确已知，且不需要额外模拟。

下一节：用于收益率的 scikit-learn 流水线，在收益率数据上拟合并评估预测模型。

实验：评分内容

• 用 control_coefficient 选取合适的倍数25 分
• 用 control_variate_estimate 修正估计25 分
• 用 control_variate_error 缩小误差条25 分
• 用 error_reduction 计算改进倍数25 分