Michael Brenndoerfer | Data & AI, Private Equity, Technology

Michael Brenndoerfer

Integral Calculus and Differential Equations › Numerical Integration Methods

Michael BrenndoerferOctober 28, 20256 min read

Quantitative FinanceAll chapters

Part 1 of 15Mathematical Foundations

Chapter 7Integral Calculus and Differential Equations › Numerical Integration Methods

Reading section 10 of 14 · Numerical Integration Methods · 6 minSection 10 of 14 · 6 min

Section 10 of 14 · 6 min read

Integral Calculus and Differential Equations › Numerical Integration MethodsOne section · Sections in this chapter10 / 14 ↓

Sections in this chapter14

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

Many complex, path-dependent, or high-dimensional pricing and risk problems require numerical integration or simulation because a convenient closed form is unavailable. The appropriate method depends on the payoff, state dimension, model dynamics, and accuracy target.

Quadrature MethodsLink Copied

Numerical integration approximates definite integrals using weighted sums of function values.

\int_a^b f(x) \, dx \approx \sum_{i=1}^{n} w_i f(x_i)

where:

$a, b$ : integration bounds
$f(x)$ : the function being integrated
$n$ : number of evaluation points
$x_i$ : nodes (evaluation points) within $[a, b]$
$w_i$ : weights determining the contribution of each function value

The idea is to sample the function at carefully chosen points and combine these samples with appropriate weights. Different methods choose different nodes $x_i$ and weights $w_i$ :

The displayed composite trapezoidal rule uses equally spaced nodes and piecewise-linear interpolation; supplied coordinates also allow nonuniform nodes.
Simpson's rule uses quadratic interpolation over pairs of subintervals.
Gaussian quadrature chooses nodes and weights so that an $n$ -node rule integrates polynomials through degree $2n-1$ exactly under its specified weight function.

In[34]:

Code

import numpy as np
from scipy import integrate


# Compare integration methods on the same truncated lognormal expectation
def option_density(s, S0=100, r=0.05, sigma=0.2, T=1):
    """Log-normal density for stock price under GBM."""
    if s <= 0:
        return 0
    mu = np.log(S0) + (r - 0.5 * sigma**2) * T
    var = sigma**2 * T
    return np.exp(-((np.log(s) - mu) ** 2) / (2 * var)) / (s * np.sqrt(2 * np.pi * var))


# All numerical methods integrate over [0, upper_bound]. The analytic benchmark
# is the untruncated GBM expectation, so reported errors also include the tiny
# omitted tail beyond the upper bound.
S0, r, T = 100, 0.05, 1
upper_bound = 500.0
analytical_expected = S0 * np.exp(r * T)


def expectation_integrand(s):
    return s * option_density(s)


n_grid = 4001
grid = np.linspace(0, upper_bound, n_grid)
grid_values = np.array([expectation_integrand(s) for s in grid])
trapezoid_estimate = np.trapezoid(grid_values, grid)
simpson_estimate = integrate.simpson(grid_values, x=grid)

n_gauss = 64
gauss_nodes, gauss_weights = np.polynomial.legendre.leggauss(n_gauss)
mapped_nodes = 0.5 * upper_bound * (gauss_nodes + 1)
gauss_estimate = (
    0.5
    * upper_bound
    * np.sum(gauss_weights * np.array([expectation_integrand(s) for s in mapped_nodes]))
)

adaptive_estimate, adaptive_error, adaptive_info = integrate.quad(
    expectation_integrand, 0, upper_bound, full_output=True
)

quadrature_results = [
    ("Trapezoid", trapezoid_estimate, n_grid),
    ("Simpson", simpson_estimate, n_grid),
    ("Gauss-Legendre", gauss_estimate, n_gauss),
    ("Adaptive quad", adaptive_estimate, adaptive_info["neval"]),
]

Out[35]:

Console

Expected stock price integrated on [0, 500]
Analytic untruncated benchmark: $105.12710964

Method                   Estimate     Abs. error  Evaluations
--------------------------------------------------------------
Trapezoid          $105.12710964 $   6.395e-13         4001
Simpson            $105.12710964 $   6.537e-13         4001
Gauss-Legendre     $105.12710964 $   6.963e-13           64
Adaptive quad      $105.12710964 $   6.679e-13          147

Adaptive quad error estimate on the truncated integral: 2.66e-07

Monte Carlo IntegrationLink Copied

For high-dimensional integrals involving multiple assets or path-dependent options, Monte Carlo integration is often a practical and scalable approach.

\int f(x) p(x) \, dx \approx \frac{1}{N} \sum_{i=1}^{N} f(X_i)

where:

$f(x)$ : the function being integrated (e.g., option payoff)
$p(x)$ : probability density function
$N$ : number of random samples
$X_i$ : independent random samples drawn from distribution $p(x)$
$\frac{1}{N} \sum_{i=1}^{N} f(X_i)$ : sample average approximating the expected value

Under the Law of Large Numbers conditions, the sample average converges to the expected value. Drawing from $p(x)$ represents each region according to its probability mass. Accuracy depends on the variance of $f(X)$ , not density alone; for rare but influential tail payoffs, a different sampling density $q$ can reduce variance only with the corresponding likelihood-ratio weights $p(X)/q(X)$ and coverage wherever the original weighted payoff is nonzero. Changing the sampling density alone changes the target expectation.

In[36]:

Code

import numpy as np

# Monte Carlo integration for expected call payoff
S0, K, r, sigma, T = 100, 100, 0.05, 0.2, 1
n_simulations = 100000

# Generate log-normal terminal prices
Z = np.random.standard_normal(n_simulations)
ST = S0 * np.exp((r - 0.5 * sigma**2) * T + sigma * np.sqrt(T) * Z)

# Calculate call payoffs and present value
payoffs = np.maximum(ST - K, 0)
mc_price = np.exp(-r * T) * np.mean(payoffs)
mc_std_error = np.exp(-r * T) * np.std(payoffs, ddof=1) / np.sqrt(n_simulations)

Out[37]:

Console

European Call Option Pricing:
  Parameters: S0=$100, K=$100, r=5.0%, sigma=20.0%, T=1 year

Monte Carlo (100,000 simulations):
  Price: $10.4739
  Std Error: $0.0466

Black-Scholes Analytical:
  Price: $10.4506

Difference: $0.0233

Out[38]:

Visualization

Nested-sample Monte Carlo call-price estimates with pointwise normal-approximation intervals and the Black-Scholes benchmark. — Nested-sample Monte Carlo estimates for the displayed European call. For independent finite-variance payoffs, the estimator's standard error and asymptotic RMSE scale as 1/√N; the band is a pointwise normal-approximation interval.

For independent samples with finite payoff variance, the Monte Carlo sample mean has standard error $O(N^{-1/2})$ and, when unbiased, RMSE of the same order. Quadrupling $N$ therefore halves those asymptotic measures, not necessarily the realized absolute error in one run. Monte Carlo remains useful for complex and path-dependent payoffs without closed forms.