Monte Carlo Simulation in Python
Monte Carlo simulation uses repeated random sampling to estimate numerical results that would be difficult or impossible to calculate analytically. In Python, you can build a working Monte Carlo simulation in under 20 lines of code using NumPy - and that same approach scales to pricing exotic derivatives, estimating portfolio risk, and solving problems across quantitative finance.
This tutorial walks through every step, from generating random numbers to building production-quality simulations. Each section includes complete, runnable Python code. If you want the theoretical background first, see our Monte Carlo simulation overview. This article is purely hands-on.
You'll need Python 3.9+ and three libraries:
pip install numpy scipy matplotlib
All the code below uses these imports:
import numpy as np from scipy import stats import matplotlib.pyplot as plt
Random Number Generation and Reproducibility
Every Monte Carlo simulation starts with random numbers. NumPy's numpy.random module provides the generators you'll use constantly - normal distributions for stock returns, uniform distributions for sampling, and the seed mechanism that makes your results reproducible.
# Set a seed for reproducibility rng = np.random.default_rng(seed=42) # Draw 10 standard normal random numbers z = rng.standard_normal(10) print(z) # [-0.1397, 0.6474, 1.0826, -1.1528, -0.2646, ...] # Draw from N(mu, sigma^2) mu, sigma = 0.05, 0.2 returns = rng.normal(loc=mu, scale=sigma, size=10_000) print(f"Mean: {returns.mean():.4f}, Std: {returns.std():.4f}") # Mean: 0.0498, Std: 0.2002 # Uniform random numbers on [0, 1) u = rng.uniform(size=10_000)
A few things to note. First, always use np.random.default_rng() rather than the legacy np.random.seed() API. The new generator API (introduced in NumPy 1.17) is faster, statistically superior, and avoids global state issues. Second, always set a seed when developing and testing - it means you and anyone reviewing your code get identical results on every run. Third, use size= to generate arrays in one call rather than looping. This vectorised approach is orders of magnitude faster and is the foundation of performant Monte Carlo in Python.
Simulating Stock Prices with Geometric Brownian Motion
The standard model for simulating stock prices is geometric Brownian motion (GBM). Under GBM, the stock price at each time step evolves as:
S(t + dt) = S(t) * exp((mu - sigma^2 / 2) * dt + sigma * sqrt(dt) * Z)
where mu is the expected annual return, sigma is the annual volatility, dt is the time step, and Z is a standard normal random variable.
Here's a complete Monte Carlo simulation that generates 10,000 stock price paths over one year:
def simulate_gbm(S0, mu, sigma, T, dt, n_paths, seed=42): """Simulate stock price paths using geometric Brownian motion.""" rng = np.random.default_rng(seed) n_steps = int(T / dt) # Generate all random shocks at once (n_steps x n_paths) Z = rng.standard_normal((n_steps, n_paths)) # Calculate log returns for each step drift = (mu - 0.5 * sigma**2) * dt diffusion = sigma * np.sqrt(dt) * Z log_returns = drift + diffusion # Build price paths by cumulative sum of log returns log_paths = np.vstack([np.zeros(n_paths), np.cumsum(log_returns, axis=0)]) paths = S0 * np.exp(log_paths) return paths # Parameters S0 = 100 # Initial stock price mu = 0.08 # 8% expected annual return sigma = 0.25 # 25% annual volatility T = 1.0 # 1 year dt = 1/252 # Daily steps paths = simulate_gbm(S0, mu, sigma, T, dt, n_paths=10_000) print(f"Shape: {paths.shape}") # (253, 10000) print(f"Final prices - Mean: {paths[-1].mean():.2f}, Std: {paths[-1].std():.2f}")
The key to performance is generating the entire matrix of random numbers in one call (rng.standard_normal((n_steps, n_paths))) and using np.cumsum to build paths without loops. This Monte Carlo simulation python example runs 10,000 paths with daily steps in under a second on any modern machine.
To visualise a sample of paths:
time_grid = np.linspace(0, T, paths.shape[0]) plt.figure(figsize=(10, 6)) plt.plot(time_grid, paths[:, :50], alpha=0.3, linewidth=0.5) plt.xlabel("Time (years)") plt.ylabel("Stock Price (£)") plt.title("Monte Carlo Stock Price Simulation - 50 of 10,000 Paths") plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()
The fan shape shows the range of possible outcomes widening over time - a direct consequence of the sigma * sqrt(dt) term compounding across steps. If you're interested in the mathematical underpinnings, our guide to Python for finance covers the numerical computing fundamentals in more detail.
Pricing a European Call Option
Monte Carlo option pricing works by simulating many stock price paths, computing the payoff of each path at expiry, and discounting the average payoff back to today. For a European call option the payoff is max(S_T - K, 0).
Under risk-neutral pricing, we replace the real-world drift mu with the risk-free rate r:
def price_european_call_mc(S0, K, r, sigma, T, n_paths=100_000, seed=42): """Price a European call option using Monte Carlo simulation.""" rng = np.random.default_rng(seed) # Simulate terminal stock prices (no need for full paths) Z = rng.standard_normal(n_paths) ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z) # Compute discounted payoffs payoffs = np.maximum(ST - K, 0) price = np.exp(-r * T) * payoffs.mean() # Standard error se = np.exp(-r * T) * payoffs.std() / np.sqrt(n_paths) return price, se # Option parameters S0, K, r, sigma, T = 100, 105, 0.05, 0.25, 1.0 mc_price, mc_se = price_european_call_mc(S0, K, r, sigma, T) print(f"Monte Carlo price: £{mc_price:.4f} (SE: {mc_se:.4f})")
Notice we don't need intermediate time steps here - only the terminal price S_T matters for a European option. This makes the simulation much faster: one random draw per path instead of 252.
Let's compare our Monte Carlo estimate against the exact Black-Scholes formula:
def black_scholes_call(S0, K, r, sigma, T): """Exact Black-Scholes price for a European call.""" d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T)) d2 = d1 - sigma * np.sqrt(T) return S0 * stats.norm.cdf(d1) - K * np.exp(-r * T) * stats.norm.cdf(d2) bs_price = black_scholes_call(S0, K, r, sigma, T) print(f"Black-Scholes price: £{bs_price:.4f}") print(f"MC error: £{abs(mc_price - bs_price):.4f}")
With 100,000 paths, the Monte Carlo price typically lands within a few pence of the exact answer. The standard error tells you the precision of your estimate - it shrinks proportionally to 1 / sqrt(n_paths), so quadrupling the number of paths halves the error.
Pricing Exotic Options
Monte Carlo simulation really earns its keep with exotic options - instruments whose payoffs depend on the entire price path, not just the terminal value. These options rarely have closed-form solutions, making simulation the primary pricing method.
Asian Option (Arithmetic Average)
An Asian call option pays max(A - K, 0) where A is the arithmetic average stock price over the option's life. Because the average of log-normal variables isn't log-normal, there's no exact analytical formula.
def price_asian_call_mc(S0, K, r, sigma, T, n_steps=252, n_paths=100_000, seed=42): """Price an arithmetic average Asian call option.""" rng = np.random.default_rng(seed) dt = T / n_steps Z = rng.standard_normal((n_steps, n_paths)) log_returns = (r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z log_paths = np.cumsum(log_returns, axis=0) paths = S0 * np.exp(log_paths) # Arithmetic average across all time steps for each path avg_prices = paths.mean(axis=0) payoffs = np.maximum(avg_prices - K, 0) price = np.exp(-r * T) * payoffs.mean() se = np.exp(-r * T) * payoffs.std() / np.sqrt(n_paths) return price, se asian_price, asian_se = price_asian_call_mc(S0=100, K=105, r=0.05, sigma=0.25, T=1.0) print(f"Asian call price: £{asian_price:.4f} (SE: {asian_se:.4f})")
Up-and-Out Barrier Option
A barrier option becomes worthless (or activates) if the stock price crosses a specified level during the option's life. An up-and-out call pays max(S_T - K, 0) only if the stock never exceeds the barrier B:
def price_up_and_out_call_mc(S0, K, B, r, sigma, T, n_steps=252, n_paths=100_000, seed=42): """Price an up-and-out barrier call option.""" rng = np.random.default_rng(seed) dt = T / n_steps Z = rng.standard_normal((n_steps, n_paths)) log_returns = (r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z log_paths = np.cumsum(log_returns, axis=0) paths = S0 * np.exp(log_paths) # Check if any path crossed the barrier knocked_out = np.any(paths >= B, axis=0) # Terminal payoff, zeroed out for knocked-out paths ST = paths[-1] payoffs = np.maximum(ST - K, 0) payoffs[knocked_out] = 0.0 price = np.exp(-r * T) * payoffs.mean() se = np.exp(-r * T) * payoffs.std() / np.sqrt(n_paths) return price, se barrier_price, barrier_se = price_up_and_out_call_mc( S0=100, K=100, B=130, r=0.05, sigma=0.25, T=1.0 ) print(f"Up-and-out call price: £{barrier_price:.4f} (SE: {barrier_se:.4f})")
One subtlety with barrier options: discrete monitoring (checking the barrier only at each time step) underestimates the true continuous barrier crossing probability. More time steps improve accuracy but cost more computation. For production pricing, continuity corrections or Brownian bridge techniques address this bias.
Estimating Value at Risk with Monte Carlo
Monte Carlo Value at Risk (VaR) works by simulating thousands of possible portfolio returns, then reading off the loss at the desired percentile. It's more flexible than parametric VaR because you can model any return distribution and capture non-linear risk from options or other complex instruments.
def monte_carlo_var(S0, mu, sigma, T, confidence=0.95, n_paths=100_000, seed=42): """Estimate VaR and CVaR for a single-asset portfolio.""" rng = np.random.default_rng(seed) # Simulate terminal portfolio values Z = rng.standard_normal(n_paths) ST = S0 * np.exp((mu - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z) # Calculate P&L pnl = ST - S0 # VaR: the loss at the (1 - confidence) percentile var = -np.percentile(pnl, (1 - confidence) * 100) # CVaR (Expected Shortfall): average loss beyond VaR tail_losses = pnl[pnl <= -var] cvar = -tail_losses.mean() return var, cvar, pnl # £100,000 equity position, 10-day horizon portfolio_value = 100_000 daily_vol = 0.02 T_days = 10 var_95, cvar_95, pnl = monte_carlo_var( S0=portfolio_value, mu=0.0, sigma=daily_vol * np.sqrt(252), T=T_days / 252, confidence=0.95 ) print(f"95% 10-day VaR: £{var_95:,.2f}") print(f"95% 10-day CVaR: £{cvar_95:,.2f}")
You can also visualise the P&L distribution with the VaR threshold:
plt.figure(figsize=(10, 6)) plt.hist(pnl, bins=200, density=True, alpha=0.7, color="steelblue", edgecolor="none") plt.axvline(-var_95, color="red", linestyle="--", linewidth=2, label=f"95% VaR: £{var_95:,.0f}") plt.axvline(-cvar_95, color="darkred", linestyle="--", linewidth=2, label=f"95% CVaR: £{cvar_95:,.0f}") plt.xlabel("Profit / Loss (£)") plt.ylabel("Density") plt.title("Monte Carlo P&L Distribution with VaR and CVaR") plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.show()
The red dashed line marks the VaR threshold - 95% of simulated outcomes fall to the right of it. The darker line shows CVaR (also called Expected Shortfall), which averages all losses beyond VaR and gives a better picture of tail risk.
Portfolio Simulation with Correlated Assets
Real portfolios contain multiple assets whose returns are correlated. To simulate correlated random variables, you use the Cholesky decomposition of the correlation (or covariance) matrix. The Cholesky decomposition finds a lower triangular matrix L such that the covariance matrix C = L * L^T. Multiplying independent standard normal draws by L produces correlated draws with the desired covariance structure.
def simulate_correlated_portfolio( S0_vec, mu_vec, sigma_vec, corr_matrix, T, dt, n_paths=10_000, seed=42 ): """Simulate correlated GBM paths for multiple assets.""" rng = np.random.default_rng(seed) n_assets = len(S0_vec) n_steps = int(T / dt) # Build covariance matrix from correlation and volatilities D = np.diag(sigma_vec) cov_matrix = D @ corr_matrix @ D # Cholesky decomposition L = np.linalg.cholesky(cov_matrix) # Generate independent normals, then correlate Z_indep = rng.standard_normal((n_steps, n_paths, n_assets)) Z_corr = Z_indep @ L.T # (n_steps, n_paths, n_assets) # Simulate each asset drift = (mu_vec - 0.5 * sigma_vec**2) * dt diffusion = np.sqrt(dt) * Z_corr log_returns = drift + diffusion log_paths = np.cumsum(log_returns, axis=0) paths = S0_vec * np.exp(log_paths) return paths # Three-asset portfolio S0 = np.array([100, 50, 200]) mu = np.array([0.08, 0.06, 0.10]) sigma = np.array([0.20, 0.15, 0.30]) corr = np.array([ [1.0, 0.5, 0.3], [0.5, 1.0, 0.2], [0.3, 0.2, 1.0] ]) weights = np.array([0.4, 0.3, 0.3]) paths = simulate_correlated_portfolio(S0, mu, sigma, corr, T=1.0, dt=1/252) print(f"Shape: {paths.shape}") # (252, 10000, 3) # Portfolio value over time initial_investment = 100_000 shares = (weights * initial_investment) / S0 portfolio_paths = np.sum(paths * shares, axis=2)
Let's compute VaR for this multi-asset portfolio:
# Terminal portfolio P&L terminal_values = portfolio_paths[-1] pnl = terminal_values - initial_investment var_99 = -np.percentile(pnl, 1) cvar_99 = -pnl[pnl <= -var_99].mean() print(f"Portfolio 99% 1-year VaR: £{var_99:,.2f}") print(f"Portfolio 99% 1-year CVaR: £{cvar_99:,.2f}") print(f"Mean return: {pnl.mean() / initial_investment:.2%}")
The Cholesky approach preserves the correlation structure across all simulated paths. If you set all off-diagonal correlations to zero, the assets move independently. Set them to 1.0, and they move in lockstep. Real correlation matrices estimated from market data sit somewhere in between and often shift during crises - something to bear in mind when interpreting the results.
Variance Reduction Techniques
Standard Monte Carlo converges slowly - halving the standard error requires four times as many paths. Variance reduction techniques improve accuracy without increasing computation. Here are the two most practical methods.
Antithetic Variates
For every random draw Z, also use -Z. Because Z and -Z are negatively correlated, the average of their payoffs has lower variance than two independent draws.
def price_call_antithetic(S0, K, r, sigma, T, n_paths=50_000, seed=42): """European call with antithetic variance reduction.""" rng = np.random.default_rng(seed) Z = rng.standard_normal(n_paths) # Original paths ST_pos = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z) payoff_pos = np.maximum(ST_pos - K, 0) # Antithetic paths (using -Z) ST_neg = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * (-Z)) payoff_neg = np.maximum(ST_neg - K, 0) # Average paired payoffs, then take the mean paired_payoffs = 0.5 * (payoff_pos + payoff_neg) price = np.exp(-r * T) * paired_payoffs.mean() se = np.exp(-r * T) * paired_payoffs.std() / np.sqrt(n_paths) return price, se # Compare: standard MC with 100k paths vs antithetic with 50k pairs (same total draws) std_price, std_se = price_european_call_mc(S0, K, r, sigma, T, n_paths=100_000) anti_price, anti_se = price_call_antithetic(S0, K, r, sigma, T, n_paths=50_000) print(f"Standard MC: £{std_price:.4f} (SE: {std_se:.4f})") print(f"Antithetic: £{anti_price:.4f} (SE: {anti_se:.4f})") print(f"SE reduction: {(1 - anti_se/std_se):.1%}")
You'll typically see a 20 - 40% reduction in standard error for free - no extra random numbers, no extra computation.
Control Variates
Control variates use a related quantity with a known expected value to reduce variance. For option pricing, a natural control variate is the stock price itself, whose expected value under the risk-neutral measure is S0 * exp(r * T).
def price_call_control_variate(S0, K, r, sigma, T, n_paths=100_000, seed=42): """European call with control variate variance reduction.""" rng = np.random.default_rng(seed) Z = rng.standard_normal(n_paths) ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z) payoffs = np.maximum(ST - K, 0) discount = np.exp(-r * T) # Control variate: the discounted stock price # Known mean of ST under risk-neutral measure expected_ST = S0 * np.exp(r * T) # Estimate optimal coefficient beta cov_matrix = np.cov(payoffs, ST) beta = cov_matrix[0, 1] / cov_matrix[1, 1] # Adjusted payoffs adjusted = payoffs - beta * (ST - expected_ST) price = discount * adjusted.mean() se = discount * adjusted.std() / np.sqrt(n_paths) return price, se cv_price, cv_se = price_call_control_variate(S0, K, r, sigma, T) print(f"Control var: £{cv_price:.4f} (SE: {cv_se:.4f})") print(f"SE reduction vs standard: {(1 - cv_se/std_se):.1%}")
Control variates are especially powerful when you have a highly correlated analytical benchmark. For European options, combining both antithetic variates and control variates can reduce standard error by 70% or more compared to naive Monte Carlo - the equivalent of running ten times as many simulations.
Convergence and Accuracy
How many simulations do you actually need? The answer depends on the precision you require and the variance of the quantity you're estimating.
The standard error of a Monte Carlo estimate decreases as 1 / sqrt(N), where N is the number of paths. This means:
| Paths | Relative SE | Typical Price Error |
|---|---|---|
| 1,000 | 1.000 | £0.30 - £0.50 |
| 10,000 | 0.316 | £0.10 - £0.15 |
| 100,000 | 0.100 | £0.03 - £0.05 |
| 1,000,000 | 0.032 | £0.01 - £0.02 |
You can plot convergence empirically to see this in action:
def convergence_plot(S0, K, r, sigma, T, max_paths=200_000, seed=42): """Show how the MC estimate converges as paths increase.""" rng = np.random.default_rng(seed) Z = rng.standard_normal(max_paths) ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z) payoffs = np.maximum(ST - K, 0) * np.exp(-r * T) # Cumulative mean and SE at various sample sizes checkpoints = np.arange(100, max_paths + 1, 100) means = np.array([payoffs[:n].mean() for n in checkpoints]) ses = np.array([payoffs[:n].std() / np.sqrt(n) for n in checkpoints]) bs_price = black_scholes_call(S0, K, r, sigma, T) plt.figure(figsize=(10, 6)) plt.plot(checkpoints, means, linewidth=0.8, color="steelblue", label="MC estimate") plt.fill_between( checkpoints, means - 1.96 * ses, means + 1.96 * ses, alpha=0.2, color="steelblue", label="95% CI" ) plt.axhline(bs_price, color="red", linestyle="--", label=f"Black-Scholes: £{bs_price:.4f}") plt.xlabel("Number of Paths") plt.ylabel("Option Price (£)") plt.title("Monte Carlo Convergence") plt.legend() plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() convergence_plot(S0, K, r, sigma, T)
The plot shows the Monte Carlo estimate bouncing around the true Black-Scholes price, with the 95% confidence interval (blue band) narrowing as N increases. By 50,000 paths the estimate is stable to within a few pence for a standard European option.
For confidence intervals, the 95% CI for your Monte Carlo estimate is:
n_paths = 100_000 mc_price, mc_se = price_european_call_mc(S0, K, r, sigma, T, n_paths=n_paths) ci_lower = mc_price - 1.96 * mc_se ci_upper = mc_price + 1.96 * mc_se print(f"Price: £{mc_price:.4f}, 95% CI: [£{ci_lower:.4f}, £{ci_upper:.4f}]")
If you need tighter bounds, either increase N or apply the variance reduction techniques from the previous section.
Performance Tips
Monte Carlo simulations are embarrassingly parallel and respond well to optimisation. Here are the techniques that make the biggest difference in Python.
Vectorise Everything with NumPy
The single most important rule: never loop over paths in Python. NumPy operations execute in compiled C, making them 100 - 1000x faster than equivalent Python loops.
import time n = 1_000_000 rng = np.random.default_rng(42) # Slow: Python loop start = time.perf_counter() results_loop = [] for _ in range(n): z = rng.standard_normal() st = 100 * np.exp(0.05 * 1.0 + 0.25 * z) results_loop.append(max(st - 105, 0)) loop_time = time.perf_counter() - start # Fast: vectorised NumPy rng2 = np.random.default_rng(42) start = time.perf_counter() Z = rng2.standard_normal(n) ST = 100 * np.exp(0.05 * 1.0 + 0.25 * Z) results_vec = np.maximum(ST - 105, 0) vec_time = time.perf_counter() - start print(f"Loop: {loop_time:.3f}s") print(f"Vectorised: {vec_time:.3f}s") print(f"Speedup: {loop_time / vec_time:.0f}x")
On a typical machine in 2026, you'll see a 200 - 500x speedup. This isn't a micro-optimisation - it's the difference between a simulation taking 10 minutes and one taking under a second.
Use Multiprocessing for Large Simulations
When you genuinely need millions of paths with many time steps, split the work across CPU cores:
from multiprocessing import Pool from functools import partial def _simulate_chunk(chunk_seed, S0, K, r, sigma, T, n_paths): """Worker function for parallel MC pricing.""" rng = np.random.default_rng(chunk_seed) Z = rng.standard_normal(n_paths) ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z) payoffs = np.maximum(ST - K, 0) return payoffs.sum(), (payoffs**2).sum(), n_paths def price_call_parallel(S0, K, r, sigma, T, total_paths=1_000_000, n_workers=4): """Price a European call using multiple CPU cores.""" paths_per_worker = total_paths // n_workers seeds = [1000 + i for i in range(n_workers)] worker_fn = partial( _simulate_chunk, S0=S0, K=K, r=r, sigma=sigma, T=T, n_paths=paths_per_worker ) with Pool(n_workers) as pool: results = pool.map(worker_fn, seeds) total_sum = sum(r[0] for r in results) total_sq_sum = sum(r[1] for r in results) total_n = sum(r[2] for r in results) mean_payoff = total_sum / total_n var_payoff = total_sq_sum / total_n - mean_payoff**2 price = np.exp(-r * T) * mean_payoff se = np.exp(-r * T) * np.sqrt(var_payoff / total_n) return price, se par_price, par_se = price_call_parallel(S0, K, r, sigma, T) print(f"Parallel MC: £{par_price:.4f} (SE: {par_se:.4f})")
Each worker gets a different seed to ensure independent random streams. The results are aggregated using running sums to avoid transferring large arrays between processes. On a 4-core machine, this gives close to a 4x speedup for computationally heavy simulations.
Memory Considerations
For path-dependent options that require full price histories, memory can become the bottleneck before CPU time. A simulation with 1,000,000 paths and 252 time steps uses about 2 GB of memory (float64). If you're hitting memory limits:
- Process paths in batches and accumulate results
- Use float32 instead of float64 (halves memory, negligible impact on accuracy for most applications)
- Only store what you need - for a barrier option, you can track the running maximum instead of storing the full path
def price_barrier_memory_efficient(S0, K, B, r, sigma, T, n_steps=252, n_paths=1_000_000, batch_size=100_000, seed=42): """Memory-efficient barrier option pricing using batches.""" rng = np.random.default_rng(seed) dt = T / n_steps total_payoff = 0.0 total_payoff_sq = 0.0 for batch_start in range(0, n_paths, batch_size): actual_batch = min(batch_size, n_paths - batch_start) # Track running price and max (no full path storage) S = np.full(actual_batch, S0) S_max = np.full(actual_batch, S0) for t in range(n_steps): Z = rng.standard_normal(actual_batch) S = S * np.exp((r - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * Z) S_max = np.maximum(S_max, S) knocked_out = S_max >= B payoffs = np.maximum(S - K, 0) payoffs[knocked_out] = 0.0 total_payoff += payoffs.sum() total_payoff_sq += (payoffs**2).sum() mean_payoff = total_payoff / n_paths var_payoff = total_payoff_sq / n_paths - mean_payoff**2 price = np.exp(-r * T) * mean_payoff se = np.exp(-r * T) * np.sqrt(var_payoff / n_paths) return price, se eff_price, eff_se = price_barrier_memory_efficient(S0=100, K=100, B=130, r=0.05, sigma=0.25, T=1.0) print(f"Barrier (1M paths, batched): £{eff_price:.4f} (SE: {eff_se:.4f})")
This approach uses under 10 MB of memory regardless of how many total paths you simulate.
Putting It All Together
Here's a complete, self-contained script that prices a European call using all the techniques from this tutorial - GBM simulation, antithetic variates, control variates, and confidence intervals:
import numpy as np from scipy import stats def monte_carlo_pricer(S0, K, r, sigma, T, n_paths=200_000, seed=42): """ Full-featured European call pricer with variance reduction. Returns price, standard error, and 95% confidence interval. """ rng = np.random.default_rng(seed) Z = rng.standard_normal(n_paths) # Antithetic variates Z_all = np.concatenate([Z, -Z]) ST_all = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z_all) payoffs_all = np.maximum(ST_all - K, 0) # Pair antithetic paths payoff_pairs = 0.5 * (payoffs_all[:n_paths] + payoffs_all[n_paths:]) # Control variate adjustment ST_pairs = 0.5 * (ST_all[:n_paths] + ST_all[n_paths:]) expected_ST = S0 * np.exp(r * T) cov_est = np.cov(payoff_pairs, ST_pairs) beta = cov_est[0, 1] / cov_est[1, 1] adjusted = payoff_pairs - beta * (ST_pairs - expected_ST) discount = np.exp(-r * T) price = discount * adjusted.mean() se = discount * adjusted.std() / np.sqrt(n_paths) ci = (price - 1.96 * se, price + 1.96 * se) return price, se, ci # Run it price, se, ci = monte_carlo_pricer(S0=100, K=105, r=0.05, sigma=0.25, T=1.0) bs = black_scholes_call(100, 105, 0.05, 0.25, 1.0) print(f"MC Price: £{price:.4f}") print(f"Std Error: £{se:.6f}") print(f"95% CI: [£{ci[0]:.4f}, £{ci[1]:.4f}]") print(f"Black-Scholes: £{bs:.4f}") print(f"Absolute Error: £{abs(price - bs):.6f}")
Frequently Asked Questions
How many Monte Carlo simulations should I run?
It depends on the precision you need. For rough estimates (within 1% of the true value), 10,000 paths usually suffice. For pricing standard options to within a penny, 100,000 - 500,000 paths is typical. For exotic options with discontinuous payoffs (barrier options, digital options), you may need 1,000,000+ paths. The standard error formula SE = std(payoffs) / sqrt(N) tells you exactly how many you need: decide your target SE, estimate the payoff standard deviation from a pilot run, and solve for N.
Is Monte Carlo simulation slower than analytical formulas?
Yes, significantly. A Black-Scholes calculation takes microseconds. A Monte Carlo estimate with 100,000 paths takes milliseconds. For instruments with closed-form solutions, analytical formulas are always preferable for speed. Monte Carlo earns its place when no closed-form solution exists - exotic options, path-dependent payoffs, complex correlation structures, and multi-factor models. In practice, most quant teams use analytical formulas where they can and fall back to Monte Carlo for everything else.
Can I use Monte Carlo simulation for American options?
Standard Monte Carlo can't directly price American options because the optimal exercise decision at each time step depends on the continuation value, which is what you're trying to estimate. The Longstaff-Schwartz algorithm (Least Squares Monte Carlo, or LSMC) solves this by regressing continuation values onto basis functions of the current stock price at each time step. It works well in practice and is the standard approach for pricing American and Bermudan options with Monte Carlo. That said, for single-asset American options, binomial trees or finite difference methods are usually faster and more accurate.
What are the limitations of Monte Carlo simulation in finance?
Monte Carlo simulation is only as good as the model driving it. If you assume GBM but the real stock has stochastic volatility, jumps, or regime changes, your prices will be biased regardless of how many paths you run. Convergence is slow - the 1 / sqrt(N) rate means going from 1% precision to 0.1% precision requires 100x more computation. High-dimensional problems (large portfolios, many risk factors) amplify this. And for risk management, Monte Carlo estimates of tail quantiles (99.9% VaR, for example) are inherently noisy because only a tiny fraction of paths fall into the tail.
Should I use NumPy or a library like QuantLib for Monte Carlo?
For learning and custom implementations, NumPy is the right choice. You understand every line, you can modify anything, and the performance is excellent for vectorised code. For production systems that need tested, validated implementations of complex models (stochastic volatility, multi-curve discounting, exotic payoffs with early exercise), QuantLib provides battle-tested engines out of the box. Many quant teams prototype in NumPy and deploy with QuantLib or an internal C++ library. In 2026, the NumPy ecosystem (including JAX and CuPy for GPU acceleration) has closed much of the performance gap for pure Monte Carlo workloads.
Want to go deeper on Monte Carlo Simulation in Python: Step-by-Step Tutorial 2026?
This article covers the essentials, but there's a lot more to learn. Inside Quantt, you'll find hands-on coding exercises, interactive quizzes, and structured lessons that take you from fundamentals to production-ready skills — across 50+ courses in technology, finance, and mathematics.
Free to get started · No credit card required