Numerical Integration

Motivation: why this matters in quant finance

Most integrals in quantitative finance do not have closed-form solutions. The Black-Scholes formula is a celebrated exception — the log-normal assumption produces an integral that can be evaluated analytically. But the moment you move beyond vanilla European options on a single log-normal underlying, closed forms disappear: stochastic volatility models (Heston), local volatility models, path-dependent options (Asian, lookback, barrier), multi-asset options (baskets, rainbows), and American options all require numerical integration.

Numerical integration — also called quadrature — approximates

\int_a^b f(x)\,dx

by a weighted sum of function evaluations. The choice of method depends on the dimension of the integral, the smoothness of the integrand, and the accuracy required:

Low-dimensional, smooth integrands (1D–3D): deterministic quadrature (trapezoidal, Simpson's, Gaussian) is fast and accurate.
High-dimensional integrands (5D+): Monte Carlo simulation scales gracefully with dimension and is the workhorse of derivatives pricing.
Special structures (oscillatory integrands, semi-infinite domains): Fourier-based methods (FFT) exploit the characteristic function directly.

This page covers the main methods, their error properties, and when to use each.

Deterministic quadrature

General framework

A quadrature rule approximates $\int_a^b f(x)\,dx$ by:

Q_n[f] = \sum_{i=0}^{n} w_i\,f(x_i)

where

x_0, \ldots, x_n

are nodes (the points where

f

is evaluated) and

w_0, \ldots, w_n

are weights. Different choices of nodes and weights give different methods with different accuracy.

Trapezoidal rule

Approximate $f$ by a piecewise linear function on each subinterval. For $n$ equal subintervals of width $h = (b-a)/n$ :

T_n = h\left[\frac{f(a)}{2} + \sum_{i=1}^{n-1} f(a + ih) + \frac{f(b)}{2}\right]

Error:

|T_n - I| = O(h^2) = O(1/n^2)

. Doubling the number of points quadruples the accuracy.

When to use: Simple to implement, good for smooth functions, and surprisingly effective for periodic functions on a full period (where the error drops to

O(e^{-cn})

— exponential convergence).

Simpson's rule

Approximate $f$ by a piecewise quadratic (parabola) on each pair of subintervals. For $n$ subintervals ( $n$ even):

S_n = \frac{h}{3}\left[f(a) + 4\sum_{i \text{ odd}} f(a + ih) + 2\sum_{i \text{ even}, i \neq 0, n} f(a + ih) + f(b)\right]

Error:

|S_n - I| = O(h^4) = O(1/n^4)

. Much faster convergence than the trapezoidal rule for smooth functions.

When to use: The default choice for low-dimensional smooth integrals. Most implementations of numerical CDF evaluation (e.g., computing

\Phi(z)

for the normal distribution) use Simpson's rule or a variant internally.

Gaussian quadrature

Choose both the nodes

x_i

and the weights

w_i

optimally. Gauss-Legendre quadrature uses

n+1

nodes to integrate polynomials of degree up to

2n+1

exactly:

\int_{-1}^{1} f(x)\,dx \approx \sum_{i=0}^{n} w_i\,f(x_i)

where $x_i$ are the roots of the Legendre polynomial $P_{n+1}(x)$ and the weights are computed from these roots.

Error:

O(1/n^{2n+2})

for smooth functions — much faster than Simpson's rule.

Variants for specific domains and weight functions:

Variant	Domain	Weight function	Use case
Gauss-Legendre	$[-1, 1]$	1	Standard bounded integrals
Gauss-Hermite	$(-\infty, \infty)$	$e^{-x^2}$	Expectations under the normal distribution
Gauss-Laguerre	$[0, \infty)$	$e^{-x}$	Discounted payoff integrals

Gauss-Hermite is particularly relevant in quant finance: it is tailor-made for computing

\mathbb{E}[f(Z)]

where

Z \sim \mathcal{N}(0,1)

, because the weight function matches the Gaussian density. With only 10–20 nodes, it can achieve machine precision for smooth payoffs.

When to use: When high accuracy is needed for low-dimensional integrals and the integrand is smooth. Gauss-Hermite is the method of choice for semi-analytical pricing formulas that reduce to 1D or 2D integrals over normal variables.

Monte Carlo integration

The basic idea

Monte Carlo integration estimates $\int_a^b f(x)\,dx$ by random sampling. If $X \sim \text{Uniform}(a, b)$ , then:

\int_a^b f(x)\,dx = (b-a)\,\mathbb{E}[f(X)] \approx (b-a) \cdot \frac{1}{N}\sum_{i=1}^{N} f(X_i)

where $X_1, \ldots, X_N$ are i.i.d. draws. The law of large numbers guarantees convergence as $N \to \infty$ .

More generally, for any probability density $p(x)$ :

\mathbb{E}[f(X)] = \int f(x)\,p(x)\,dx \approx \frac{1}{N}\sum_{i=1}^{N} f(X_i), \qquad X_i \sim p

Error and convergence

The Monte Carlo estimator has standard error:

\text{SE} = \frac{\sigma_f}{\sqrt{N}}

where $\sigma_f^2 = \text{Var}(f(X))$ . The convergence rate is $O(1/\sqrt{N})$ — independent of the dimension of the integral.

This is the key advantage: For a 1D integral, deterministic quadrature converges at

O(1/n^4)

(Simpson's) or better, vastly outperforming Monte Carlo. But for a 100-dimensional integral (e.g., pricing a basket option on 100 assets), the convergence of deterministic methods degrades exponentially with dimension (the "curse of dimensionality"), while Monte Carlo remains at

O(1/\sqrt{N})

. In high dimensions, Monte Carlo wins.

The cost:

O(1/\sqrt{N})

convergence is slow. To halve the error, you need four times as many samples. For 2 digits of accuracy you need ~

10^4

samples; for 4 digits, ~

10^8

Pricing options by Monte Carlo

The risk-neutral pricing formula $V_0 = e^{-rT}\mathbb{E}^{\mathbb{Q}}[\text{payoff}(S_T)]$ is directly amenable to Monte Carlo:

Simulate $N$ independent paths of $S_T$ under $\mathbb{Q}$ :

S_T^{(i)} = S_0\exp\left(\left(r - \frac{\sigma^2}{2}\right)T + \sigma\sqrt{T}\,Z_i\right), \qquad Z_i \sim \mathcal{N}(0,1)

Compute the average discounted payoff:

\hat{V}_0 = e^{-rT}\frac{1}{N}\sum_{i=1}^{N}\text{payoff}(S_T^{(i)})

For a European call: $\text{payoff}(S_T^{(i)}) = (S_T^{(i)} - K)^+$ .

The standard error provides a confidence interval: $\hat{V}_0 \pm 1.96 \cdot \text{SE}$ is a 95% CI.

Variance reduction techniques

Since Monte Carlo convergence is $O(1/\sqrt{N})$ , reducing the variance $\sigma_f^2$ is as effective as increasing $N$ . The main techniques are:

Antithetic variates. For each sample

Z_i

, also use

-Z_i

. Since

Z

and

-Z

have the same distribution, both are valid samples, but they are negatively correlated, which reduces the variance of the average. For a call option on GBM, antithetic variates typically reduce variance by 50–80%.

Control variates. Use a correlated variable with known expectation to adjust the estimate. For example, the stock price itself has known expectation

\mathbb{E}^{\mathbb{Q}}[S_T] = S_0 e^{rT}

, so:

\hat{V}_{\text{CV}} = \hat{V}_0 - \beta\left(\frac{1}{N}\sum S_T^{(i)} - S_0 e^{rT}\right)

where $\beta$ is chosen to minimise variance. This exploits the correlation between the payoff and the control to cancel out noise.

Importance sampling. Sample from a distribution

q(x)

that concentrates probability where the integrand is large, and reweight by

p(x)/q(x)

. This is closely related to the change of measure framework: importance sampling is a numerical change of measure, and the optimal importance distribution is the one that minimises variance.

Stratified sampling. Divide the sampling space into strata and sample from each stratum separately, ensuring uniform coverage. Reduces variance by eliminating "clumping" of random samples.

Quasi-Monte Carlo

Replace the pseudo-random samples with low-discrepancy sequences (Sobol, Halton, Niederreiter) that fill the sample space more uniformly than random points. The convergence rate improves to approximately

O(1/N)

(up from

O(1/\sqrt{N})

) for sufficiently smooth integrands in moderate dimensions.

In practice, quasi-Monte Carlo with Sobol sequences is widely used in derivatives pricing desks, often providing 10–100x speedup over standard Monte Carlo for the same accuracy.

Fourier-based methods

The Carr-Madan FFT approach

For European options under models where the characteristic function of

\ln S_T

is known analytically (GBM, Heston, variance gamma, etc.), the option price can be written as:

C(K) = \frac{e^{-\alpha \ln K}}{\pi}\int_0^{\infty} e^{-iv\ln K}\,\psi(v)\,dv

where

\psi(v)

involves the characteristic function and

\alpha

is a damping parameter. This integral is a Fourier transform evaluated at

\ln K

, and it can be computed for many strikes simultaneously using the Fast Fourier Transform (FFT) in

O(N\log N)

operations.

Advantages: Extremely fast (prices at hundreds of strikes in milliseconds), exact up to discretisation error, and works for any model with a known characteristic function — including stochastic volatility and jump-diffusion models where the density itself has no closed form.

When to use: Calibrating models to the implied volatility surface, where you need option prices at many strikes and maturities. The FFT approach is the industry standard for Heston model calibration.

COS method

An alternative Fourier method that expands the density in a cosine series. Often faster and more accurate than FFT for individual option prices, with error control that is easier to analyse.

Choosing the right method

Scenario	Recommended method	Reason
1D smooth integral	Gaussian quadrature	Exponential convergence, machine precision with few nodes
$\Phi(z)$ evaluation	Simpson's or series expansion	Standard, fast, well-implemented in libraries
European option, known char. function	FFT / COS	$O(N\log N)$ , many strikes simultaneously
European option, no char. function	Monte Carlo with variance reduction	Flexible, dimension-independent
Path-dependent option (Asian, barrier)	Monte Carlo	Must simulate paths; no analytical shortcuts
High-dimensional (basket, multi-asset)	Monte Carlo / quasi-MC	Curse of dimensionality kills deterministic methods
American option	Least-squares Monte Carlo (Longstaff-Schwartz) or tree/PDE	Early exercise requires backward induction

Examples and applications

Example 1: computing $\Phi(z)$ by Simpson's rule

The standard normal CDF $\Phi(z) = \int_{-\infty}^{z} \phi(t)\,dt$ has no elementary antiderivative. In practice, it is computed by transforming to a finite interval and applying quadrature. For $z > 0$ , use the symmetry $\Phi(z) = 0.5 + \int_0^z \phi(t)\,dt$ and apply Simpson's rule to the finite integral $\int_0^z \phi(t)\,dt$ with $n = 20$ subintervals. Since $\phi$ is smooth and rapidly decaying, the error is negligible for all practical $z$ values.

Example 2: Monte Carlo pricing of an Asian call

An Asian call has payoff

(\bar{S} - K)^+

where

\bar{S} = \frac{1}{m}\sum_{j=1}^{m} S_{t_j}

is the average price over

m

observation dates. There is no closed-form price because the sum of log-normal random variables is not log-normal.

Monte Carlo approach:

For each path $i = 1, \ldots, N$ : simulate $S_{t_1}^{(i)}, \ldots, S_{t_m}^{(i)}$ under $\mathbb{Q}$ using the GBM recursion $S_{t_{j+1}} = S_{t_j}\exp((r - \sigma^2/2)\Delta t + \sigma\sqrt{\Delta t}\,Z_{ij})$ .
Compute the average $\bar{S}^{(i)} = \frac{1}{m}\sum_j S_{t_j}^{(i)}$ .
Compute the payoff $(\bar{S}^{(i)} - K)^+$ .
Estimate price $= e^{-rT}\frac{1}{N}\sum_i (\bar{S}^{(i)} - K)^+$ .

With $N = 100{,}000$ paths and $m = 252$ (daily observations), the standard error is typically a few cents on a $100 option — sufficient for trading purposes.

Example 3: Gauss-Hermite for a 1D Heston-style integral

After Fourier inversion, many stochastic volatility models reduce the option price to a 1D integral of the form:

C = \int_{-\infty}^{\infty} h(v)\,e^{-v^2/2}\,dv

where $h(v)$ involves the characteristic function. Gauss-Hermite quadrature with $n = 32$ nodes evaluates this to 12+ digits of accuracy, vastly outperforming Monte Carlo for this particular structure.

Common confusions and pitfalls

"Monte Carlo is always the best method." For low-dimensional smooth integrals, deterministic quadrature is faster by orders of magnitude. Monte Carlo's advantage is in high dimensions and for non-smooth or path-dependent integrands. Use the right tool for the problem.

"More random samples always improves accuracy." The standard error decreases as

1/\sqrt{N}

, so going from

10^4

10^5

samples only improves accuracy by a factor of

\sqrt{10} \approx 3.2

. Variance reduction (antithetic variates, control variates) is often more effective than brute-force sampling.

"Numerical integration is exact." All numerical methods have error. The question is whether the error is acceptable for the application. For trading decisions, a few cents of error in a Monte Carlo estimate may be fine. For model calibration, you may need the precision of Gaussian quadrature or FFT.

"The FFT works for all models." The FFT/COS methods require a known characteristic function. Models defined only by their SDE (e.g., local volatility with a general surface) do not have analytical characteristic functions and must use Monte Carlo or PDE methods instead.

Where this goes next

Numerical integration is the computational engine that makes the analytical framework practical:

The Riemann Integral: Numerical integration approximates the Riemann integral when no antiderivative exists.
Normal Distribution: The CDF $\Phi(z)$ is computed by numerical integration, typically Simpson's rule or a rational approximation.
Log-Normal Distribution: Option pricing integrals under the log-normal model are solved analytically (Black-Scholes) or numerically (Monte Carlo for exotics).
Change of Variables: The change-of-measure framework connects importance sampling (a Monte Carlo technique) to the Radon-Nikodym derivative and risk-neutral pricing.
Brownian Motion and GBM: Monte Carlo pricing simulates paths of these processes and averages discounted payoffs.

References

Stewart, J. (2008). Single Variable Calculus: Early Transcendentals (6th ed.). Thomson Brooks/Cole. Ch. 7 Section 7.7 (Approximate Integration) for midpoint, trapezoidal, Simpson, and error-estimate ideas.

Numerical Integration

Motivation: why this matters in quant finance

Deterministic quadrature

General framework

Trapezoidal rule

Simpson's rule

Gaussian quadrature

Monte Carlo integration

The basic idea

Error and convergence

Pricing options by Monte Carlo

Variance reduction techniques

Quasi-Monte Carlo

Fourier-based methods

The Carr-Madan FFT approach

COS method

Choosing the right method

Examples and applications

Example 1: computing Φ(z)\Phi(z)Φ(z) by Simpson's rule

Example 2: Monte Carlo pricing of an Asian call

Example 3: Gauss-Hermite for a 1D Heston-style integral

Common confusions and pitfalls

Where this goes next

References

Example 1: computing $\Phi(z)$ by Simpson's rule