Monte Carlo methods price derivatives by simulating thousands of random paths for the underlying asset, then averaging the discounted payoffs. The law of large numbers guarantees convergence; variance reduction makes it practical.
Simulating random paths
Under geometric Brownian motion, a stock price evolves as S(t+dt) = S(t) * exp((mu - sigma^2/2)*dt + sigma*sqrt(dt)*Z), where Z is a standard normal draw. Each simulation traces one possible future. Run enough of them and you map the distribution.
Scheme
; Geometric Brownian Motion path simulation; Using Box-Muller transform for normal random numbers
(define (box-muller u1 u2)
(* (sqrt (* -2 (log u1)))
(cos (* 23.141592653589793 u2))))
; Simulate one GBM step; S(t+dt) = S(t) * exp((mu - sigma^2/2)*dt + sigma*sqrt(dt)*Z)
(define (gbm-step s mu sigma dt z)
(* s (exp (+ (* (- mu (* 0.5 sigma sigma)) dt)
(* sigma (sqrt dt) z)))))
; Simulate a path of n steps
(define (gbm-path s0 mu sigma dt n seed)
(define (go s i acc seed-val)
(if (= i n) (reverse acc)
(let* ((u1 (/ (modulo (* seed-val 16807) 2147483647) 2147483647.0))
(u2 (/ (modulo (* (+ seed-val 12345) 16807) 2147483647) 2147483647.0))
(z (box-muller u1 u2))
(s-new (gbm-step s mu sigma dt z)))
(go s-new (+ i 1) (cons s-new acc) (modulo (* seed-val 16807) 2147483647)))))
(go s0 0 (list s0) seed))
(define path (gbm-path 1000.050.20.25442))
(display "GBM path (S0=100, mu=5%, sigma=20%, quarterly): ")
(newline)
(for-each (lambda (p) (display (/ (round (* p 100)) 100)) (display " ")) path)
Python
importmath, random
def gbm_path(s0, mu, sigma, dt, n, seed=42):
random.seed(seed)
path = [s0]
s = s0
for _ inrange(n):
z = random.gauss(0, 1)
s = s * math.exp((mu - 0.5*sigma**2)*dt + sigma*math.sqrt(dt)*z)
path.append(s)
return path
path = gbm_path(100, 0.05, 0.2, 0.25, 4)
print("GBM path (S0=100, mu=5%, sigma=20%, quarterly):")
print([round(p, 2) for p in path])
Pricing European options by simulation
A European call pays max(S_T - K, 0) at expiry. Simulate many terminal prices under the risk-neutral measure (replace mu with r), compute each payoff, discount back, and average. With enough paths, the Monte Carlo estimate converges to the Black-Scholes price.
Scheme
; Monte Carlo European call pricing; Risk-neutral: use r instead of mu
(define (box-muller u1 u2)
(* (sqrt (* -2 (log u1)))
(cos (* 23.141592653589793 u2))))
(define (mc-european-call s0 k r sigma t n-paths seed)
(define (simulate i sum-payoff seed-val)
(if (= i n-paths)
(* (exp (* (- r) t)) (/ sum-payoff n-paths))
(let* ((u1 (/ (modulo (* seed-val 16807) 2147483647) 2147483647.0))
(u2 (/ (modulo (* (+ seed-val 99991) 16807) 2147483647) 2147483647.0))
(z (box-muller u1 u2))
(st (* s0 (exp (+ (* (- r (* 0.5 sigma sigma)) t)
(* sigma (sqrt t) z)))))
(payoff (max 0 (- st k))))
(simulate (+ i 1) (+ sum-payoff payoff)
(modulo (* seed-val 16807) 2147483647)))))
(simulate 00.0 seed))
(define price (mc-european-call 1001050.050.21.01000042))
(display "MC European call price (10k paths): ")
(display (/ (round (* price 100)) 100))
(newline)
(display "Black-Scholes analytical: ~8.02")
Python
importmath, random
def mc_european_call(s0, k, r, sigma, t, n_paths, seed=42):
random.seed(seed)
payoffs = []
for _ inrange(n_paths):
z = random.gauss(0, 1)
st = s0 * math.exp((r - 0.5*sigma**2)*t + sigma*math.sqrt(t)*z)
payoffs.append(max(0, st - k))
returnmath.exp(-r*t) * sum(payoffs) / n_paths
price = mc_european_call(100, 105, 0.05, 0.2, 1.0, 10000)
print(f"MC European call price (10k paths): {price:.2f}")
print("Black-Scholes analytical: ~8.02")
Variance reduction techniques
Antithetic variates: for every normal draw Z, also use -Z. The two paths are negatively correlated, so their average has lower variance than two independent draws. This cuts the standard error roughly in half for free.
An Asian option pays based on the average price along the path, not just the terminal price. A barrier option knocks in or out if the price crosses a threshold. Both require simulating the full path, not just the endpoint.
Scheme
; Asian call option: payoff = max(0, avg(S) - K); Must simulate the full path to compute the average
(define (box-muller u1 u2)
(* (sqrt (* -2 (log u1)))
(cos (* 23.141592653589793 u2))))
(define (asian-call-path s0 k r sigma t n-steps seed)
(let ((dt (/ t n-steps)))
(define (go i s sum-s seed-val)
(if (= i n-steps)
(let ((avg (/ sum-s n-steps)))
(max 0 (- avg k)))
(let* ((u1 (/ (modulo (* seed-val 16807) 2147483647) 2147483647.0))
(u2 (/ (modulo (* (+ seed-val 77773) 16807) 2147483647) 2147483647.0))
(z (box-muller u1 u2))
(s-new (* s (exp (+ (* (- r (* 0.5 sigma sigma)) dt)
(* sigma (sqrt dt) z))))))
(go (+ i 1) s-new (+ sum-s s-new)
(modulo (* seed-val 16807) 2147483647)))))
(go 0 s0 0.0 seed)))
; Single path example
(define payoff (asian-call-path 1001000.050.211242))
(display "Asian call payoff (one path, 12 monthly steps): ")
(display (/ (round (* payoff 100)) 100))
(newline)
(display "Asian options are cheaper than European โ averaging reduces volatility")
Python
importmath, random
def mc_asian_call(s0, k, r, sigma, t, n_steps, n_paths, seed=42):
random.seed(seed)
dt = t / n_steps
payoffs = []
for _ inrange(n_paths):
s = s0
path_sum = 0for _ inrange(n_steps):
z = random.gauss(0, 1)
s = s * math.exp((r - 0.5*sigma**2)*dt + sigma*math.sqrt(dt)*z)
path_sum += s
avg = path_sum / n_steps
payoffs.append(max(0, avg - k))
returnmath.exp(-r*t) * sum(payoffs) / n_paths
price = mc_asian_call(100, 100, 0.05, 0.2, 1.0, 12, 10000)
print(f"Asian call price (10k paths, 12 steps): {price:.2f}")
print("Asian < European because averaging reduces effective volatility")
Convergence and standard error
Monte Carlo converges at rate 1/sqrt(N) regardless of dimension. Double the paths, cut the error by sqrt(2). The standard error of the estimate is the sample standard deviation divided by sqrt(N). This tells you how many paths you need for a given precision.