Antithetic variates

Accuracy is expensive

The standard error of a Monte Carlo estimate falls like 1 / sqrt(N): ten times more accuracy needs a hundred times more paths. When one simulation is costly, a better lever than more paths is less noise per path. These are variance reduction methods, and the simplest is antithetic variates. They use the same random numbers twice: with a draw z we also run the path for -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

Why pairing helps

If z is large and positive, the path ends high and a call pays a lot; with -z it ends low and the call pays little. Their payoffs are negatively correlated, so the pair's average is much steadier than either one. For a monotone payoff such as a call the pair average (f(z) + f(-z)) / 2 has a variance (1 + corr) / 2 times that of one payoff, with corr negative. The price is still the average of an unbiased quantity, so there is no bias.

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

The estimator and its error bar

Treat each pair as one observation: pair = (payoff(z) + payoff(-z)) / 2. The estimate is the discounted mean of the pair payoffs, and the standard error is the discounted standard deviation (ddof=1) of the pairs over sqrt(number of pairs). Note that the count is the number of *pairs*, since the two members of a pair are not independent. To judge the gain, compare with a plain estimate that uses the same number of payoffs from independent draws.

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)

• Use each draw z together with -z; average the two payoffs into one pair payoff, which is much less noisy for a monotone payoff.
• The standard error is computed over the number of pairs, with ddof=1: discount * std(pairs) / sqrt(pairs).
• Compare with a plain estimate of the same cost; no help for payoffs symmetric in z, and it combines with control variates.

Next up: control variates, using a quantity whose exact value you know to cancel noise.

Lab · what's graded

• Turn draws into prices with terminal_price25 pts
• Report the estimate and its error with mc_estimate25 pts
• Pair each draw with its opposite via antithetic_call_price25 pts
• Measure the gain with variance_ratio25 pts