Numerical Methods and Algorithms in Finance › Numerical Integration

Numerical integration approximates definite integrals when analytical solutions are unavailable. In finance, this arises when computing expected values, pricing path-dependent options, or evaluating probability distributions. Integration is the continuous analog of summation. Many financial quantities are expressed as continuous sums over possible outcomes, weighted by their probabilities.

The connection between integration and expected values makes numerical integration useful in derivatives pricing. In a specified arbitrage-free model, let $N_t$ be the value of the chosen positive numeraire and $\mathbb{Q}^N$ its associated equivalent martingale measure. For a terminal payoff $H$ with $H/N_T$ integrable under that measure, the conditional-expectation price is $N_t\mathbb{E}^{\mathbb{Q}^N}[H/N_T\mid\mathcal{F}_t]$, where $\mathcal{F}_t$ represents information available at time $t$. Integrability concerns the numeraire-discounted payoff, not just $H$ itself. For an admissibly replicable claim, no-arbitrage fixes its replication price; in an incomplete market it need not select a unique value for a non-replicable claim. Quadrature and Monte Carlo are numerical routes to a specified expectation.

Consider the expected payoff of an option under a probability distribution $p(S_T)$:

$$\mathbb{E}[\text{Payoff}] = \int_0^{\infty} \text{Payoff}(S_T) \cdot p(S_T) \, dS_T$$

where:

• $\mathbb{E}[\text{Payoff}]$: expected value of the option payoff
• $S_T$: stock price at expiration time $T$
• $\text{Payoff}(S_T)$: option payoff as a function of terminal stock price
• $p(S_T)$: probability density function of the stock price at expiration

This integral accumulates the payoff at each possible terminal price, weighted by the probability density at that price. The integral runs from zero to infinity because the lognormal Black-Scholes model assigns stock prices positive, unbounded support. For simple European options under Black-Scholes, this integral has a closed form. For more complex payoffs or distributions, quadrature or simulation can approximate the expectation.

The trapezoidal rule approximates the area under a curve by dividing the interval into trapezoids. On each subinterval, it replaces the function with the straight line joining the two endpoint values and integrates that line exactly.

For an integral $\int_a^b f(x) \, dx$ with $n$ equal subintervals and $n+1$ equally spaced points $x_0 = a, x_1, ..., x_n = b$:

$$\int_a^b f(x) \, dx \approx \frac{h}{2}\left[f(x_0) + 2\sum_{i=1}^{n-1}f(x_i) + f(x_n)\right]$$

where:

• $h = \frac{b-a}{n}$: step size (width of each trapezoid)
• $f(x_0)$ and $f(x_n)$: function values at endpoints (each effectively weighted by $\frac{1}{2}$)
• $f(x_i)$ for $i = 1, \ldots, n-1$: function values at interior points (coefficient 2 appears because each interior point is shared by two adjacent trapezoids)

The factor of 2 for interior points arises because each interior point is shared by two adjacent trapezoids. The area of a single trapezoid with parallel sides $f(x_i)$ and $f(x_{i+1})$ is $\frac{h}{2}(f(x_i) + f(x_{i+1}))$; summing all trapezoids and collecting terms yields the formula above. The endpoints appear only once (each is the edge of only one trapezoid), while interior points appear twice (as the right edge of one trapezoid and the left edge of the next).

For a sufficiently smooth integrand with a bounded second derivative, the composite trapezoidal rule has asymptotic error $O(h^2)$. Once the grid is in that asymptotic regime, halving $h$ (approximately doubling the number of intervals) reduces the leading error by about a factor of four. This result follows from a Taylor expansion: the rule is exact for linear functions, while its leading error depends on the integrand's second derivative. Payoff kinks and other nonsmooth points require separate treatment, such as splitting the interval at the kink.

For sufficiently smooth integrands in the asymptotic grid regime, Simpson's rule can achieve higher accuracy than the trapezoidal rule by integrating a quadratic interpolant over each pair of panels. Payoff kinks can remove that advantage unless the interval is split at the nonsmooth point. Three values at distinct abscissas uniquely determine a polynomial of degree at most two, so we fit that interpolant over each pair of panels.

For $n$ intervals (where $n$ must be even):

$$\int_a^b f(x) \, dx \approx \frac{h}{3}\left[f(x_0) + 4\sum_{\substack{i=1,3,\ldots,n-1}}f(x_i) + 2\sum_{\substack{i=2,4,\ldots,n-2}}f(x_i) + f(x_n)\right]$$

where:

• Odd-indexed interior points ($i = 1, 3, \ldots, n-1$): weighted by 4 (midpoints of parabolic segments)
• Even-indexed interior points ($i = 2, 4, \ldots, n-2$): weighted by 2 (shared endpoints between segments); the sum is empty for $n=2$
• Endpoints $f(x_0)$ and $f(x_n)$: weighted by 1

The 1-4-2-4-2-...-4-1 weighting pattern comes from integrating the quadratic interpolant exactly, equivalently matching polynomial moments through degree 3 on each two-panel block. Odd-indexed points carry weight 4 because of this moment-matching calculation; even-indexed interior endpoints receive weight 2 because adjacent blocks share them.

To understand why Simpson's rule uses these specific weights, consider integrating a parabola $p(x) = ax^2 + bx + c$ over an interval $[-h, h]$ centered at the origin. The exact integral equals $\frac{h}{3}[p(-h) + 4p(0) + p(h)]$, which is precisely the Simpson's rule formula applied to three points. By choosing weights that make the formula exact for parabolas, we ensure high accuracy for any function that is approximately quadratic over each pair of sub-intervals.

For a sufficiently smooth integrand, composite Simpson's rule has error $O(h^4)$. In the asymptotic regime, halving the step size reduces that error by about a factor of 16, compared with about 4 for the trapezoidal rule. On a compatible even grid, both rules can reuse the same function evaluations with different weights.

Gauss-Legendre quadrature chooses both the evaluation points (called nodes or abscissas) and their weights rather than fixing equally spaced points.

With $n$ function evaluations, a quadrature rule has $n$ node positions and $n$ weights. Gauss-Legendre quadrature chooses both so that the rule is exact for every polynomial through degree $2n-1$. By contrast, an interpolatory rule on fixed generic nodes is guaranteed exact only through degree $n-1$, although symmetry can increase that degree.

For Gauss-Legendre quadrature on $[-1, 1]$:

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

where:

• $n$: number of quadrature points
• $x_i$: nodes (abscissas), roots of the $n$-th degree Legendre polynomial $P_n(x)$
• $w_i$: weights, computed to make the formula exact for polynomials up to degree $2n-1$

Example values:

• For $n = 2$: nodes at $x = \pm 1/\sqrt{3} \approx \pm 0.577$, weights $w = 1$
• For $n = 3$: nodes at $x = 0, \pm\sqrt{3/5} \approx \pm 0.775$, weights $w = 8/9, 5/9, 5/9$

The two- and three-node examples above are equally spaced. For $n \geq 4$, the nodes are not equally spaced; as $n$ increases, they become denser toward the interval endpoints. Their placement is determined by Legendre orthogonality and the requirement that an $n$-point Gauss-Legendre rule be exact through degree $2n-1$. Describing the endpoint density as boundary-error control can be a useful heuristic, but it is not the derivation of the nodes.

An $n$-point Gauss-Legendre rule is exact for polynomials of degree up to $2n-1$. For a nonpolynomial integrand, that degree of exactness does not by itself specify the approximation error; high accuracy follows when the integrand is sufficiently well approximated by polynomials on the interval.

Integration arises throughout quantitative finance in several key areas:

Under a specified equivalent martingale measure $\mathbb{Q}$ for the money-market numeraire and a constant deterministic short rate, the assigned price of a terminal payoff that is integrable under $\mathbb{Q}$ can be written as a discounted expectation:

$$V_0 = e^{-rT} \mathbb{E}^{\mathbb{Q}}[\text{Payoff}]$$

where:

• $V_0$: current value of the derivative
• $e^{-rT}$: discount factor from expiration to today
• $\mathbb{E}^{\mathbb{Q}}$: expectation under the risk-neutral probability measure
• $\text{Payoff}$: derivative's payoff at expiration, typically a function of the underlying asset price $S_T$

For a replicable claim, every admissible pricing measure agrees with its replication price, so no-arbitrage makes this value unique. In an incomplete market, different admissible measures can assign different expectation prices to a non-replicable payoff; the pricing model or measure must then be specified.

With stochastic rates or intermediate cash flows, the corresponding expression uses a stochastic discount factor or a chosen numeraire and conditional expectations. When the resulting expectations lack closed forms, numerical integration is one possible computational method.

Value at Risk (VaR) is a quantile of a loss distribution:

$$\text{VaR}_\alpha = \inf\{x : P(L > x) \leq 1 - \alpha\}$$

where:

• $\text{VaR}_\alpha$: Value at Risk at confidence level $\alpha$
• $L$: portfolio loss (positive values = losses)
• $\alpha$: confidence level (typically 95% or 99%)
• $\inf$: infimum (smallest value satisfying the condition)

This definition gives the lower $\alpha$-quantile of the loss distribution for $0<\alpha<1$. The loss distribution must correspond to a specified time horizon; VaR is not the maximum possible loss or an expected tail loss.

For a specified parametric distribution, a model quantile can be evaluated with its inverse CDF; for example, scipy.stats.norm.ppf evaluates a normal quantile. For historical or simulated loss samples, numpy.quantile computes a sample quantile directly from the supplied values rather than by numerically integrating a density.

Probability transforms in copula models and risk aggregation can involve integrating joint densities to compute marginal distributions or dependence measures. Analytic marginal distributions and sample-based methods are alternatives that need not perform numerical density integration.