Monte Carlo Option Pricing in Python: A Practical Guide

You can derive Black-Scholes on a whiteboard and still freeze the first time someone asks you to put a number on the screen. That gap — between the formula and the working price — is where most quant work actually lives. Monte Carlo is the tool that closes it.

This is a hands-on walkthrough of Monte Carlo option pricing in Python. We’ll price a European call, check it against the Black-Scholes closed form, and — the part most tutorials skip — put an honest error bar on the result. By the end you’ll have code you can run, and a mental model you can extend to the payoffs that don’t have a formula.

Why Monte Carlo at all?

Black-Scholes gives you a European call price in one line. So why simulate anything?

Because the closed form is the exception, not the rule. The moment a payoff depends on the path — an Asian option averaging over time, a barrier that knocks out, a basket of correlated names — the neat formula disappears. Monte Carlo doesn’t care. It prices whatever you can simulate. Learning it on a vanilla call, where you can check your answer, is how you earn the right to trust it on the exotics, where you can’t.

The one idea everything rests on

Risk-neutral pricing says something almost suspiciously simple: the price of an option today is the discounted expected value of its payoff, taken under a specific probability measure — the risk-neutral measure, usually written .

That’s it. Every line of code below is just a way of computing that expectation when we can’t do the integral by hand. (Why that measure and not the real-world one? That’s the change of measure — the single most misunderstood step in the whole subject.)

The model: a stock under the risk-neutral measure

Under , the standard assumption is that the stock follows geometric Brownian motion with drift equal to the risk-free rate:

To price a European option we don’t need the whole path — only the terminal price matters. And GBM has an exact solution for it, so there’s no discretisation error to worry about:

Read that off: draw a standard normal, plug it in, and you have one possible value of the stock at expiry. Do it a million times and you’ve sampled the distribution of under .

The algorithm, in four lines of thought

1. Draw standard-normal samples.
2. Turn each into a terminal price with the formula above.
3. Compute each discounted payoff: for a call.
4. Average them. That average is your price.

The code: Monte Carlo option pricing in Python

Here’s the whole thing in NumPy. It’s deliberately plain — no classes, no cleverness, just the four steps.

import numpy as np

# Market and contract
S0    = 100.0    # spot
K     = 100.0    # strike
r     = 0.05     # risk-free rate (continuously compounded)
sigma = 0.20     # volatility
T     = 1.0      # maturity in years
n_sims = 1_000_000

# 1. draw standard normals   2. terminal prices under Q
np.random.seed(42)
Z  = np.random.standard_normal(n_sims)
ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)

# 3. discounted call payoffs   4. average
discounted = np.exp(-r * T) * np.maximum(ST - K, 0.0)
price   = discounted.mean()
std_err = discounted.std(ddof=1) / np.sqrt(n_sims)

print(f"Monte Carlo price: {price:.4f}")
print(f"Standard error:    {std_err:.4f}")
print(f"95% CI: [{price - 1.96*std_err:.4f}, {price + 1.96*std_err:.4f}]")

Run it and you’ll see:

Monte Carlo price: 10.4342
Standard error:    0.0147
95% CI: [10.4053, 10.4630]

A one-line price, and — just as important — a standard error and a confidence interval. We’ll come back to why those two extra lines matter more than the price itself.

Does it actually work? Benchmark against Black-Scholes

A number on its own means nothing. The discipline that separates a practitioner from someone who “ran a simulation” is always checking against something you trust. For a vanilla European call, we’re lucky: the exact answer exists.

from scipy.stats import norm

def black_scholes_call(S0, K, r, sigma, T):
    d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)

print(f"Black-Scholes price: {black_scholes_call(S0, K, r, sigma, T):.4f}")

Which prints:

Black-Scholes price: 10.4506

The Monte Carlo estimate lands right on top of the closed form — 10.4342 against 10.4506, comfortably inside the confidence interval we computed a moment ago. If it didn’t, something would be wrong, and now you’d know before it reached a book. That’s the habit worth building: for every pricer you write, find a benchmark and hold it to that benchmark.

How much should you trust the number? Error bars

Here’s what tutorials tend to gloss over: a Monte Carlo price is a random number. Run it again with a different seed and you’ll get a slightly different answer. So “10.43” isn’t really the output — the output is “10.43, give or take.”

That “give or take” is the standard error, and it shrinks like . Which is both reassuring and annoying: to halve your error you need four times the simulations. Throwing paths at the problem gets expensive fast.

The confidence interval in the code makes this concrete. With a million paths, ours is about — so we can honestly say the price is 10.43, and we know to the cent how sharp that claim is. Reporting a price without it is like quoting a measurement with no units. On a desk, nobody would take you seriously.

Making it tighter for free: antithetic variates

Since more paths get expensive, the smarter move is to make each path count for more. The simplest trick — and a properly useful one — is antithetic variates: for every random draw , also use . The two are negatively correlated, so their errors partly cancel, and your estimate tightens without a single extra payoff evaluation.

np.random.seed(42)
Z = np.random.standard_normal(n_sims // 2)

def discounted_payoff(z):
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * z)
    return np.exp(-r * T) * np.maximum(ST - K, 0.0)

# average each Z with its mirror -Z, then average across pairs
pair_avg = 0.5 * (discounted_payoff(Z) + discounted_payoff(-Z))
price    = pair_avg.mean()
std_err  = pair_avg.std(ddof=1) / np.sqrt(n_sims // 2)

print(f"Antithetic price:  {price:.4f}")
print(f"Standard error:    {std_err:.4f}")

Output:

Antithetic price:  10.4568
Standard error:    0.0104

Same number of payoff evaluations, essentially the same runtime — but the error bar dropped from 0.0147 to 0.0104. That’s the variance roughly halved, for free. There’s a whole toolkit of these variance-reduction techniques — control variates, importance sampling — and knowing which to reach for is part of what makes a pricer production-grade rather than merely correct.

Where this goes next

What you’ve built is the skeleton of every Monte Carlo pricer on a real desk. The bones are the same; production just adds muscle:

• Discounting done properly — not one flat rate, but an OIS curve.
• Path-dependent payoffs — simulate the whole path, not just , and you can price Asians, barriers, lookbacks.
• Calibration — where does even come from? Not thin air; it’s backed out from market prices.
• Validation — the martingale test, convergence checks, put-call parity: the suite that proves the number before it’s trusted with real money.

Each of those is a step from a formula toward a system you’d actually run. That path — from replication and no-arbitrage, through the change of measure, all the way to a validated Monte Carlo pricer with error bars and a full test suite — is exactly what we build, line by line, in Risk-Neutral Valuation. If this walkthrough clicked, the course is the same mindset applied end to end.

Until then: take the code above, change the payoff, break it, benchmark it. That’s how the intuition sticks.

Get the Quant Reference Card — four pages covering stochastic calculus, measure change, the models you’ll actually meet, discretisation, rates, credit and XVA, with the practitioner notes that don’t make it into textbooks. Free.