Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Variance Reduction Techniques in Simulation › Stratified Sampling

Michael BrenndoerferDecember 3, 202510 min read

Quantitative FinanceAll chapters

Part 6 of 15Numerical Methods for Derivative Pricing

Chapter 35Variance Reduction Techniques in Simulation › Stratified Sampling

Reading section 5 of 12 · Stratified Sampling · 10 minSection 5 of 12 · 10 min

Section 5 of 12 · 10 min read

Variance Reduction Techniques in Simulation › Stratified SamplingOne section · Sections in this chapter5 / 12 ↓

Sections in this chapter12

Reading tools Definitions & inline links

Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.

Article links

Make inline references clickable

Stratified sampling improves Monte Carlo estimates by ensuring systematic coverage of the random input space. This section explains the approach and demonstrates its application to option pricing.

The Core IntuitionLink Copied

Standard Monte Carlo generates random samples that may cluster in some regions while leaving other regions undersampled. Stratified sampling forces the samples to cover the input space more evenly by dividing it into strata and sampling from each stratum. This ensures that extreme outcomes in both tails are represented, reducing variance caused by sampling fluctuations.

Consider throwing darts randomly at a target. Some areas will have many darts while others have few, purely by chance. If instead you divide the target into regions and ensure exactly one dart lands in each region, the coverage is more uniform and your estimate of "average distance from center" would be more precise. The random clustering that occurs with purely random sampling adds variance to our estimates; stratified sampling eliminates this source of variance.

To understand why clustering causes problems, imagine estimating the expected payoff of a call option using only 100 random samples. By chance, these 100 standard normal random numbers might all fall between -2 and +2, leaving the tails underrepresented. But the tails matter for option pricing because extreme moves determine whether deep out-of-the-money options finish in the money. Stratified sampling ensures that we sample from the entire distribution, including the tails, in a controlled and representative manner.

Mathematical FoundationLink Copied

The mathematical analysis of stratified sampling reveals precisely why forcing uniform coverage reduces variance. Divide the sample space into $K$ non-overlapping strata with probabilities $p_1, p_2, \ldots, p_K$ where $\sum_{k=1}^K p_k = 1$ . The stratified estimator allocates $n_k$ samples to stratum $k$ and computes:

\hat{\mu}_{\text{strat}} = \sum_{k=1}^K p_k \bar{X}_k

where:

$\hat{\mu}_{\text{strat}}$ : stratified estimator
$p_k$ : probability of stratum $k$
$\bar{X}_k$ : sample mean within stratum $k$
$K$ : total number of strata

This formula is a weighted average of the sample means within each stratum, with weights equal to the stratum probabilities. The stratified estimator reconstructs the overall expectation by combining the conditional expectations within each region of the sample space.

The variance of the stratified estimator (with proportional allocation $n_k = n \cdot p_k$ ) is:

\text{Var}(\hat{\mu}_{\text{strat}}) = \frac{1}{n}\sum_{k=1}^K p_k \sigma_k^2

where:

$\text{Var}(\hat{\mu}_{\text{strat}})$ : variance of the stratified estimator
$n$ : total number of samples
$K$ : total number of strata
$p_k$ : probability of stratum $k$
$\sigma_k^2$ : variance of $X$ within stratum $k$

This variance expression contains only the within-stratum variances. To understand the improvement, we must compare this to the variance of standard Monte Carlo. The standard Monte Carlo variance can be decomposed as:

\sigma_X^2 = \sum_{k=1}^K p_k \sigma_k^2 + \sum_{k=1}^K p_k (\mu_k - \mu)^2

where:

$\sigma_X^2$ : total variance of a single sample
$K$ : total number of strata
$p_k$ : probability of stratum $k$
$\sigma_k^2$ : variance of $X$ within stratum $k$
$\mu_k$ : expected value within stratum $k$
$\mu$ : overall expected value

This decomposition, known as the law of total variance separates the total variance into two components. The first term is the within-stratum variance (what stratified sampling has). The second term is the between-stratum variance (what stratified sampling eliminates). Thus, stratified sampling always reduces variance, with the reduction equal to the between-stratum variance.

The between-stratum variance measures how much the stratum means differ from each other. When the function we're estimating varies significantly across the sample space, meaning different regions have very different expected values, the between-stratum variance is large. Stratified sampling eliminates this variation by guaranteeing representation from each region, removing the variance that would otherwise arise from random fluctuations in how many samples fall in each region.

Implementation for Option PricingLink Copied

For option pricing, we stratify the uniform random variable $U$ that generates our normal random variables via the inverse CDF transform. By ensuring $U$ values are spread evenly across $[0,1]$ , we ensure the corresponding normal random variables cover the distribution evenly. This approach is particularly elegant because it works with the standard machinery of random number generation while adding the stratification benefit.

The implementation divides the unit interval into $K$ equal strata, each with probability $1/K$ . Within each stratum, we generate uniform random numbers that fall only within that stratum's portion of the interval. When we apply the inverse normal CDF transformation to these stratified uniform values, the resulting normal values are guaranteed to cover the entire normal distribution evenly.

In[18]:

Code

def monte_carlo_call_stratified(
    S0, K, T, r, sigma, n_simulations, n_strata=100
):
    """
    Price a European call option using stratified sampling.
    """
    np.random.seed(42)

    # Number of samples per stratum
    samples_per_stratum = n_simulations // n_strata

    # Generate stratified uniform samples
    U_stratified = np.zeros(n_strata * samples_per_stratum)
    for k in range(n_strata):
        # Sample uniformly within stratum k
        U_stratum = (k + np.random.uniform(size=samples_per_stratum)) / n_strata
        U_stratified[
            k * samples_per_stratum : (k + 1) * samples_per_stratum
        ] = U_stratum

    # Transform to standard normal
    Z_stratified = norm.ppf(U_stratified)

    # Generate terminal stock prices
    ST = S0 * np.exp(
        (r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z_stratified
    )

    # Compute payoffs
    payoffs = np.exp(-r * T) * np.maximum(ST - K, 0)

    price = np.mean(payoffs)
    std_error = np.std(payoffs) / np.sqrt(len(payoffs))

    return price, std_error

Let's compare stratified sampling to standard Monte Carlo.

In[19]:

Code

# Standard Monte Carlo for comparison
np.random.seed(42)
n_sims_compare = 100000
Z_standard = np.random.standard_normal(n_sims_compare)
ST_standard = S0 * np.exp(
    (r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z_standard
)
payoffs_standard = np.exp(-r * T) * np.maximum(ST_standard - K, 0)
mc_price_std = np.mean(payoffs_standard)
mc_se_std = np.std(payoffs_standard) / np.sqrt(n_sims_compare)

# Stratified sampling
strat_price, strat_se = monte_carlo_call_stratified(
    S0, K, T, r, sigma, n_sims_compare, n_strata=1000
)

Out[20]:

Console

Stratified Sampling Results
==================================================
Black-Scholes Price:      \$10.4506

Standard Monte Carlo:
  Price:                  \$10.4739
  Standard Error:         \$0.0466

Stratified Sampling:
  Price:                  \$10.4500
  Standard Error:         \$0.0465

Variance Reduction Factor: 1.00x

Out[21]:

Visualization

Histogram and rug plot of standard Monte Carlo samples from a normal distribution showing irregular spacing with clusters and gaps. — Standard Monte Carlo sampling: 100 draws from a standard normal distribution showing random clustering and gaps. Some regions of the distribution are accidentally oversampled while others are undersampled, introducing variance in Monte Carlo estimates.

Histogram and rug plot of stratified samples from a normal distribution showing uniform coverage across all equal-probability strata. — Stratified sampling: 100 draws from the same standard normal distribution, allocated proportionally across 10 equal-probability strata with exactly 10 points per stratum. The uniform coverage eliminates variance caused by accidental density fluctuations.

Latin Hypercube SamplingLink Copied

Latin Hypercube Sampling (LHS) extends stratified sampling to multiple dimensions. In standard stratified sampling for multiple variables, the number of strata grows exponentially with dimension (the curse of dimensionality LHS avoids this by ensuring each variable is stratified marginally while allowing joint flexibility.

The challenge with multi-dimensional stratification becomes apparent when we consider path-dependent options. An option that depends on 50 time steps has 50 random inputs. Full stratification would require dividing each dimension into strata, and with even 10 strata per dimension, we would need $10^{50}$ cells to cover the space, an impossibly large number. LHS provides a practical alternative that maintains the benefits of stratification without this exponential explosion.

Latin Hypercube Sampling

LHS divides each input dimension into $n$ equally probable intervals and ensures exactly one sample falls in each interval for each dimension. This guarantees good marginal coverage without the exponential growth in strata that full stratification would require.

The name "Latin Hypercube" comes from the connection to Latin squares in combinatorics. A Latin square is an arrangement where each row and column contains each symbol exactly once. LHS generalizes this to higher dimensions, so when we project the samples onto any single dimension, they cover that dimension evenly with exactly one sample in each stratum.

In[22]:

Code

from scipy.stats import qmc


def latin_hypercube_call(S0, K, T, r, sigma, n_simulations):
    """
    Price a European call using Latin Hypercube Sampling.
    """
    # Generate LHS samples
    sampler = qmc.LatinHypercube(d=1, seed=42)
    lhs_samples = sampler.random(n=n_simulations)

    # Transform to standard normal
    Z_lhs = norm.ppf(lhs_samples).flatten()

    # Generate terminal prices
    ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z_lhs)

    # Compute payoffs
    payoffs = np.exp(-r * T) * np.maximum(ST - K, 0)

    price = np.mean(payoffs)
    std_error = np.std(payoffs) / np.sqrt(n_simulations)

    return price, std_error


lhs_price, lhs_se = latin_hypercube_call(
    S0, K, T, r, sigma, n_simulations=100000
)

Out[23]:

Console

Latin Hypercube Sampling Results
==================================================
Black-Scholes Price:      \$10.4506
LHS Price:                \$10.4505
LHS Standard Error:       \$0.0465
Variance Reduction vs MC: 1.00x

Visualizing Sampling CoverageLink Copied

The following visualization compares how different sampling methods cover the random input space.

In[24]:

Code

n_points = 500

# Standard Monte Carlo
np.random.seed(42)
mc_x = np.random.uniform(0, 1, n_points)
mc_y = np.random.uniform(0, 1, n_points)

# Latin Hypercube Sampling
sampler = qmc.LatinHypercube(d=2, seed=42)
lhs_samples_2d = sampler.random(n=n_points)

Out[25]:

Visualization

Scatter plot of 500 random Monte Carlo points in the unit square showing irregular clustering and gaps. — Standard Monte Carlo sampling of 500 points over the unit square [0,1] × [0,1]. Random clustering leaves some regions overrepresented while others have sparse coverage, introducing variance in multi-dimensional integration estimates.

Scatter plot of 500 Latin Hypercube Sampling points in the unit square showing uniform space-filling coverage. — Latin Hypercube Sampling of 500 points over the same unit square. Each row and column interval contains exactly one sample, producing a more uniform space-filling distribution that eliminates variance caused by accidental clustering.

Notice how LHS provides much more uniform coverage of the unit square. This uniformity translates directly into reduced variance for estimating integrals over this space.