Quasi-Monte Carlo. Low-Discrepancy Sequences for Faster Convergence

Plain Monte Carlo converges at $O(1/\sqrt{N})$ no matter how cleverly you sample. Halving the error takes four times the samples, which gets expensive when you need accurate prices.

But that rate has a second half people skip: it is independent of dimension. A hundred dimensions converge at the same $1/\sqrt{N}$ as one. That is MC’s guarantee and equally its ceiling, because a method that knows nothing about the structure of the integrand cannot exploit any structure either.

Quasi-Monte Carlo (QMC) swaps pseudo-random samples for deterministic low-discrepancy sequences, and the error for a $d$-dimensional integral becomes $O((\log N)^d / N)$. Note where $d$ sits: in the exponent. This rate does depend on dimension. In low dimensions it leaves MC far behind; in high dimensions it collapses first.

So every practical question about QMC, whether to use it, how, and when not to, reduces to the same one: what is your $d$? This article keeps coming back to that.

Discrepancy, and Why It Controls the Error

A sequence of points $x_1, x_2, \ldots, x_N$ in $[0, 1]^d$ has star discrepancy

where the sup runs over all axis-aligned boxes anchored at the origin, $B = [0, u_1] \times \cdots \times [0, u_d]$. It measures how uniformly the points fill the space: low discrepancy means well distributed, high discrepancy means clustered or full of gaps.

What connects it to integration error is the Koksma-Hlawka inequality:

The right side splits cleanly in two. $V_{\text{HK}}(f)$ is the variation of $f$ in the sense of Hardy and Krause, a property of the integrand alone, and it is bounded for functions smooth enough, which covers most option payoffs. $D_N^*$ is a property of the point set alone and knows nothing about $f$. Which makes the route to better accuracy explicit: use points with lower discrepancy.

The two families are far apart. Pseudo-random points give $D_N^* = O(\sqrt{\log\log N / N})$, which substituted back yields the familiar $O(1/\sqrt{N})$. Low-discrepancy sequences give $D_N^* = O((\log N)^d / N)$, which is much better.

Read the bound with two caveats attached. First, it is often loose: actual QMC error is usually far smaller than the bound suggests, which makes it good for understanding when QMC helps and useless as an error estimate. Second, the $d$ in $(\log N)^d$ is real. Let $d$ grow and that factor eats the $1/N$ advantage. Where the crossover falls depends on $d$ and the smoothness of the integrand, and in moderate dimensions it happens at quite small $N$.

Sobol Sequences

The Sobol sequence is the most widely used low-discrepancy sequence in finance, a $(t, s)$-sequence in base 2 that gets its low discrepancy from a carefully chosen set of direction numbers.

Four properties explain its standing. It is deterministic: given the direction numbers, the sequence is fully determined, with no randomness. It is low-discrepancy: $D_N^* = O((\log N)^d / N)$ for fixed $d$. It is fast to generate, being a linear recurrence in base 2 that stays cheap at large $N$. And it tolerates dimension reasonably well, staying usable up to roughly 20-50.

The construction uses direction numbers $v_k$ at each bit position, themselves binary fractions. The $i$-th point is

where $i_k$ is the $k$-th bit of $i$ and $\oplus$ is bitwise XOR. The direction numbers are chosen so the points fill the unit cube evenly, and they are also the one part of Sobol that genuinely needs careful implementation: get one entry of the table wrong and the low-discrepancy property is gone. Most numerical libraries (scipy.stats.qmc.Sobol in Python) get this right, so unless you have a specific reason, do not write your own.

In use you specify $d$ and $N$ and the library hands back an $N \times d$ matrix in $[0, 1]^d$. One detail people skip: skip the first points. The opening stretch of a Sobol sequence has poor discrepancy, so discard a few hundred to a thousand points before using it. And since Sobol’s balance properties are defined on powers of two, taking $N = 2^m$ gets you the uniformity it actually promises, where 10,000 does not.

Halton, Faure, and the Rest

The Halton sequence generalizes van der Corput to multiple dimensions: dimension $j$ uses base $p_j$, the $j$-th prime (2, 3, 5, 7, 11, …). It is simple enough to need nothing but base conversion, no direction numbers, and its discrepancy is also $O((\log N)^d / N)$. As Sobol’s historical predecessor it remains useful for low-dimensional problems, roughly 1-5 dimensions, and it is the easiest to implement from scratch.

Halton’s problem is high dimensions. The primes get large and visible correlations appear between the early dimensions. Past $d > 10$, Sobol typically wins.

Faure sequences use one prime base $p \geq d$ for every dimension, forming a $(t, s)$-sequence with $t = 0$ in base $p$, with discrepancy comparable to Sobol. Their advantage is easier construction, since they need no primitive polynomials over $\mathbb{F}_2$. They are less common in finance and mostly serve as a fallback when you have a specific base requirement.

Two other names worth knowing: Niederreiter sequences generalize Sobol to arbitrary bases, more flexible and more complex; Latin hypercube is not strictly low-discrepancy but achieves comparable variance reduction.

As for randomized QMC, adding a random scramble to a low-discrepancy sequence, it is not a fourth entry on this list. It is a required accessory to any of them, for reasons below.

Why Finance Integrals Are Unusually Suited to This

If the QMC rate depends on dimension, and finance problems routinely have nominal dimensions in the dozens or hundreds, QMC should have been ruled out already. It is not, because the nominal dimension is not the $d$ that goes in the exponent.

Most finance integrands are flat in most dimensions, with the variation concentrated in a few important directions. For a European option price, the integrand is a smooth function of the Brownian increments, depends strongly on the terminal value, and barely notices the specific ordering of the increments. For a path-dependent option, the integrand nominally depends on every time step but mostly depends on the first few principal components of the path.

That is the effective dimension: the number of dimensions that actually contribute to the variance. An Asian option with 100 time steps has a nominal dimension of 100, but most of the variation lives in the first 3-5 principal components, so the effective dimension is 3-5. The QMC convergence rate runs on the effective dimension, $O((\log N)^{d_{\text{eff}}} / N)$, not the nominal one.

Typical magnitudes in finance: 1-3 for European options, 3-7 for path-dependent ones (Asian, lookback, barrier), and 5-20 for multi-asset options. Every one of those ranges sits comfortably inside the region where QMC beats MC. That is not a coincidence but a structural property of the problem class: it has a low-dimensional dominant subspace, and exploiting that subspace is precisely how QMC earns its keep.

Which also means reducing the effective dimension is itself a technique. Identify the highest-variance directions and sample them first, which is principal component analysis, or use the more common construction below.

The Brownian Bridge Brings $d$ Down

The standard Brownian construction generates increments in time order: $W_{t_1}$, $W_{t_2} - W_{t_1}$, through $W_{t_n} - W_{t_{n-1}}$. QMC fills the coordinates with those increments in sequence, so the effective dimension is the number of time steps. Every step is equally important and no coordinate especially deserves the low-discrepancy treatment.

The Brownian bridge construction reorders the work: sample the terminal value $W_T$ first, then the midpoint $W_{T/2}$ given $W_0$ and $W_T$, then the quarter points, recursively. The result is QMC samples ordered by importance, with the terminal value, the most important quantity, landing in the first coordinate, which is exactly where a Sobol sequence is most uniform.

The effect is to drop the effective dimension from $n$ to something on the order of $\log_2 n$. With $n = 100$ time steps, that is 100 down to single digits. It is an empirical regularity rather than a theorem, but for path-dependent options the Brownian bridge is among the most effective variance reduction techniques available, and it layers on top of QMC for a further improvement.

Scrambling: The Cost of Determinism, and the Cure

QMC has an immediate problem. The samples are deterministic, so you cannot estimate the error from a single run, and without an error estimate you don’t know when to stop.

Randomized QMC fixes it. Add a random shift or scramble that preserves the low-discrepancy property on average, run several independent scrambles, and use the standard deviation across runs as the error estimate. The estimate is unbiased and the convergence rate survives.

The most popular version is Owen scrambling (1995), which randomizes each digit of each point independently: generate Sobol points in base $b$, add a uniformly random digit from $\{0, 1, \ldots, b-1\}$ at each digit position, and the result is a scrambled Sobol point. It is what scramble=True gives you in scipy.stats.qmc.Sobol, and it is the standard randomization for finance QMC.

In practice run 10-100 scrambles, each with its own random seed, taking the mean as the price estimate and the standard deviation as the error:

import numpy as np
from scipy.stats import norm, qmc

def qmc_call_price(S0, K, r, sigma, T, n_paths=2**13, n_runs=20, seed=42):
    """European call price via randomized QMC (Owen-scrambled Sobol)."""
    rng = np.random.default_rng(seed)
    prices = []
    for _ in range(n_runs):
        # A fresh scramble per run: the runs must be independent randomizations
        # of the same sequence, not consecutive blocks of one sequence.
        sampler = qmc.Sobol(d=1, scramble=True, seed=int(rng.integers(2**32)))
        u = sampler.random(n_paths).ravel()
        # Apply inverse normal CDF
        Z = norm.ppf(u)
        ST = S0 * np.exp((r - 0.5*sigma**2)*T + sigma*np.sqrt(T)*Z)
        payoff = np.maximum(ST - K, 0)
        prices.append(np.exp(-r*T) * payoff.mean())
    return np.mean(prices), np.std(prices) / np.sqrt(n_runs)

mean, se = qmc_call_price(100, 100, 0.05, 0.20, 1.0)
print(f"QMC price: {mean:.4f} ± {se:.4f}")

One line there is easy to get wrong and worth calling out. Each run must construct a new sampler with a different scramble seed. Reuse one sampler and call random() repeatedly and you get consecutive blocks of a single sequence rather than independent randomizations, at which point the standard deviation across runs is not an error estimate but a meaningless number.

For calibration: under these parameters the standard deviation of the discounted payoff is about 14.7, so naive MC with 8,192 paths has a standard error near 0.16. The 20-scramble, 8,192-point QMC above lands in the $10^{-4}$ range, two to three orders of magnitude better. That gap is typical for a one-dimensional integrand with a single kink. Raise the dimension or roughen the payoff and it narrows.

A Multi-Asset Example

Take a basket call on $d$ correlated assets. Naive MC scatters $N$ points through $d$-dimensional space; at $d = 10$ the effective dimension after PCA is roughly 5-7, comfortably inside QMC territory.

def qmc_basket_call(S0_vec, K, r, sigma_vec, corr, T, n_paths=2**13, n_runs=10, seed=42):
    """Basket call via scrambled QMC with a control variate."""
    d = len(S0_vec)
    rng = np.random.default_rng(seed)
    L = np.linalg.cholesky(corr)
    expected_avg = S0_vec.mean() * np.exp(r*T)   # E[basket] under Q
    prices = []
    for _ in range(n_runs):
        sampler = qmc.Sobol(d=d, scramble=True, seed=int(rng.integers(2**32)))
        u = sampler.random(n_paths)
        Z = norm.ppf(u)
        # Apply correlation
        Z_corr = Z @ L.T
        ST = S0_vec * np.exp((r - 0.5*sigma_vec**2)*T + sigma_vec*np.sqrt(T)*Z_corr)
        basket = ST.mean(axis=1)
        payoff = np.maximum(basket - K, 0)
        # Control variate: the basket level itself, whose mean is known
        beta = np.cov(payoff, basket)[0, 1] / np.var(basket)
        cv_payoff = payoff - beta * (basket - expected_avg)
        prices.append(np.exp(-r*T) * cv_payoff.mean())
    return np.mean(prices), np.std(prices) / np.sqrt(n_runs)

Three things are stacked here: QMC sampling, Cholesky correlation, and a control variate. Total variance reduction over naive MC typically runs 100-1000x.

How QMC Relates to Other Variance Reduction

QMC counts as variance reduction, but it works on a different part of the machine than antithetic variates, control variates, and importance sampling. Those three change the estimator and leave the samples random. QMC changes the samples and leaves the estimator untouched.

Because they act on different parts, they compose, and the effect is multiplicative. A well-tuned combination of QMC plus control variates plus antithetic variates routinely gives three-figure variance reduction over naive MC.

The standard production recipe just bolts these on in order: Sobol for the base samples; a control variate with a known expectation; antithetic variates where appropriate; and a Brownian bridge to squeeze the effective dimension for path-dependent payoffs. None of the four steps is hard, and with all four in place the accuracy your compute budget buys is not in the same league as naive MC.

When It Doesn’t Work

Back to that $d$ in the exponent. QMC is not a silver bullet, and every one of its failure modes can be read straight off the rate.

When the effective dimension is too high. $(\log N)^{d}$ grows with $d$ while MC’s $1/\sqrt{N}$ does not move at all. Hold $N$ fixed, let $d$ grow, and QMC’s error climbs until it meets MC. Where they cross depends on $N$ and the smoothness of the integrand: for a typical finance problem with $N = 10{,}000$ and effective dimension 10, QMC is still far ahead; at 50 the two may be comparable or reversed; past 100, MC is usually better. The working boundary is around 30. Below it use QMC, above it use MC or switch to importance sampling or Latin hypercube. Which means the first step is always to know your effective dimension, not to install Sobol.

When the integrand is not smooth enough. $V_{\text{HK}}(f)$ in Koksma-Hlawka blows up for discontinuous functions. Digital options, and barrier options right at the boundary, see their QMC advantage shrink substantially. Smoothing techniques recover part of it.

When the sample size is small. QMC’s edge shows up at large $N$. Below a thousand points it is comparable to MC; above ten thousand the gap is wide.

And that error estimate you have to supply yourself. A deterministic sequence carries no measure of uncertainty. Either you randomize or you chase analytical bounds, and the analytical bounds are loose. In practice randomization is the only road.

For most finance MC problems none of these bind: the integrands are smooth, the effective dimension is moderate, and the libraries are well-tested.

Randomized QMC Theory, in One Line

Koksma-Hlawka is a deterministic worst-case bound, too loose to serve as an error estimate. Randomized QMC has a result closer to observed behavior: the expected error satisfies

where $c(d)$ depends on the smoothness of $f$ and the dimension $d$, and is small for smooth $f$ in moderate $d$. The theoretical foundation is in Owen (1995, 1998).

The practitioner’s takeaway is one sentence: use Owen scrambling, run 10-100 scrambles, average, report the standard deviation.

Where It Sits Among Integration Methods

Step back and QMC is a deterministic method for high-dimensional integration, wedged between two families. Grid-based methods (Gaussian quadrature, tensor-product grids) are accurate and limited by the curse of dimensionality. Monte Carlo ignores dimension and converges slowly.

The division of labor: at 1-3 dimensions, Gaussian quadrature, tensor-product grids, or solving the PDE outright all beat QMC; at 4-20, QMC is the best choice; above 20, Monte Carlo is often the only practical option. Solving the PDE is hard to beat on one to three underlyings, but it does not scale with dimension while QMC does, which is how high-dimensional problems end up with QMC by default.

Finance usually lands squarely in the middle band. It is also one of the few places where a sophisticated numerical technique has a clear practical edge over the simple alternative, since most of the time the gain from sophistication does not cover the bugs it introduces.

Parallelism and Production

QMC parallelizes readily: each scramble is independent, so you can run thousands at once on a cluster.

Do the arithmetic carefully, though. Say one machine runs 100 scrambles per second at $10^4$ points each, so 100 machines run $10^4$ scrambles per second for $10^8$ sample points total. The error reduction from averaging across scrambles goes as $1/\sqrt{K}$ with $K$ the number of scrambles, so 100 machines cut the standard error to a tenth of one machine’s, not to a ten-thousandth. That larger number comes from conflating variance reduction with sample count. The real gain is a product of two things: within each scramble you get the QMC rate over MC, and across scrambles you get $1/\sqrt{K}$.

Even so, for pricing that genuinely needs precision, a few seconds on commodity hardware gets you there. The cost is low and the benefit is real.

Why Most Quants Still Use Naive MC

The industry’s reasons for reaching for QMC go past speed and accuracy. Better prices for the same compute budget matter directly on tightly quoted instruments. Without randomization QMC is deterministic, so the same code with the same parameters returns the same result, which makes backtesting and validation easier (with randomization, record the seed and you get the same property). And it is indifferent to model, payoff, and underlying, as long as the effective dimension stays under roughly 20.

Given all that, adoption is still thin, and the reasons are enumerable. QMC is not taught in most finance courses, so plenty of people have never met it. Naive MC is one line of code, while QMC wants Sobol sequences, scrambling, and a Brownian bridge, which is genuinely more engineering. Randomization then demands multiple runs and error estimation, making the code longer still. Not every numerical library has decent QMC support (scipy.stats.qmc only does in recent versions). And there is cultural inertia, since “we have always done MC” is a powerful argument inside an established organization.

None of those reasons is technical. If you are starting a new MC project, use QMC from the beginning: the benefits are large, the cost is moderate, and library support is there. If you are maintaining an existing MC system, the migration usually pays.

For further reading, Niederreiter (1992), Random Number Generation and Quasi-Monte Carlo Methods, is the rigorous mathematical treatment; Glasserman (2004), Chapter 5, covers the finance-specific applications; Owen (1995, 1998) are the standard references for randomized QMC; and Sobol (1967) is the original paper on the sequence itself.

And back to the question running through all of it: before you install Sobol, estimate your effective dimension. That number decides everything downstream.