Case study: the error of a Monte Carlo price

The price is a mean

The inputs are those of case studies 2.4.6, 2.6.5 and 5.6.5. The stock costs $S_0 = 100$ dollars, the strike is $K = 100$, the expiry $T = 0.5$ years. The interest rate is $r = 0.03$ and the volatility $\sigma = 0.2$, both annual decimals.

Case study 5.6.5 wrote the price at expiry as $S_T = S_0e^{m + vZ}$, with $Z$ standard normal, $m = (r - \frac{\sigma^2}{2})T$ and $v = \sigma\sqrt{T}$. Case study 2.6.5 priced the call as an integral over $z$ weighted by $\varphi(z)$. That integral is an expected value: the mean of the discounted payoff

So the call's price is $c = \mathbb{E}[Y]$. Draw $N$ independent copies of $Z$, compute each $Y_i$, and average. By the law of large numbers (lesson 5.7.2) the average settles on $c$.

The payoff is lopsided

Before measuring the error of the average, look at what is being averaged. The call pays nothing when $S_T \le K$, which has probability $1 - N(d_2) \approx 0.4859$. When it pays, the payoff can be large.

The same integrals give its spread. $\mathbb{E}[Y^2] \approx 130.2379$, so the standard deviation of one discounted payoff is

That is about 1.49 times the price. The skewness (lesson 5.3.3) is 1.83: a long right tail.

The seeded sample looks the same. Of its 10,000 payoffs, 0.4903 are zero, the median is $0.3174, the largest is $66.43 and the sample skewness is 1.78. One payoff is nothing like normal; the central limit theorem (lesson 5.7.3) is about their mean.

The standard error and a 95% interval

The average of $N$ independent payoffs has standard deviation $\frac{\sigma_Y}{\sqrt{N}}$ (lesson 5.7.2): $\frac{9.4683}{100} \approx 0.0947$ for $N = 10{,}000$. On a real problem $\sigma_Y$ is unknown, since a formula for it would usually give the price too.

So estimate it from the same draws, with the sample standard deviation $s$, dividing by $N - 1$ (lesson 5.5.2). By the central limit theorem (lesson 5.7.3), $\hat{c}_N$ is close to normal. So it lies within 1.96 standard errors of $c$ about 95% of the time.

For the seeded run, $s \approx 9.3272$, so $\text{SE} \approx 0.0933$. The interval runs from 6.0885 to 6.4541 and covers 6.3710. The miss, $-0.0998$, is 1.07 standard errors below the price: ordinary, as the hook said.

What one more digit costs

Carry the same seeded run on to a million draws. The first 10,000 are the ones above.

Every interval covers 6.3710. Each tenfold step divides the standard error by about $\sqrt{10} \approx 3.16$ (lesson 5.7.3); case study 4.7.5 saw the same slope for a probability.

To pin the price to the cent with 95% confidence, a half-width of 0.01, needs $N = \left(\frac{1.96 \times 9.4683}{0.01}\right)^2$: about 3.44 million draws. A tenth of a cent needs about 344 million. Each extra digit costs a hundred times the draws.

When the 95% can be trusted

One run covering the price proves little. Repeat the whole simulation many times and count how often the interval covers 6.3710.

With 2,000 draws the coverage is close to 95%, and the misses fall about evenly on each side. With 20 draws it is only 90%, and most misses are low.

The lopsided payoff is the reason. The mean of 20 payoffs still has skewness about $\frac{1.83}{\sqrt{20}} \approx 0.41$ (lesson 5.7.3).

A small sample often misses the rare large payoffs. Then both its average and its $s$ come out too small: the interval sits low and is too narrow.

The interval needs enough draws for the theorem to apply, and more when the payoff is lopsided.

Antithetic draws

The standard error is $\frac{\sigma_Y}{\sqrt{N}}$, so there are two ways to shrink it: more draws, or a smaller $\sigma_Y$. Methods that take the second way are called variance reduction. The simplest uses antithetic draws.

If $Z$ is standard normal, so is $-Z$. Draw $\frac{N}{2}$ values of $Z$ and use each twice, as $Z$ and as $-Z$.

A high $Z$ gives a large payoff and its mirror a small one, so the two payoffs of a pair are negatively correlated. Their average varies less than two independent payoffs would.

With 10,000 payoffs in 5,000 pairs, one seeded run gives $6.3704 with standard error 0.0694. Plain draws give about 0.0947 for the same 10,000 payoffs. The correlation within a pair was $-0.4569$; the variance of a pair's average is $\frac{\sigma_Y^2(1 + \rho)}{2}$, with $\rho$ that correlation.

Two traps come with the method.

• The standard error must come from the pairs. The 10,000 payoffs are not independent, so $\frac{s}{\sqrt{10{,}000}}$ over all of them, 0.0942 here, is wrong. Use the 5,000 pair averages: their sample standard deviation over $\sqrt{5{,}000}$.
• It helps only when the payoff moves one way with $Z$. A straddle pays $|S_T - K|$, high at both ends. Its two payoffs in a pair have correlation $+0.8915$: near duplicates. An estimate from $\frac{N}{2}$ pairs then has $1 + \rho \approx 1.89$ times the variance of one from $N$ plain draws, so antithetic draws make it worse.

On the desk

import math

BS = 6.371028  # Black–Scholes price of the shared call (case study 2.4.6)


def test_seeded_run():
    """Seed 575 with 10,000 draws gives 6.2713 with standard error 0.0933, as in the case study"""
    price, se = mc_call(100, 100, 0.5, 0.03, 0.2, 10_000, 575)
    expect_close(price, 6.271264328239705, what="price", rel=1e-9)
    expect_close(se, 0.09327231109608519, what="standard error (sample sd with n - 1, over sqrt(n))", rel=1e-9)


def test_closes_on_black_scholes():
    """With 400,000 draws the estimate is within 3 standard errors of 6.3710"""
    price, se = mc_call(100, 100, 0.5, 0.03, 0.2, 400_000, 11)
    assert abs(price - BS) <= 3 * se, f"expected a price within 3 standard errors ({3 * se:.4f}) of {BS}, got {price:.4f}"


def test_standard_error_size():
    """The standard error is close to 9.4683 / sqrt(n), the payoff's standard deviation over sqrt(n)"""
    _, se = mc_call(100, 100, 0.5, 0.03, 0.2, 400_000, 12)
    expect_close(se, 9.4683 / math.sqrt(400_000), what="standard error with 400,000 draws", rel=0.02)


def test_four_times_the_draws():
    """Four times the draws halve the standard error"""
    small = mc_call(100, 100, 0.5, 0.03, 0.2, 10_000, 3)[1]
    large = mc_call(100, 100, 0.5, 0.03, 0.2, 40_000, 3)[1]
    expect_close(large / small, 0.5, what="ratio of standard errors for 40,000 and 10,000 draws", rel=0.03)


def test_antithetic():
    """Antithetic pairs, seed 5753 and 10,000 payoffs: 6.3704 with standard error 0.0694, from the 5,000 pair averages"""
    price, se = mc_call(100, 100, 0.5, 0.03, 0.2, 10_000, 5753, antithetic=True)
    expect_close(price, 6.370379263137437, what="antithetic price", rel=1e-9)
    expect_close(se, 0.0694189616918534, what="antithetic standard error", rel=1e-9)
    plain = mc_call(100, 100, 0.5, 0.03, 0.2, 10_000, 5753)[1]
    assert se < plain, f"expected the antithetic standard error {se:.4f} to be below the plain one {plain:.4f}"


def test_no_volatility():
    """With sigma = 0 every payoff is 100 - 100e^(-0.015) = 1.4888, and the standard error is 0"""
    price, se = mc_call(100, 100, 0.5, 0.03, 0.0, 1_000, 1)
    expect_close(price, 100 - 100 * math.exp(-0.015), what="price with sigma = 0", rel=1e-9)
    assert abs(se) < 1e-9, f"expected a standard error of 0, got {se}"


def test_at_expiry():
    """With T = 0 the call is worth max(S0 - K, 0): 5 for S0 = 105 and K = 100"""
    price, se = mc_call(105, 100, 0, 0.03, 0.2, 1_000, 1)
    expect_close(price, 5.0, what="price at expiry", rel=1e-12)
    assert abs(se) < 1e-9, f"expected a standard error of 0, got {se}"


def test_deep_out_of_the_money():
    """Strike 150: the standard error is more than 30% of the price with 10,000 draws"""
    price, se = mc_call(100, 150, 0.5, 0.03, 0.2, 10_000, 5754)
    assert price > 0, f"expected a positive price, got {price}"
    assert se / price > 0.3, f"expected a relative error above 0.3, got {se / price:.3f}"


def test_bad_inputs():
    """S0 <= 0, a negative sigma or T, n below 2, or an odd n with antithetic draws raises ValueError"""
    cases = [
        ("S0 = 0", lambda: mc_call(0, 100, 0.5, 0.03, 0.2, 100, 1)),
        ("sigma = -0.2", lambda: mc_call(100, 100, 0.5, 0.03, -0.2, 100, 1)),
        ("T = -1", lambda: mc_call(100, 100, -1, 0.03, 0.2, 100, 1)),
        ("n = 1", lambda: mc_call(100, 100, 0.5, 0.03, 0.2, 1, 1)),
        ("odd n with antithetic", lambda: mc_call(100, 100, 0.5, 0.03, 0.2, 101, 1, antithetic=True)),
    ]
    for name, call in cases:
        try:
            call()
        except ValueError:
            continue
        raise AssertionError(f"expected ValueError for {name}")