There is a class of problem in mathematics and applied science where the analytical solution either does not exist, or exists but is so computationally intractable that it provides no practical value. High-dimensional integrals. Stochastic differential equations with complex boundary conditions. Option pricing under jump-diffusion. Insurance loss distributions with fat tails and correlated perils.
For this class of problem, Monte Carlo simulation is often not just the best available tool it is the only practical one.
I've used Monte Carlo methods throughout my career: in actuarial pricing models for takaful products, in quantitative trading strategy evaluation, in flood risk modelling for insurance premium calibration, and in stress-testing data pipelines under synthetic load. In each case, the core idea is the same approximate a difficult deterministic quantity by averaging over a large number of random samples but the details of how you sample, how many samples you need, and how you reduce variance differ substantially across domains.
This post is my attempt to write the guide I wished existed when I first started building these models seriously. It's long. It covers the theory, the mathematics, the Python, and the domain applications. Get a coffee.
Part 1: What Monte Carlo Actually Is
The name comes from the Monaco casino district, coined by Stanisław Ulam and John von Neumann during the Manhattan Project in the 1940s. Ulam, recovering from an illness, was playing solitaire and found himself wondering: what is the probability that a random deal of cards produces a successful game? He couldn't calculate it analytically. So he thought about playing hundreds of games and tracking the fraction that succeeded.
The same instinct estimate a probability or expectation by simulating many trials and averaging the outcomes underlies every Monte Carlo method.
Formally, suppose we want to compute:
θ = E[f(X)] = ∫ f(x) · p(x) dx
where X is a random variable with density p(x) and f is some function. If we can draw N independent samples x₁, x₂, ..., xN from p(x), then the Monte Carlo estimator is:
θ̂ = (1/N) · Σᵢ f(xᵢ)
By the Law of Large Numbers, θ̂ → θ as N → ∞. By the Central Limit Theorem, the error is:
θ̂ - θ ~ N(0, σ²/N)
where σ² = Var[f(X)]. The standard error of the estimator is σ/√N.
This convergence rate O(1/√N) is simultaneously Monte Carlo's greatest weakness and greatest strength. It's slow: to halve the error, you need four times as many samples. But crucially, it's independent of dimension. Numerical integration methods like Gaussian quadrature converge much faster in one dimension, but their convergence degrades exponentially as dimensions increase. Monte Carlo's O(1/√N) rate holds regardless of whether you're integrating over 2 dimensions or 2,000. That's why it dominates high-dimensional problems.
Part 2: The Foundation Generating Random Samples
Everything in Monte Carlo depends on the quality and structure of your random samples. Let's build from the ground up.
Uniform Random Numbers
The foundation is the uniform distribution on [0, 1]. NumPy's np.random.default_rng() uses the PCG64 generator a modern, statistically robust pseudorandom number generator that passes all standard randomness tests.
import numpy as np
rng = np.random.default_rng(seed=42) # always seed for reproducibility
u = rng.uniform(0, 1, size=1_000_000)
Always use default_rng() rather than the legacy np.random.rand(). The new interface is faster, more statistically sound, and supports explicit seeding without polluting global state.
The Inverse Transform Method
If we can draw uniform samples, we can draw samples from any distribution whose CDF we can invert.
The Inverse Transform Theorem states: if U ~ Uniform(0,1) and F is a CDF, then X = F⁻¹(U) has distribution F.
import numpy as np
from scipy import stats
rng = np.random.default_rng(42)
N = 1_000_000
# Exponential distribution via inverse transform
# F(x) = 1 - exp(-λx) → F⁻¹(u) = -ln(1-u) / λ
lam = 2.0
u = rng.uniform(0, 1, N)
x_exp = -np.log(1 - u) / lam
# Verify: should match scipy's exponential
x_scipy = stats.expon(scale=1/lam).rvs(N, random_state=42)
print(f"Inverse transform mean: {x_exp.mean():.4f} (expected {1/lam:.4f})")
print(f"Inverse transform std: {x_exp.std():.4f} (expected {1/lam:.4f})")
The Box-Muller Transform (Normal Samples)
For Gaussian samples, the Box-Muller transform converts two uniform samples into two independent standard normals:
Z₁ = √(-2 ln U₁) · cos(2π U₂)
Z₂ = √(-2 ln U₁) · sin(2π U₂)
u1 = rng.uniform(0, 1, N)
u2 = rng.uniform(0, 1, N)
z1 = np.sqrt(-2 * np.log(u1)) * np.cos(2 * np.pi * u2)
z2 = np.sqrt(-2 * np.log(u1)) * np.sin(2 * np.pi * u2)
# z1 and z2 are independent N(0,1)
print(f"z1 mean: {z1.mean():.4f}, z1 std: {z1.std():.4f}")
In practice, use rng.standard_normal() which is optimised internally. But understanding Box-Muller matters when you need to generate correlated normals from a known correlation structure.
Generating Correlated Multivariate Normals (Cholesky Decomposition)
This is essential in any multi-asset or multi-risk simulation. Given a correlation matrix Σ, the Cholesky decomposition Σ = LLᵀ lets you transform independent standard normals into correlated ones:
import numpy as np
rng = np.random.default_rng(42)
N = 100_000
# Suppose we're simulating 3 correlated asset returns
# Correlation matrix
corr = np.array([
[1.00, 0.65, 0.30],
[0.65, 1.00, 0.45],
[0.30, 0.45, 1.00]
])
# Volatilities (annualised)
vols = np.array([0.18, 0.25, 0.35]) # 18%, 25%, 35%
# Covariance matrix
cov = np.diag(vols) @ corr @ np.diag(vols)
# Cholesky decomposition
L = np.linalg.cholesky(cov)
# Generate N independent standard normal vectors
Z = rng.standard_normal((3, N))
# Apply Cholesky: X = L @ Z gives correlated samples
X = L @ Z # shape (3, N) - each column is one scenario
print("Empirical correlation matrix:")
print(np.corrcoef(X).round(3))
# Should be close to the input correlation matrix
The Cholesky approach is the workhorse of multi-risk simulation. It's used in everything from equity portfolio VaR to multi-peril catastrophe models in reinsurance.
Part 3: Option Pricing The Classic Application
Let me ground the theory with the problem that cemented Monte Carlo's place in quantitative finance: pricing options under the Black-Scholes model.
Geometric Brownian Motion
Under risk-neutral measure, a stock price S follows Geometric Brownian Motion (GBM):
dS = r·S·dt + σ·S·dW
where r is the risk-free rate, σ is volatility, and W is a Wiener process. The exact solution over time horizon T is:
S(T) = S(0) · exp[(r - σ²/2)·T + σ·√T·Z]
where Z ~ N(0,1). This is the log-normal stock price model.
European Call Option
A European call with strike K pays max(S(T) K, 0) at maturity. Its price is the discounted expected payoff under the risk-neutral measure:
C = e^(-rT) · E[max(S(T) - K, 0)]
import numpy as np
from scipy import stats
def black_scholes_call(S0, K, r, sigma, T):
"""Analytical Black-Scholes formula for reference."""
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)
def mc_european_call(S0, K, r, sigma, T, N=1_000_000, seed=42):
"""Monte Carlo pricing of a European call option."""
rng = np.random.default_rng(seed)
# Simulate terminal stock prices
Z = rng.standard_normal(N)
ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
# Payoffs
payoffs = np.maximum(ST - K, 0)
# Discounted expected payoff
price = np.exp(-r*T) * payoffs.mean()
std_error = np.exp(-r*T) * payoffs.std() / np.sqrt(N)
return price, std_error
# Parameters
S0 = 100 # spot price
K = 105 # strike price
r = 0.05 # risk-free rate
sigma = 0.20 # volatility
T = 1.0 # 1 year to maturity
# Price
mc_price, mc_se = mc_european_call(S0, K, r, sigma, T)
bs_price = black_scholes_call(S0, K, r, sigma, T)
print(f"Black-Scholes: {bs_price:.4f}")
print(f"Monte Carlo: {mc_price:.4f} ± {mc_se:.4f} (95% CI: [{mc_price-1.96*mc_se:.4f}, {mc_price+1.96*mc_se:.4f}])")
With 1 million paths, the MC estimate will typically be within 1-2 cents of the analytical price. But the power of Monte Carlo is not in pricing vanilla options where the analytical formula exists it's in pricing the ones where no formula does.
Asian Option (Path-Dependent)
An Asian option pays based on the average price over the option's life, not just the terminal price. No closed-form solution exists for arithmetic Asian options under GBM. Monte Carlo handles it naturally:
def mc_asian_call(S0, K, r, sigma, T, n_steps=252, N=100_000, seed=42):
"""
Monte Carlo pricing of an arithmetic Asian call option.
Payoff: max(mean(S(t₁), ..., S(tₙ)) - K, 0)
"""
rng = np.random.default_rng(seed)
dt = T / n_steps
# Simulate full price paths: shape (N, n_steps)
Z = rng.standard_normal((N, n_steps))
# Incremental log-returns
log_returns = (r - 0.5*sigma**2)*dt + sigma*np.sqrt(dt)*Z
# Cumulative sum gives log(S(t)/S0), then exponentiate
log_paths = np.cumsum(log_returns, axis=1)
paths = S0 * np.exp(log_paths) # shape (N, n_steps)
# Arithmetic average price for each path
avg_price = paths.mean(axis=1) # shape (N,)
# Payoffs
payoffs = np.maximum(avg_price - K, 0)
price = np.exp(-r*T) * payoffs.mean()
std_error = np.exp(-r*T) * payoffs.std() / np.sqrt(N)
return price, std_error
price, se = mc_asian_call(S0=100, K=100, r=0.05, sigma=0.20, T=1.0)
print(f"Asian Call Price: {price:.4f} ± {se:.4f}")
The same framework extends to barrier options, lookback options, Bermudan options (with least-squares regression for early exercise), and any payoff you can compute from a simulated price path.
Part 4: Variance Reduction The Real Engineering
Raw Monte Carlo with N=1,000,000 samples gives standard error O(1/1000). For many applications, that's adequate. But when you need tighter estimates pricing exotic derivatives, estimating tail risk at 99.97% confidence, calibrating a reinsurance treaty you need variance reduction techniques to extract more accuracy from fewer samples.
This is where the real engineering of Monte Carlo lives.
4.1 Antithetic Variates
The simplest and most widely applicable technique. For every random draw Z, also simulate with -Z. Because Z and -Z are negatively correlated, their average payoffs are also negatively correlated, which reduces variance.
def mc_call_antithetic(S0, K, r, sigma, T, N=500_000, seed=42):
"""
European call with antithetic variates.
Uses N/2 pairs of (Z, -Z) paths - same total paths as standard with N.
"""
rng = np.random.default_rng(seed)
n_pairs = N // 2
Z = rng.standard_normal(n_pairs)
# Original paths
ST_pos = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
# Antithetic paths
ST_neg = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*(-Z))
payoffs_pos = np.maximum(ST_pos - K, 0)
payoffs_neg = np.maximum(ST_neg - K, 0)
# Average the paired payoffs
payoffs = 0.5 * (payoffs_pos + payoffs_neg)
price = np.exp(-r*T) * payoffs.mean()
std_error = np.exp(-r*T) * payoffs.std() / np.sqrt(n_pairs)
return price, std_error
# Compare
mc_price, mc_se = mc_european_call(S0=100, K=105, r=0.05, sigma=0.20, T=1.0, N=1_000_000)
av_price, av_se = mc_call_antithetic(S0=100, K=105, r=0.05, sigma=0.20, T=1.0, N=1_000_000)
print(f"Standard MC: {mc_price:.4f} ± {mc_se:.5f}")
print(f"Antithetic: {av_price:.4f} ± {av_se:.5f}")
print(f"Variance reduction: {((mc_se/av_se)**2 - 1)*100:.1f}%")
For at-the-money options, antithetic variates typically achieves 30-70% variance reduction equivalent to running 1.5-3x as many standard paths.
4.2 Control Variates
A more powerful technique. If you have a control variate Y a quantity correlated with your target f(X) whose expectation E[Y] is known analytically you can use it to correct the MC estimate.
The estimator is:
θ̂_cv = (1/N) Σᵢ [f(Xᵢ) - c·(Yᵢ - E[Y])]
The optimal coefficient is c* = Cov[f(X), Y] / Var[Y], which minimises variance.
For an Asian option, the geometric average Asian option has a closed-form price and is highly correlated with the arithmetic average. Use it as a control variate:
from scipy import stats
def geometric_asian_call_price(S0, K, r, sigma, T, n):
"""
Analytical price of geometric average Asian call.
Rogers & Shi (1995) formula.
"""
sigma_hat = sigma * np.sqrt((2*n + 1) / (6*(n + 1)))
mu_hat = 0.5 * (r - 0.5*sigma**2 + sigma_hat**2)
d1 = (np.log(S0/K) + (mu_hat + 0.5*sigma_hat**2)*T) / (sigma_hat*np.sqrt(T))
d2 = d1 - sigma_hat*np.sqrt(T)
return np.exp(-r*T) * (S0 * np.exp(mu_hat*T) * stats.norm.cdf(d1)
- K * stats.norm.cdf(d2))
def mc_asian_call_cv(S0, K, r, sigma, T, n_steps=252, N=100_000, seed=42):
"""
Arithmetic Asian call using geometric Asian as control variate.
"""
rng = np.random.default_rng(seed)
dt = T / n_steps
Z = rng.standard_normal((N, n_steps))
log_returns = (r - 0.5*sigma**2)*dt + sigma*np.sqrt(dt)*Z
log_paths = np.cumsum(log_returns, axis=1)
paths = S0 * np.exp(log_paths)
# Arithmetic payoffs (target)
arith_avg = paths.mean(axis=1)
payoffs_arith = np.maximum(arith_avg - K, 0)
# Geometric payoffs (control variate)
geo_avg = np.exp(np.log(paths).mean(axis=1))
payoffs_geo = np.maximum(geo_avg - K, 0)
# Known expectation of geometric payoff
E_geo = geometric_asian_call_price(S0, K, r, sigma, T, n_steps) * np.exp(r*T)
# Optimal coefficient c*
cov_matrix = np.cov(payoffs_arith, payoffs_geo)
c_star = cov_matrix[0, 1] / cov_matrix[1, 1]
# Control variate estimate
payoffs_cv = payoffs_arith - c_star * (payoffs_geo - E_geo)
price = np.exp(-r*T) * payoffs_cv.mean()
std_error = np.exp(-r*T) * payoffs_cv.std() / np.sqrt(N)
return price, std_error, c_star
price_raw, se_raw = mc_asian_call(S0=100, K=100, r=0.05, sigma=0.20, T=1.0, N=100_000)
price_cv, se_cv, c = mc_asian_call_cv(S0=100, K=100, r=0.05, sigma=0.20, T=1.0, N=100_000)
print(f"Raw MC: {price_raw:.4f} ± {se_raw:.5f}")
print(f"Control Variate: {price_cv:.4f} ± {se_cv:.5f}")
print(f"Variance reduction: {((se_raw/se_cv)**2 - 1)*100:.1f}%")
print(f"Optimal c*: {c:.4f}")
Control variates can achieve variance reductions of 90%+ in favourable cases effectively multiplying your sample efficiency by 10x or more.
4.3 Importance Sampling
When estimating rare event probabilities tail risk, default probability, catastrophic loss most standard Monte Carlo samples contribute nothing because the rare event almost never occurs in normal sampling.
Importance sampling changes the sampling distribution to oversample the region of interest, then corrects for the change using likelihood ratios.
If we sample X from a proposal distribution q(x) instead of p(x), the estimator becomes:
θ̂_IS = (1/N) · Σᵢ f(xᵢ) · [p(xᵢ) / q(xᵢ)]
The ratio w(x) = p(x)/q(x) is the likelihood ratio or importance weight.
For tail risk estimation say, estimating P(loss > L) for a large loss threshold we shift the mean of the sampling distribution toward the tail:
def tail_probability_importance_sampling(threshold, mu=0, sigma=1, shift=None, N=100_000, seed=42):
"""
Estimate P(X > threshold) for X ~ N(mu, sigma²) using importance sampling.
Optimal shift: sample from N(threshold, sigma²) instead of N(mu, sigma²).
"""
rng = np.random.default_rng(seed)
# Optimal shift: centre proposal at threshold
if shift is None:
shift = threshold
# Sample from shifted distribution
X_shifted = rng.normal(shift, sigma, N)
# Indicator: did the event occur?
indicator = (X_shifted > threshold).astype(float)
# Likelihood ratio: p(x) / q(x) = N(x; mu, sigma) / N(x; shift, sigma)
log_lr = (-(X_shifted - mu)**2 / (2*sigma**2)
+ (X_shifted - shift)**2 / (2*sigma**2))
lr = np.exp(log_lr)
# Importance-weighted estimate
estimate = (indicator * lr).mean()
std_error = (indicator * lr).std() / np.sqrt(N)
return estimate, std_error
# Estimate P(X > 4) for X ~ N(0,1) - true value ≈ 3.167e-5
threshold = 4.0
# Standard MC - needs enormous N to see any events
rng = np.random.default_rng(42)
x_standard = rng.standard_normal(100_000)
p_standard = (x_standard > threshold).mean()
se_standard = np.sqrt(p_standard * (1 - p_standard) / 100_000)
# Importance sampling
p_is, se_is = tail_probability_importance_sampling(threshold, N=100_000)
from scipy.stats import norm
true_prob = 1 - norm.cdf(threshold)
print(f"True probability: {true_prob:.2e}")
print(f"Standard MC: {p_standard:.2e} ± {se_standard:.2e}")
print(f"Importance Sampling: {p_is:.2e} ± {se_is:.2e}")
print(f"IS variance reduction: {((se_standard/se_is)**2):.0f}x")
Importance sampling is indispensable for credit risk (estimating default probability at high-confidence levels), reinsurance (estimating catastrophic loss tail), and operational risk capital calculation.
4.4 Quasi-Monte Carlo (Low-Discrepancy Sequences)
Instead of pseudo-random numbers, use low-discrepancy sequences deterministic sequences designed to fill the sample space more uniformly than random samples. The most common are Sobol sequences and Halton sequences.
from scipy.stats.qmc import Sobol, Halton
from scipy.stats import norm as scipy_norm
def mc_call_sobol(S0, K, r, sigma, T, N=65536, seed=42):
"""European call using Sobol quasi-random sequence."""
# Sobol sampler - dimension 1 (one risk factor)
sampler = Sobol(d=1, scramble=True, seed=seed)
u = sampler.random(N).flatten()
# Transform uniform Sobol to standard normal
Z = scipy_norm.ppf(u)
ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
payoffs = np.maximum(ST - K, 0)
price = np.exp(-r*T) * payoffs.mean()
std_error = np.exp(-r*T) * payoffs.std() / np.sqrt(N)
return price, std_error
bs_price = black_scholes_call(S0=100, K=105, r=0.05, sigma=0.20, T=1.0)
mc_price, mc_se = mc_european_call(S0=100, K=105, r=0.05, sigma=0.20, T=1.0, N=65536)
qmc_price, qmc_se = mc_call_sobol(S0=100, K=105, r=0.05, sigma=0.20, T=1.0, N=65536)
print(f"Black-Scholes: {bs_price:.6f}")
print(f"Standard MC: {mc_price:.6f} (error: {abs(mc_price-bs_price):.6f})")
print(f"Sobol QMC: {qmc_price:.6f} (error: {abs(qmc_price-bs_price):.6f})")
QMC can achieve convergence rates approaching O(1/N) rather than O(1/√N) a massive improvement. In practice, for smooth integrands in moderate dimensions (up to ~30), Sobol sequences consistently outperform pseudo-random Monte Carlo.
Part 5: Value at Risk and Expected Shortfall
In market risk, the two headline metrics are Value at Risk (VaR) and Expected Shortfall (ES).
- VaR at confidence level α: the loss
Lsuch thatP(loss > L) = 1 - α. At 99% VaR, you lose more thanLonly 1% of the time. - Expected Shortfall (ES): the expected loss given that the loss exceeds VaR. Also called Conditional VaR (CVaR). ES is preferred by regulators (Basel III/IV) because it accounts for the severity of tail losses, not just the threshold.
import numpy as np
from scipy import stats
def portfolio_var_es(
weights, # portfolio weights
mu, # expected annual returns
cov, # covariance matrix
horizon=10, # holding period in days
confidence=0.99,
N=500_000,
seed=42
):
"""
Monte Carlo VaR and Expected Shortfall for a multi-asset portfolio.
Returns VaR and ES as positive numbers (loss amounts).
"""
rng = np.random.default_rng(seed)
n_assets = len(weights)
dt = horizon / 252 # convert days to years
# Cholesky decomposition of covariance matrix
L = np.linalg.cholesky(cov)
# Simulate correlated returns over holding period
Z = rng.standard_normal((n_assets, N))
# Correlated returns: shape (n_assets, N)
returns = mu[:, None] * dt + L @ (np.sqrt(dt) * Z)
# Portfolio P&L for each scenario (as fraction of portfolio value)
portfolio_pnl = weights @ returns # shape (N,)
# Convert to losses (negative returns = losses)
losses = -portfolio_pnl
# VaR: quantile of loss distribution
var = np.quantile(losses, confidence)
# ES: expected loss beyond VaR
es = losses[losses >= var].mean()
# Standard errors via bootstrap
var_se = stats.sem(losses) / stats.norm.pdf(stats.norm.ppf(confidence))
return {
'VaR': var,
'ES': es,
'confidence': confidence,
'horizon_days': horizon,
'loss_distribution': losses
}
# Example: 3-asset portfolio
weights = np.array([0.50, 0.30, 0.20]) # 50% equity, 30% bonds, 20% commodities
mu = np.array([0.08, 0.03, 0.06]) # expected annual returns
corr = np.array([
[1.00, 0.20, 0.40],
[0.20, 1.00, 0.10],
[0.40, 0.10, 1.00]
])
vols = np.array([0.18, 0.06, 0.25])
cov = np.diag(vols) @ corr @ np.diag(vols)
result = portfolio_var_es(weights, mu, cov)
print(f"10-day 99% VaR: {result['VaR']*100:.2f}% of portfolio value")
print(f"10-day 99% ES: {result['ES']*100:.2f}% of portfolio value")
Plotting the Loss Distribution
import matplotlib.pyplot as plt
losses = result['loss_distribution']
var = result['VaR']
es = result['ES']
fig, ax = plt.subplots(figsize=(12, 5))
ax.hist(losses, bins=200, density=True, alpha=0.6, color='steelblue', label='Loss distribution')
ax.axvline(var, color='red', lw=2, linestyle='--', label=f'99% VaR = {var*100:.2f}%')
ax.axvline(es, color='darkred', lw=2, linestyle=':', label=f'99% ES = {es*100:.2f}%')
# Shade the tail
x_tail = losses[losses >= var]
ax.hist(x_tail, bins=50, density=True, alpha=0.4, color='red')
ax.set_xlabel('Portfolio Loss (fraction of portfolio value)')
ax.set_ylabel('Density')
ax.set_title('Monte Carlo Portfolio Loss Distribution - 10-day Holding Period')
ax.legend()
plt.tight_layout()
plt.savefig('loss_distribution.png', dpi=150, bbox_inches='tight')
Part 6: Insurance Aggregate Loss Modelling
This is the domain where Monte Carlo is perhaps most practically useful in financial services. Actuarial science has used simulation for decades, precisely because insurance loss distributions are complicated: they're mixtures, they have heavy tails, they have correlation structures across lines of business, and they need to be aggregated across many policies with varying exposures.
The Collective Risk Model
The collective risk model represents aggregate losses as:
S = X₁ + X₂ + ... + X_N
where N is the (random) number of claims and each Xᵢ is a (random) claim severity. N and the Xᵢ are assumed independent. The distribution of S is a compound distribution.
For most combinations of claim frequency and severity distributions, the aggregate loss distribution has no closed form. Monte Carlo is the standard approach.
import numpy as np
from scipy import stats
def simulate_aggregate_losses(
# Frequency parameters (Poisson claim count)
expected_claims,
# Severity parameters (Lognormal claim size)
severity_mu, # log-mean of claim size
severity_sigma, # log-std of claim size
# Simulation
N=100_000,
seed=42
):
"""
Simulate aggregate insurance losses under the collective risk model.
N: number of claims ~ Poisson(expected_claims)
X_i: claim size ~ LogNormal(severity_mu, severity_sigma)
S = X_1 + ... + X_N
"""
rng = np.random.default_rng(seed)
# Step 1: Simulate claim counts for each policy year
n_claims = rng.poisson(lam=expected_claims, size=N)
# Step 2: For each simulated year, sum that many claim severities
# Vectorised: draw total_claims severities, then split by year
total_claims = n_claims.sum()
all_severities = rng.lognormal(mean=severity_mu, sigma=severity_sigma, size=total_claims)
# Assign severities to policy years using cumulative splits
cumulative = np.concatenate([[0], np.cumsum(n_claims)])
aggregate_losses = np.array([
all_severities[cumulative[i]:cumulative[i+1]].sum()
for i in range(N)
])
return aggregate_losses
def loss_statistics(losses, confidence_levels=[0.95, 0.99, 0.995, 0.999]):
"""Compute key risk metrics from simulated aggregate losses."""
results = {
'mean': losses.mean(),
'std': losses.std(),
'cv': losses.std() / losses.mean(), # coefficient of variation
}
for alpha in confidence_levels:
var = np.quantile(losses, alpha)
es = losses[losses >= var].mean()
results[f'VaR_{int(alpha*100)}'] = var
results[f'ES_{int(alpha*100)}'] = es
return results
# Takaful motor portfolio example:
# - Expected 850 claims per year
# - Lognormal severity with mean RM 8,500 and CV of 1.2
# Lognormal parameters: mu = log(mean) - sigma²/2
mean_severity = 8500
cv_severity = 1.2
sigma_sev = np.sqrt(np.log(1 + cv_severity**2))
mu_sev = np.log(mean_severity) - 0.5*sigma_sev**2
losses = simulate_aggregate_losses(
expected_claims=850,
severity_mu=mu_sev,
severity_sigma=sigma_sev,
N=500_000
)
stats_out = loss_statistics(losses)
print(f"Expected aggregate loss: RM {stats_out['mean']:,.0f}")
print(f"Std deviation: RM {stats_out['std']:,.0f}")
print(f"Coefficient of variation: {stats_out['cv']:.3f}")
print()
print(f"95th percentile (VaR): RM {stats_out['VaR_95']:,.0f}")
print(f"99th percentile (VaR): RM {stats_out['VaR_99']:,.0f}")
print(f"99.5th percentile (VaR): RM {stats_out['VaR_995']:,.0f}")
print(f"99.9th percentile (VaR): RM {stats_out['VaR_999']:,.0f}")
print()
print(f"99% Expected Shortfall: RM {stats_out['ES_99']:,.0f}")
Reinsurance Pricing: Excess-of-Loss Treaties
An excess-of-loss (XL) reinsurance treaty pays losses above a retention R up to a limit L:
Reinsurer pays: max(0, min(S - R, L))
Pricing this correctly requires the full aggregate loss distribution, which Monte Carlo provides directly:
def price_xl_reinsurance(losses, retention, limit):
"""
Price an excess-of-loss reinsurance layer.
losses: array of simulated aggregate losses
retention: cedant retains losses up to this amount
limit: reinsurer covers losses from retention to retention+limit
"""
reinsurer_payments = np.clip(losses - retention, 0, limit)
pure_premium = reinsurer_payments.mean()
loss_on_line = pure_premium / limit # pure premium as % of limit
prob_attachment = (losses > retention).mean()
prob_exhaustion = (losses > retention + limit).mean()
return {
'pure_premium': pure_premium,
'loss_on_line': loss_on_line,
'prob_attachment': prob_attachment,
'prob_exhaustion': prob_exhaustion,
'expected_cedant_loss': np.minimum(losses, retention).mean()
}
# Price an XL layer: RM 1M xs RM 5M (reinsurer covers losses between RM 5M and RM 6M)
xl_result = price_xl_reinsurance(
losses=losses,
retention=5_000_000,
limit=1_000_000
)
print(f"XL Layer: RM 1M xs RM 5M")
print(f"Pure premium: RM {xl_result['pure_premium']:,.0f}")
print(f"Loss on line: {xl_result['loss_on_line']*100:.2f}%")
print(f"Probability of attach:{xl_result['prob_attachment']*100:.2f}%")
print(f"Probability of exhaust:{xl_result['prob_exhaustion']*100:.2f}%")
Part 7: Stochastic Differential Equations Beyond GBM
Real financial and actuarial models rarely follow pure GBM. Here are the most important extensions.
Heston Stochastic Volatility Model
The Heston model allows volatility to be random it follows its own mean-reverting process correlated with the stock price:
dS = r·S·dt + √V·S·dW₁
dV = κ(θ - V)·dt + ξ·√V·dW₂
corr(dW₁, dW₂) = ρ
Parameters: κ (mean reversion speed), θ (long-run variance), ξ (vol-of-vol), ρ (correlation).
def simulate_heston(S0, V0, r, kappa, theta, xi, rho, T, n_steps=252, N=50_000, seed=42):
"""
Simulate stock price paths under the Heston stochastic volatility model.
Uses the Euler-Maruyama discretisation with full truncation scheme.
"""
rng = np.random.default_rng(seed)
dt = T / n_steps
sqrt_dt = np.sqrt(dt)
# Storage
S = np.full(N, S0, dtype=float)
V = np.full(N, V0, dtype=float)
for _ in range(n_steps):
# Correlated Brownian motions
Z1 = rng.standard_normal(N)
Z2 = rng.standard_normal(N)
W1 = Z1
W2 = rho*Z1 + np.sqrt(1 - rho**2)*Z2
# Full truncation: use max(V, 0) in diffusion terms
sqrt_V = np.sqrt(np.maximum(V, 0))
# Update variance (full truncation scheme)
V = (np.maximum(V, 0)
+ kappa*(theta - np.maximum(V, 0))*dt
+ xi*sqrt_V*sqrt_dt*W2)
# Update stock price
S = S * np.exp((r - 0.5*np.maximum(V, 0))*dt + sqrt_V*sqrt_dt*W1)
return S
# Price a call under Heston
S0, K, r, T = 100, 100, 0.05, 1.0
V0 = 0.04 # initial variance (vol = 20%)
kappa = 2.0 # mean reversion speed
theta = 0.04 # long-run variance
xi = 0.3 # vol-of-vol
rho = -0.7 # typical negative correlation (leverage effect)
ST = simulate_heston(S0, V0, r, kappa, theta, xi, rho, T, N=200_000)
payoffs = np.maximum(ST - K, 0)
heston_price = np.exp(-r*T) * payoffs.mean()
bs_price = black_scholes_call(S0, K, r, np.sqrt(theta), T)
print(f"Heston price: {heston_price:.4f}")
print(f"BS (flat vol=20%): {bs_price:.4f}")
print(f"Difference: {heston_price - bs_price:.4f} (Heston smile effect)")
The difference between Heston and Black-Scholes reflects the volatility smile the market observation that implied volatilities are not constant across strikes, which GBM cannot capture.
The Cox-Ingersoll-Ross (CIR) Model for Interest Rates
Interest rate modelling uses mean-reverting processes with non-negativity constraints. The CIR model is:
dr = κ(θ - r)·dt + σ·√r·dW
Used for pricing bonds, interest rate derivatives, and for discount rate simulation in insurance liability valuation.
def simulate_cir(r0, kappa, theta, sigma, T, n_steps=252, N=50_000, seed=42):
"""
Simulate short rate paths under the CIR model.
Full truncation ensures non-negative rates.
"""
rng = np.random.default_rng(seed)
dt = T / n_steps
sqrt_dt = np.sqrt(dt)
rates = np.full((N, n_steps + 1), r0, dtype=float)
for t in range(n_steps):
r = rates[:, t]
Z = rng.standard_normal(N)
rates[:, t+1] = np.maximum(
r + kappa*(theta - r)*dt + sigma*np.sqrt(np.maximum(r, 0))*sqrt_dt*Z,
0
)
return rates
# Zero-coupon bond price via Monte Carlo
r0, kappa, theta, sigma_cir, T = 0.03, 0.5, 0.04, 0.1, 5.0
rate_paths = simulate_cir(r0, kappa, theta, sigma_cir, T, n_steps=252*5, N=100_000)
dt = T / (252*5)
# Bond price = E[exp(-∫r dt)]
integrated_rates = rate_paths[:, :-1].sum(axis=1) * dt
bond_price_mc = np.exp(-integrated_rates).mean()
print(f"5-year ZCB price (MC): {bond_price_mc:.4f}")
Part 8: Convergence Diagnostics and Practical Considerations
Running a Monte Carlo simulation is easy. Running one you can trust requires discipline.
Convergence Monitoring
def mc_convergence_plot(payoffs, true_value=None):
"""Plot running mean and confidence bands to assess convergence."""
N = len(payoffs)
running_mean = np.cumsum(payoffs) / np.arange(1, N+1)
# Running standard error
running_var = np.array([payoffs[:i].var() for i in range(1, N+1, N//100)])
n_vals = np.arange(1, N+1, N//100)
running_se = np.sqrt(running_var / n_vals)
# 95% confidence bands
upper = running_mean[::N//100] + 1.96*running_se
lower = running_mean[::N//100] - 1.96*running_se
# Print summary table
print(f"{'N':>12} {'Estimate':>12} {'Std Error':>12} {'95% CI Width':>14}")
print("-" * 56)
for n, m, se in zip([1000, 10000, 100000, N],
running_mean[[999, 9999, 99999, -1]],
[payoffs[:n].std()/np.sqrt(n) for n in [1000, 10000, 100000, N]]):
print(f"{n:>12,} {m:>12.4f} {se:>12.6f} {3.92*se:>14.6f}")
if true_value:
print(f"\nTrue value: {true_value:.4f}")
print(f"Final error: {abs(running_mean[-1] - true_value):.6f}")
# Demonstrate convergence for European call
S0, K, r, sigma, T = 100, 105, 0.05, 0.20, 1.0
rng = np.random.default_rng(42)
Z = rng.standard_normal(1_000_000)
ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
payoffs = np.exp(-r*T) * np.maximum(ST - K, 0)
mc_convergence_plot(payoffs, true_value=black_scholes_call(S0, K, r, sigma, T))
Sample Size Planning
Before running a simulation, decide how precise you need the estimate to be. If you want a 95% confidence interval of width w, you need:
N ≥ (1.96 · σ / (w/2))²
But you typically don't know σ in advance. Run a pilot simulation with a small N (e.g., 10,000 paths), estimate σ, then compute the required N for your target precision.
def required_sample_size(pilot_std, target_half_width, confidence=0.95):
"""
Compute required N for Monte Carlo to achieve target precision.
pilot_std: std of f(X) from pilot run
target_half_width: desired half-width of confidence interval
"""
z = stats.norm.ppf(0.5 + confidence/2)
N = (z * pilot_std / target_half_width)**2
return int(np.ceil(N))
# Example: how many paths to price the Asian option within 0.01 (1 cent)?
price, se = mc_asian_call(S0=100, K=100, r=0.05, sigma=0.20, T=1.0, N=10_000)
pilot_std = se * np.sqrt(10_000) # recover std from standard error
N_required = required_sample_size(pilot_std, target_half_width=0.01)
print(f"Pilot estimate: {price:.4f} ± {se:.4f}")
print(f"Pilot std: {pilot_std:.4f}")
print(f"Required N for ±0.01 at 95% confidence: {N_required:,}")
Part 9: Why This Matters to Me
I'll close with something more personal.
My day job sits in institutional data pipelines, governance, platform architecture. But the quantitative side of what I do through Guinevere Analytics and in the actuarial and risk work I've done across takaful and financial services keeps pulling me back to simulation.
The reason is this: the hard problems in financial risk are almost always distributional problems. Not "what is the expected value?" but "what does the full distribution of outcomes look like, and what happens in the tail?" Analytical tools answer the first question well. Monte Carlo answers the second.
In takaful, this matters acutely. The aggregate loss distribution of a takaful fund underpins the adequacy of the tabarru' contribution rate. The retakaful structuring decision how much risk to cede, at what retention, with what treaty structure is a question about tail risk. You cannot answer it without simulating the full distribution.
In quantitative investing, strategy evaluation is a Monte Carlo problem. Backtest results are a single realised path from a distribution of possible paths. Running a strategy's logic over thousands of simulated market regimes stress scenarios, synthetic data, bootstrapped historical returns gives you a far more honest picture of what the strategy actually does than any single historical backtest.
And in data engineering, Monte Carlo has an underappreciated application: synthetic data generation for pipeline testing. If you can model the distribution of your input data (claim counts, transaction amounts, sensor readings), you can generate realistic synthetic datasets to stress-test your pipelines before production data arrives. You can test edge cases extreme values, correlated spikes, missing data bursts that may not appear in historical data but should be handled correctly.
The technique is 80 years old. The mathematics is settled. The computational cost has collapsed to near-zero with modern hardware and vectorised NumPy. There is almost no domain in quantitative analytics where Monte Carlo doesn't have a legitimate application.
Learn it properly. The investment compounds.
Appendix: Full Reference Implementation
"""
monte_carlo_reference.py
Complete reference implementation of Monte Carlo tools covered in this post.
"""
import numpy as np
from scipy import stats
from scipy.stats.qmc import Sobol
# ── Random Sampling ────────────────────────────────────────────────────────────
def cholesky_correlated_normals(corr, vols, N, rng):
cov = np.diag(vols) @ corr @ np.diag(vols)
L = np.linalg.cholesky(cov)
Z = rng.standard_normal((len(vols), N))
return L @ Z
# ── Option Pricing ─────────────────────────────────────────────────────────────
def gbm_terminal(S0, r, sigma, T, Z):
return S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
def gbm_paths(S0, r, sigma, T, n_steps, Z):
dt = T / n_steps
log_returns = (r - 0.5*sigma**2)*dt + sigma*np.sqrt(dt)*Z
return S0 * np.exp(np.cumsum(log_returns, axis=1))
def mc_price(payoffs, r, T, N):
discount = np.exp(-r*T)
return discount*payoffs.mean(), discount*payoffs.std()/np.sqrt(N)
# ── Variance Reduction ─────────────────────────────────────────────────────────
def antithetic_pairs(rng, shape):
Z = rng.standard_normal(shape)
return Z, -Z
def control_variate_estimate(target, control, E_control):
cov_mat = np.cov(target, control)
c_star = cov_mat[0,1] / cov_mat[1,1]
return target - c_star*(control - E_control), c_star
# ── Risk Metrics ───────────────────────────────────────────────────────────────
def var_es(losses, confidence):
var = np.quantile(losses, confidence)
es = losses[losses >= var].mean()
return var, es
# ── Aggregate Loss ─────────────────────────────────────────────────────────────
def compound_poisson_lognormal(lam, mu_sev, sigma_sev, N, rng):
n = rng.poisson(lam, N)
total = n.sum()
severities = rng.lognormal(mu_sev, sigma_sev, total)
cumul = np.concatenate([[0], np.cumsum(n)])
return np.array([severities[cumul[i]:cumul[i+1]].sum() for i in range(N)])
def xl_layer(losses, retention, limit):
return np.clip(losses - retention, 0, limit)
All code in this post is written for Python 3.11+ with NumPy 2.x, SciPy 1.13+. For production implementations, consider using QuantLib for financial instruments and actuarial libraries for insurance-specific distributions.