精度は高くつく
モンテカルロ推定の標準誤差は 1 / sqrt(N) で減るため、精度を10倍にするには経路が100倍必要です。1回のシミュレーションが高価なら、経路を増やすよりも、経路あたりのノイズを減らすほうが有効です。これが分散削減の手法で、最も単純なのが対称変量法です。同じ乱数を2回使い、乱数 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 では低く終わってコールはほとんど払いません。2つのペイオフは負の相関を持つため、ペアの平均はどちらよりずっと安定します。コールのような単調なペイオフでは、ペア平均 (f(z) + f(-z)) / 2 の分散は1つのペイオフの (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推定量とその誤差
各ペアを1つの観測として扱います：pair = (payoff(z) + payoff(-z)) / 2。推定値はペアのペイオフの割引平均で、標準誤差はペアの割引標準偏差（ddof=1）を sqrt(ペア数) で割ったものです。ペアの2要素は独立でないので、数えるのは*ペア*の数です。改善幅を見るには、独立な乱数から同数のペイオフを使った通常の推定と比べます。
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 に対して対称なペイオフでは、f(z) と f(-z) が似ていてペアが情報を足さず、効果がないか逆効果になりえます。それでも追加の乱数は要らず無料で、次のトピックのコントロール変量とも組み合わせられます。- 各乱数
zを-zと一緒に使い、2つのペイオフを平均して1つのペア・ペイオフにする。単調なペイオフではノイズがずっと小さい - 標準誤差はペア数で計算する（
ddof=1）：discount * std(pairs) / sqrt(pairs) - 同じコストの通常の推定と比べる。
zに対称なペイオフには効かず、コントロール変量と併用できる
次はコントロール変量法です。正確な値が分かっている量を使ってノイズを打ち消します。
ラボ：採点される内容
terminal_priceで乱数を価格に変える25点mc_estimateで推定値と誤差を報告する25点antithetic_call_priceで乱数をその反対と組にする25点variance_ratioで改善幅を測る25点