Importance sampling for rare events, antithetic variates, control variates, stratified sampling, and stacking all three. Achieves 100-10,000x variance reduction over crude Monte Carlo. Table stakes for production simulation.
When crude Monte Carlo is too slow or imprecise. Specifically:
When a contract trades at $0.003, crude MC at 100K samples produces zero or one hit. Importance sampling tilts the distribution to oversample the rare region, then corrects with a likelihood ratio.
import numpy as np
def rare_event_IS(S0, K_crash, sigma, T, N_paths=100_000):
"""
Importance sampling for extreme downside binary contracts.
Example: P(S&P drops 20% in one week)
"""
K = S0 * (1 - K_crash)
mu_original = -0.5 * sigma**2
log_threshold = np.log(K / S0)
mu_tilt = log_threshold / T
Z = np.random.standard_normal(N_paths)
# Simulate under TILTED measure
log_returns_tilted = mu_tilt * T + sigma * np.sqrt(T) * Z
S_T_tilted = S0 * np.exp(log_returns_tilted)
# Likelihood ratio
log_LR = (
-0.5 * ((log_returns_tilted - mu_original * T) / (sigma * np.sqrt(T)))**2
+ 0.5 * ((log_returns_tilted - mu_tilt * T) / (sigma * np.sqrt(T)))**2
)
LR = np.exp(log_LR)
# IS estimator
payoffs = (S_T_tilted < K).astype(float)
is_estimates = payoffs * LR
p_IS = is_estimates.mean()
se_IS = is_estimates.std() / np.sqrt(N_paths)
# Compare with crude MC
Z_crude = np.random.standard_normal(N_paths)
S_T_crude = S0 * np.exp(mu_original * T + sigma * np.sqrt(T) * Z_crude)
p_crude = (S_T_crude < K).mean()
se_crude = np.sqrt(p_crude * (1 - p_crude) / N_paths) if p_crude > 0 else float('inf')
return {
'p_IS': p_IS, 'se_IS': se_IS,
'p_crude': p_crude, 'se_crude': se_crude,
'variance_reduction': (se_crude / se_IS)**2 if se_IS > 0 else float('inf')
}
result = rare_event_IS(S0=5000, K_crash=0.20, sigma=0.15, T=5/252)
print(f"IS estimate: {result['p_IS']:.6f} +/- {result['se_IS']:.6f}")
print(f"Crude estimate: {result['p_crude']:.6f} +/- {result['se_crude']:.6f}")
print(f"Variance reduction: {result['variance_reduction']:.1f}x")
On extreme contracts, IS achieves 100–10,000x variance reduction. 100 IS samples can beat 1,000,000 crude samples.
For monotone payoffs (all binary contracts), pair each draw Z with -Z:
import numpy as np
def antithetic_binary(S0, K, sigma, T, N_paths=100_000):
"""Antithetic variates for binary contracts. ~50-75% variance reduction, zero cost."""
N_half = N_paths // 2
Z = np.random.standard_normal(N_half)
S_T_pos = S0 * np.exp((-0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
S_T_neg = S0 * np.exp((-0.5 * sigma**2) * T + sigma * np.sqrt(T) * (-Z))
payoff_pos = (S_T_pos > K).astype(float)
payoff_neg = (S_T_neg > K).astype(float)
paired = (payoff_pos + payoff_neg) / 2
p_hat = paired.mean()
se = paired.std() / np.sqrt(N_half)
return p_hat, se
p, se = antithetic_binary(S0=100, K=105, sigma=0.20, T=30/365)
print(f"Antithetic estimate: {p:.6f} +/- {se:.6f}")
Partition the probability space into J strata, sample within each, combine.
import numpy as np
from scipy.stats import norm
def stratified_binary_mc(S0, K, sigma, T, J=10, N_total=100_000):
"""Stratified MC for binary contract pricing."""
n_per_stratum = N_total // J
estimates = []
for j in range(J):
U = np.random.uniform(j/J, (j+1)/J, n_per_stratum)
Z = norm.ppf(U)
S_T = S0 * np.exp((-0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)
stratum_mean = (S_T > K).mean()
estimates.append(stratum_mean)
p_stratified = np.mean(estimates)
se_stratified = np.std(estimates) / np.sqrt(J)
return p_stratified, se_stratified
p, se = stratified_binary_mc(S0=100, K=105, sigma=0.20, T=30/365)
print(f"Stratified estimate: {p:.6f} +/- {se:.6f}")
Use Black-Scholes digital price (closed form) as a control variate for stochastic volatility simulations:
p_cv = p_hat_SV - beta * (p_hat_BS_MC - p_BS_exact)
where beta = Cov(payoff_SV, payoff_BS) / Var(payoff_BS), estimated from the same paths.
Antithetic variates inside each stratum, with a control variate correction. Routine 100–500x variance reduction over crude MC. This is table stakes in production.