Monte Carlo Simulation for Option Pricing

Why Monte Carlo?

Black-Scholes gives you an elegant closed-form formula for European options. But most real-world derivatives are more complex: Asian options depend on the average price, barrier options knock in or out at certain levels, lookback options depend on the maximum or minimum price. None of these have simple analytical solutions.

Monte Carlo simulation is the Swiss Army knife of quant finance. The idea is brutally simple: simulate thousands of possible stock price paths, calculate the payoff for each path, and take the average. By the law of large numbers, this average converges to the true price.

Every major investment bank and hedge fund in London uses Monte Carlo methods. At firms like Barclays, HSBC, and JP Morgan, MC engines price trillions of pounds worth of exotic derivatives daily.

Simulating Stock Paths

We use the exact solution to the GBM stochastic differential equation. This avoids discretisation error that comes from Euler schemes. Each path represents one possible future for the stock price.

import numpy as np

def simulate_gbm(S0, mu, sigma, T, dt, n_paths):
    """
    Simulate stock paths using Geometric Brownian Motion.
    dS = mu*S*dt + sigma*S*dW
    
    S0: initial stock price
    mu: expected return (drift)
    sigma: volatility
    T: time horizon (years)
    dt: time step (e.g. 1/252 for daily)
    n_paths: number of simulation paths
    """
    n_steps = int(T / dt)
    
    # Generate random normal shocks
    Z = np.random.standard_normal((n_steps, n_paths))
    
    # Pre-allocate price array
    S = np.zeros((n_steps + 1, n_paths))
    S[0] = S0
    
    # Simulate each time step
    for t in range(1, n_steps + 1):
        S[t] = S[t-1] * np.exp(
            (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z[t-1]
        )
    
    return S

With 10,000 paths and daily time steps, you get a rich picture of possible price trajectories. Some paths surge upward, some crash, most meander. The distribution of terminal prices will be log-normal, as Black-Scholes assumes.

Pricing a European Call via MC

For European options, we only need the terminal stock price (not the full path). This makes it extremely efficient. The key insight: we simulate under the risk-neutral measure, using the risk-free rate r as the drift instead of the real-world expected return μ.

def mc_european_call(S0, K, T, r, sigma, n_paths=100000):
    """Price a European call option via Monte Carlo."""
    dt = T  # Single step for European (only need terminal price)
    
    # Simulate terminal prices under risk-neutral measure
    Z = np.random.standard_normal(n_paths)
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
    
    # Calculate payoffs
    payoffs = np.maximum(ST - K, 0)
    
    # Discount to present value
    price = np.exp(-r * T) * np.mean(payoffs)
    
    # Standard error
    se = np.exp(-r * T) * np.std(payoffs) / np.sqrt(n_paths)
    
    return price, se

# Price a call: S=100, K=105, T=0.5, r=5%, sigma=20%
price, se = mc_european_call(100, 105, 0.5, 0.05, 0.2, n_paths=100000)
print(f"MC Call Price: {price:.4f} (+/- {1.96*se:.4f})")
print(f"BS Call Price: 3.7724")  # From Black-Scholes for comparison

Convergence to Black-Scholes

As you increase the number of simulated paths, the MC price converges to the Black-Scholes analytical price. The error decreases proportionally to 1/√N, meaning you need 4x more paths to halve the error. This is a fundamental limitation of naive Monte Carlo.

def convergence_analysis(S0, K, T, r, sigma, bs_price):
    """Show how MC converges to BS price as paths increase."""
    path_counts = [100, 500, 1000, 5000, 10000, 50000, 100000, 500000]
    
    print(f"Paths      | MC Price   | Error      | SE")
    print("-" * 50)
    
    for n in path_counts:
        price, se = mc_european_call(S0, K, T, r, sigma, n)
        error = price - bs_price
        print(f"{n:>10} | {price:.4f}    | {error:.4f}    | {se:.4f}")

convergence_analysis(100, 105, 0.5, 0.05, 0.2, 3.7724)

Pricing an Exotic: Asian Option

This is where Monte Carlo truly shines. An Asian option payoff depends on the average stock price over the life of the option, not just the terminal price. There is no closed-form solution for arithmetic Asian options. MC handles it naturally: simulate the full path, compute the average, calculate the payoff.

def mc_asian_call(S0, K, T, r, sigma, n_steps=252, n_paths=100000):
    """
    Price an Asian call option (arithmetic average).
    Payoff: max(avg(S) - K, 0)
    No closed-form solution exists - MC is the standard approach.
    """
    dt = T / n_steps
    
    # Simulate full paths (need all prices for averaging)
    S = np.zeros((n_steps + 1, n_paths))
    S[0] = S0
    
    for t in range(1, n_steps + 1):
        Z = np.random.standard_normal(n_paths)
        S[t] = S[t-1] * np.exp(
            (r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z
        )
    
    # Average price along each path (excluding initial price)
    avg_price = np.mean(S[1:], axis=0)
    
    # Payoff based on average
    payoffs = np.maximum(avg_price - K, 0)
    
    price = np.exp(-r * T) * np.mean(payoffs)
    se = np.exp(-r * T) * np.std(payoffs) / np.sqrt(n_paths)
    
    return price, se

asian_price, asian_se = mc_asian_call(100, 105, 0.5, 0.05, 0.2)
print(f"Asian Call: {asian_price:.4f} (+/- {1.96*asian_se:.4f})")
print("Note: Asian < European (averaging reduces volatility)")