Monte Carlo estimation

Answering with averages of random draws

Many questions in finance have no tidy formula: the chance that a portfolio loses 20% within a year, the price of an option whose payoff depends on the whole path. Monte Carlo answers them with brute force. Write the quantity as an average, E[f(X)], simulate many draws of X, apply f to each and take the mean. By the law of large numbers the sample mean converges to the true expectation. A probability is just the expectation of a yes/no indicator, so "the fraction of simulations where the event happened" estimates it.

import math
import numpy as np

# Estimate pi: the quarter circle x^2 + y^2 <= 1 covers pi/4 of the unit square
rng = np.random.default_rng(42)
x, y = rng.random(100000), rng.random(100000)
inside = x ** 2 + y ** 2 <= 1                        # True/False for every point
print(round(4 * inside.mean(), 4))                    # 3.1378   close to 3.1416, with 100,000 points

Simulating prices and probabilities

To simulate a stock, draw standard normal shocks and push them through the lognormal formula from the previous topic: all n terminal prices come out of one array expression. Then any question about the distribution of prices becomes an average over the array. The probability of finishing below 80 is np.mean(prices < 80), the expected price is prices.mean(), and a quantile of the losses is np.quantile(...). Pass the generator rng into the function instead of creating it inside, so the caller controls repeatability.

def simulate_terminal_prices(rng, s0, mu, sigma, t, n):
    z = rng.standard_normal(n)
    return s0 * np.exp((mu - 0.5 * sigma ** 2) * t + sigma * math.sqrt(t) * z)

prices = simulate_terminal_prices(np.random.default_rng(1), 100, 0.05, 0.2, 1, 200000)
print(round(float(np.mean(prices < 80)), 4))          # 0.1027   P(price below 80 after a year), about 10%
print(round(float(np.mean(prices > 130)), 4))         # 0.1216   P(price above 130)

Pricing an option

A European call pays max(S - K, 0) at expiry. In the risk-neutral world every asset grows at the risk-free rate, so the fair price is the discounted expected payoff: simulate terminal prices with mu = rate, compute the payoff of each, average, and multiply by exp(-rate * t). For this simple option a closed formula exists (Black-Scholes, 10.4506 for the numbers below), and the simulation lands close to it. The point of Monte Carlo is that the same few lines still work for payoffs with no formula at all.

def mc_call_price(rng, s0, strike, rate, sigma, t, n):
    prices = simulate_terminal_prices(rng, s0, rate, sigma, t, n)    # risk-neutral: the price grows at the risk-free rate
    payoff = np.maximum(prices - strike, 0)                           # max(S - K, 0) for every simulated path
    return math.exp(-rate * t) * payoff.mean()                        # discounted average payoff

price = mc_call_price(np.random.default_rng(1), 100, 100, 0.05, 0.2, 1, 200000)
print(round(price, 4))        # 10.3993   the exact Black-Scholes price is 10.4506: the difference is Monte Carlo error

# More points, less error (1 / sqrt(n)): the estimate of pi with the same seed
for n in (100, 10000, 1000000):
    g = np.random.default_rng(42)
    x, y = g.random(n), g.random(n)
    print(n, round(4 * float((x * x + y * y <= 1).mean()), 4))
# 100 3.36
# 10000 3.1488
# 1000000 3.1424

• Monte Carlo: write the quantity as an average E[f(X)], simulate many draws, apply f, take the mean.
• A probability is the mean of a comparison: np.mean(prices < 80); pass the generator in so the run is repeatable.
• Option price = exp(-rate * t) * mean(max(S - K, 0)) with the prices simulated at mu = rate; report the error too.

Next up: Monte Carlo error and convergence, how accurate an estimate is and how to get more for less.

Lab · what's graded

• Throw darts at a square to find pi with estimate_pi25 pts
• Count how often an event happens with estimate_probability25 pts
• Simulate where a stock could end up with simulate_terminal_prices25 pts
• Price an option by simulation with mc_call_price25 pts