Pricing an option with Monte Carlo

Monte Carlo simulation is one of the most useful tools a quant has: when a quantity is hard to compute directly, you simulate the randomness many times and average the results. In this lesson you'll use it to price an option, then check your answer against an exact formula. By the end you will have the core of a portfolio project.

What we're pricing

A European call option gives its holder the right, but not the obligation, to buy a stock at a fixed strike price $K$ on a fixed expiry date $T$ years from now. If the stock ends above the strike you exercise and pocket the difference; otherwise the option expires worthless. So the payoff at expiry is

$$\text{payoff} = \max(S_T - K,\ 0),$$

where $S_T$ is the stock price at expiry. Today we only know $S_0$, the current price. The whole problem is that $S_T$ is random.

A model for the stock price

The standard starting point is geometric Brownian motion. For pricing, we use it under the risk-neutral measure, where every asset is expected to grow at the risk-free rate $r$:

$$dS_t = r\,S_t\,dt + \sigma\,S_t\,dW_t .$$

Here $\sigma$ is the volatility (how much the price wobbles, per year) and $W_t$ is a Brownian motion. Why $r$ and not the stock's real expected return? Because an option can be hedged with the stock itself, its fair price cannot depend on anyone's view of where the stock is heading. Pricing in a world where everything grows at $r$ gives the arbitrage-free answer.

This equation has an exact solution, which is what makes the simulation so simple:

$$S_T = S_0 \exp\!\Big(\big(r - \tfrac{1}{2}\sigma^2\big)T + \sigma\sqrt{T}\,Z\Big), \qquad Z \sim \mathcal{N}(0, 1).$$

One standard normal draw $Z$ gives one possible future price $S_T$. The $-\tfrac{1}{2}\sigma^2$ term is a correction that keeps the expected price at $S_0 e^{rT}$.

The Monte Carlo estimator

The fair price of the option is its expected payoff, discounted back to today:

$$C = e^{-rT}\,\mathbb{E}\big[\max(S_T - K,\ 0)\big].$$

We can't compute that expectation by hand easily, but we can estimate it. Draw $N$ independent values $Z_1, \dots, Z_N$, turn each into a terminal price $S_T^{(i)}$, and average the discounted payoffs:

$$\hat{C}_N = \frac{e^{-rT}}{N} \sum_{i=1}^{N} \max\big(S_T^{(i)} - K,\ 0\big).$$

By the law of large numbers, $\hat{C}_N \to C$ as $N \to \infty$.

How accurate is it?

An estimate without an error bar is only half an answer. If $s$ is the sample standard deviation of the discounted payoffs, the standard error of the estimate is

$$\mathrm{SE} = \frac{s}{\sqrt{N}},$$

and by the central limit theorem a 95% confidence interval is roughly $\hat{C}_N \pm 1.96\, \mathrm{SE}$. Notice the $\sqrt{N}$: to halve the error you need four times as many paths. That slow convergence is why so much quant research goes into variance reduction.

Checking against Black–Scholes

For this particular option there is an exact answer, the Black–Scholes formula:

$$C = S_0\,\Phi(d_1) - K e^{-rT}\,\Phi(d_2),$$ $$d_1 = \frac{\ln(S_0/K) + \big(r + \tfrac{1}{2}\sigma^2\big)T}{\sigma\sqrt{T}}, \qquad d_2 = d_1 - \sigma\sqrt{T},$$

where $\Phi$ is the standard normal cumulative distribution function. With $S_0 = 100$, $K = 100$, $r = 5\%$, $\sigma = 20\%$ and $T = 1$ year, it gives $C \approx 10.45$.

Having an exact answer to compare against is what makes this a great first project: you can prove your simulation works before using the same technique on options with no formula at all.

Try it

Run the simulation with the parameters above. Watch how the error bar shrinks as the number of paths grows, and how the true price stays inside it.

Build it yourself

Here is the whole method in a few lines of Python with NumPy:

import numpy as np

def monte_carlo_call(s0, k, r, sigma, t, n_paths, seed=0):
    rng = np.random.default_rng(seed)
    z = rng.standard_normal(n_paths)
    s_t = s0 * np.exp((r - 0.5 * sigma**2) * t + sigma * np.sqrt(t) * z)
    discounted = np.exp(-r * t) * np.maximum(s_t - k, 0.0)
    price = discounted.mean()
    standard_error = discounted.std(ddof=1) / np.sqrt(n_paths)
    return price, standard_error

price, se = monte_carlo_call(100, 100, 0.05, 0.2, 1.0, 100_000)
print(f"{price:.3f} ± {1.96 * se:.3f}")

Make it a portfolio project

A pricer that matches Black–Scholes is a solid start. To turn it into something worth showing, extend it and write up what you find:

• Variance reduction. Add antithetic variates (use both $Z$ and $-Z$) and a control variate, and measure how much each one shrinks the standard error for the same number of paths.
• Convergence plot. Plot the estimate and its confidence interval against $N$ on a log scale, with the exact price as a reference line.
• Options with no formula. Price an Asian option (payoff based on the average price) or a barrier option. These need the full path, not just $S_T$, so simulate the steps in between.
• Greeks. Estimate delta, the sensitivity of the price to $S_0$, with finite differences and compare it with the Black–Scholes delta $\Phi(d_1)$.
• A clear write-up. Explain the model, show your checks against the exact answer, and say honestly where the method breaks down. Clear communication is what reviewers notice.

Key takeaways

• Monte Carlo turns an expectation you can't compute into an average you can.
• Price under the risk-neutral measure: drift $r$, not the stock's real expected return.
• Always report a standard error; it shrinks like $1/\sqrt{N}$.
• Validate against a case with a known answer before trusting the method on harder ones.