Numerical Integration: From Simple Sums to Quasi-Random Sequences

Chapter 81 tackles a problem that looks simple on the surface but gets surprisingly deep: how do you calculate a multi-dimensional integral when you cannot do it with pen and paper? If you can write an option price as an integral (and for many European options on multiple assets, you can), then all you need is a good way to evaluate that integral numerically. Wilmott shows three approaches, and the last one is genuinely clever.

When Option Prices Are Integrals

For a non-path-dependent European option on multiple lognormal assets, the fair value can be written as a big integral. The payoff function gets multiplied by a multivariate normal probability density and integrated over all possible asset values at expiry, then discounted back. If you have one asset, it is a single integral. Two assets, a double integral. Five assets, a quintuple integral. And so on.

The question is: how do you evaluate a d-dimensional integral efficiently? Three methods, each with different strengths.

Method 1: The Regular Grid

The simplest approach. Lay down an evenly spaced grid in d dimensions. Evaluate the function at every grid point. Add everything up using something like the trapezium rule or the midpoint rule.

If you use N total grid points, you have N to the power of 1/d points in each direction. The error is O(N to the power of -2/d). For one dimension, this is O(N squared), which is great. For two dimensions, O(N inverse). For five or more dimensions, the error shrinks painfully slowly as you add more points.

This is the “curse of dimensionality.” In high dimensions, a regular grid needs an astronomical number of points to achieve decent accuracy. Also, the payoff function typically has a kink (like the hockey stick of a call option), so higher-order methods like Simpson’s rule do not help much unless you carefully locate the kink.

Method 2: Monte Carlo Integration

Instead of evaluating on a regular grid, just throw random points into the d-dimensional space. Evaluate the function at each point and take the average, multiplied by the volume.

The big advantage: the error is O(N to the power of -1/2), regardless of the number of dimensions. In one or two dimensions, the regular grid wins. But once you hit five or more dimensions, Monte Carlo beats the grid.

Wilmott gives a code fragment for pricing a European basket option on up to five lognormal assets using this method. The process is:

1. Generate uncorrelated normal random variables (using Box-Muller).
2. Transform them into correlated variables (using Cholesky factorization from Chapter 80).
3. Calculate the lognormal asset prices at expiry.
4. Evaluate the payoff and accumulate the sum.
5. Divide by the number of evaluations and discount.

Simple, general, and works for any number of assets. The catch is that the convergence is slow, and the random nature of the sampling means the points are not evenly distributed. Some areas get clumped together while others get left empty. This waste of sampling effort is what the next method addresses.

Method 3: Low-Discrepancy (Quasi-Random) Sequences

This is the star of the chapter. The idea starts with a simple observation about Monte Carlo: random points clump. Look at 200 random points in a unit square and you will see clusters and gaps everywhere. If you could spread the points more evenly, you would get a better estimate of the integral with the same number of points.

But you cannot just use a regular grid, because if you want to add more points later, the grid is no longer regular. You need a sequence that is well-spread at any number of points, and stays well-spread when you add more.

This is exactly what low-discrepancy sequences (also called quasi-random sequences) provide. Despite the name, there is nothing random about them. They are deterministic sequences designed to fill space as uniformly as possible.

The Halton Sequence

Wilmott focuses on the Halton sequence because it is the easiest to explain. Pick a base, say 2. Write the positive integers in binary: 1, 10, 11, 100, 101, 110, 111, and so on. Now “reflect” each number around the decimal point. The integer 1 (binary: 1) becomes 0.1 (which is 0.5). The integer 2 (binary: 10) becomes 0.01 (which is 0.25). The integer 3 (binary: 11) becomes 0.11 (which is 0.75).

The resulting sequence is: 0.5, 0.25, 0.75, 0.125, 0.625, 0.375, 0.875, …

Notice something beautiful. Each new number lands as far as possible from the numbers already placed. First you hit 0.5 (the middle). Then 0.25 and 0.75 (the quarters). Then the eighths. The sequence naturally fills in the gaps at finer and finer levels.

For multi-dimensional integration, use different prime bases for each dimension. In two dimensions, use base 2 for the x-coordinate and base 3 for the y-coordinate. Point i has coordinates (Halton(i, 2), Halton(i, 3)). The points spread beautifully across the unit square with no clumping.

Wilmott provides a short Visual Basic function to generate Halton numbers of any base. It is maybe 10 lines of code.

How Much Better Is It?

The error for quasi-random integration is O(log(N) to the power of d divided by N), which is much better than the O(N to the power of -1/2) of Monte Carlo. In practice, even when the theoretical advantage is unclear because of the log factors, quasi-random methods are at worst as good as Monte Carlo and often much better.

Wilmott shows a graph comparing the relative error for a five-dimensional contract using regular Monte Carlo versus a Halton sequence. The Halton results show a clear inverse relationship between error and number of points. The Monte Carlo results are noisier and converge more slowly.

There is also a nice bonus property: if you project the quasi-random points onto a lower dimension (collapse the 2D points onto a line, for example), the points do not pile up on top of each other. They stay well-distributed. This means that if your integrand depends strongly on only a few of the variables, the method still works great.

Other Sequences

Halton is the simplest, but not the best. Sobol sequences are generally considered superior for higher dimensions. The tradeoff is that Sobol is much harder to implement and explain. Wilmott points the reader to Press et al. (1995) for code.

Barrett, Moore, and Wilmott (1992) were the first to apply these techniques to finance, which Wilmott mentions with characteristic honesty, crediting his co-authors for the real numerical analysis work.

Advanced Variance Reduction

Wilmott briefly covers several techniques that improve convergence without changing the fundamental rate:

Antithetic variables (from Chapter 80) work just as well for integration as for simulation. Use each point and its “mirror.” Always do this because it costs almost nothing and helps.

Control variates work the same way too. If you have a similar integral with a known answer, compute both with the same points and use the known answer to calibrate your estimate.

Importance sampling transforms the integration variables so the integrand becomes as close to constant as possible. The more constant the integrand, the better any sampling method performs. Wilmott notes this is rarely used in finance.

Stratified sampling divides the integration region into subregions and allocates points optimally based on where the variance is highest. Recursive versions make the decision about subdividing each region automatically. Also rarely used in finance, partly because in high dimensions it starts looking like laying down a grid.

Why This Matters

For anyone pricing basket options, multi-asset derivatives, or anything with more than a few underlying variables, the choice of integration method matters enormously. A naive regular grid might need billions of points in 10 dimensions. Monte Carlo gets you there with millions. Quasi-random methods can get you there with hundreds of thousands.

The first application of low-discrepancy sequences to finance came from Wilmott’s own group in 1992, and by now the technique is standard on trading desks worldwide. It is one of those cases where a bit of number theory (reflecting integers around the decimal point in different bases) directly translates into faster, more accurate option prices.

Key Takeaways

1. European option prices on multiple assets can often be written as multi-dimensional integrals. The challenge is evaluating them efficiently.
2. Regular grids suffer from the curse of dimensionality. The error grows rapidly with the number of dimensions.
3. Monte Carlo integration has error O(N to the -1/2) regardless of dimension. It beats grids above about five dimensions, but it wastes effort because random points clump.
4. Low-discrepancy (quasi-random) sequences like Halton and Sobol fill space uniformly and converge faster than Monte Carlo. They are at worst as good as Monte Carlo and usually significantly better.
5. Variance reduction techniques (antithetic variables, control variates, importance sampling) reduce the constant in front of the error term. Antithetic variables are essentially free and should always be used.

Previous post: Monte Carlo Simulation: Pricing by Random Sampling

Next post: Code Programs: Finite Differences and Monte Carlo in Action