What is Monte Carlo Integration?
Monte Carlo (MC) methods estimate quantities by repeated random sampling. The name comes from the famous casino in Monaco, and was coined by mathematicians during the Manhattan Project. At its core, MC integration replaces a difficult integral with an average over random samples:
Monte Carlo estimate of integral:
I = ?_O f(x) dx V(O) (1/N) S f(x?)
where x? are N uniform random samples in domain O with volume V(O).
Error (standard error of mean):
s_MC = s_f / vN (s_f = standard deviation of f over the domain)
The key result: error drops as 1/vN regardless of dimension d.
This beats deterministic quadrature (error ? N^{-2/d}) for d = 5.
p Estimation
The classic MC demonstration: sample (x,y) uniformly in [0,1] and test whether x+y=1 (inside the quarter-circle). The fraction of accepted points converges to p/4:
p 4 (points inside quarter-circle) / (total points)
With N=10,000 samples: typical error 4s/vN 40.5/100 0.02
With N=1,000,000 samples: typical error 0.002
The 1/vN law means 100 more samples gives only 10 less error.
Variance Reduction Techniques
Stratified Sampling
Divide the domain into N equal strata and place exactly one sample per stratum. This guarantees uniform coverage and eliminates the "clumping and gap" artefacts of pure random sampling, reducing variance by up to a factor of N for smooth functions.
Quasi-Random (Low-Discrepancy) Sequences
The Halton sequence generates deterministic points that fill the unit hypercube far more uniformly than pseudo-random numbers. The discrepancy (maximum deviation from uniform) falls as O(log(N)/N) instead of O(1/vN), giving faster convergence for smooth integrands.
Halton sequence in base b:
H(n, b) = reverse_base_b_digits(n) / b^k
Halton(1,2) = 0.5, Halton(2,2) = 0.25, Halton(3,2) = 0.75
Halton(1,3) = 0.333, Halton(2,3) = 0.667, Halton(3,3) = 0.111
2D Halton uses bases (2,3) x-coordinates and y-coordinates
are independently well-distributed.
Option Pricing (Black-Scholes Monte Carlo)
Under the Black-Scholes model, a stock follows Geometric Brownian Motion. The MC price of a European call option with strike K, maturity T, risk-free rate r, volatility s is:
Stock path (one-step GBM):
S(T) = S0 exp((r - s/2)T + svT Z) where Z ~ N(0,1)
Call payoff: max(S(T) - K, 0)
MC price: C e^{-rT} (1/N) S max(S?(T) - K, 0)
Typical: S0=100, K=105, T=1yr, r=5%, s=20%
Black-Scholes closed form: C $8.02
MC converges to this as N ? 8
Applications
| Application | What is sampled | Typical N needed |
|---|
| p estimation | Points in unit square | 106 (error 0.002) |
| High-dim integral (d=10) | Uniform hypercube points | 104106 |
| Option pricing | Stock price paths (GBM) | 104105 |
| Radiation transport | Photon paths, scattering events | 106108 |
| Path tracing (rendering) | Light ray bounces per pixel | 644096/pixel |
| MCMC (Bayesian inference) | Posterior distribution samples | 10105 |
| Polymer physics (folding) | Monomer configuration space | 10610? |
| Nuclear reactor simulation | Neutron collision chains | 10810 |
Frequently Asked Questions
How does Monte Carlo estimate p?
Random points (x,y) in the unit square: if x+y=1 the point is inside the quarter-circle. Area(quarter-circle)/Area(square) = p/4, so p 4(inside count)/(total). Error follows 1/vN: 10,000 samples gives error 0.02, one million samples gives 0.002.
What does stratified sampling do?
Instead of placing points fully at random, stratified sampling divides the domain into equal cells and draws one sample per cell. This prevents clustering and guarantees coverage. For a smooth integrand the variance decreases as N? (vs N? for pure random), making it dramatically faster to converge.
What are quasi-random (Halton) sequences?
Halton sequences are deterministic but fill space uniformly by reversing the base-b digit expansion of integers. For base 2: 1/2, 1/4, 3/4, 1/8, 5/8, 3/8, 7/8, The discrepancy drops as O(log(N)/N), much faster than the O(1/vN) of pseudo-random sequences. This makes quasi-MC methods significantly better for smooth low-dimensional integrands.
How does the option pricing preset work?
Each sample simulates a stock price under geometric Brownian motion: S(T) = S0exp((r-s/2)T + svTZ) where Z~N(0,1). The call payoff max(S(T)-K,0) discounted at e^{-rT} is averaged over N samples. With S0=100, K=105, r=5%, s=20%, T=1yr, the Black-Scholes formula gives C$8.02 the MC estimate converges to this as N increases.