Numerical integration
First thing: integration is trickier than differentiation
. . .
We need to do integration for many things in economics
. . .
. . .
Integrals are effectively infinite sums
. . .
1-dimensional example:
\lim_{dx_i\rightarrow0}\sum_{i=0}^{(a-b)/dx_i} f(x_i) dx_i
where
Just like derivatives, we face an infinite limit because as {dx_i} \rightarrow 0, \frac{(a-b)}{dx_i} \rightarrow \infty
. . .
We avoid this issue in the same way as derivatives: we replace the infinite sum with something we can handle
. . .
We will see two types of methods:
It’s probably the most commonly used form in empirical economics
. . .
Key idea: approximate an integral by relying on LLN and “randomly” sampling the integration domain
. . .
We integrate \xi = \int_0^1 f(x) dx
by drawing N uniform samples, x_1,\dots,x_N over interval [0,1]
. . .
. . .
. . .
The variance of \hat{\xi} is \sigma^2_{\hat{\xi}} = Var(\frac{V}{N}\sum_{i=1}^N f(x_i)) = \frac{V^2}{N^2} \sum_{i=1}^N Var(f(X)) = \frac{V^2}{N}\sigma^2_{f(X)}
. . .
So average error is \frac{V}{\sqrt{N}}\sigma_{f(X)}. This gives us its rate of convergence: O(1/\sqrt{N})
Notes:
Suppose we want to integrate x^2 from 0 to 10, we know this is 10^3/3 = 333.333
# Package for drawing random numbers
using Distributions
# Define a function to do the integration for an arbitrary function
function integrate_mc(f, lower, upper, num_draws)
# Draw from a uniform distribution
xs = rand(Uniform(lower, upper), num_draws)
# Expectation = mean(x)*volume
expectation = mean(f(xs))*(upper - lower)
endintegrate_mc (generic function with 1 method)Suppose we want to integrate x^2 from 0 to 10, we know this is 10^3/3 = 333.333
f(x) = x.^2;
integrate_mc(f, 0, 10, 1000)320.7085212559472Quite close for a random process!
We can also approximate integrals using a technique called quadrature
. . .
With quadrature we effectively take weighted sums to approximate integrals
. . .
We will focus on two classes of quadrature for now:
Suppose we want to integrate a one dimensional function f(x) over [a,b]
How would you do it?
. . .
One answer is to replace the function with something easier to integrate: a piece-wise polynomial
. . .
Key things to define up front:
Most basic Newton-Cotes method:
. . .
\int_{x_i}^{x_{i+1}} f(x) dx \approx hf(\frac{1}{2}(x_{i+1}+x_i))
. . .
Approximates f by a step function
Let’s integrate x^2 again from 0 to 100 (the answer is 333.33…)
# Generic function to integrate with midpoint
function integrate_midpoint(f, a, b, N)
# Calculate h given the interval [a,b] and N nodes
h = (b-a)/(N-1)
# Calculate the nodes starting from a+h/2 til b-h/2
x = collect(a+h/2:h:b-h/2)
# Calculate the expectation
expectation = sum(h*f(x))
end;f(x) = x.^2;
println("Integrating with 5 nodes:$(integrate_midpoint(f, 0, 10, 5))")
println("Integrating with 50 nodes:$(integrate_midpoint(f, 0, 10, 50))")
println("Integrating with 100 nodes:$(integrate_midpoint(f, 0, 10, 100))")Integrating with 5 nodes:328.125
Integrating with 50 nodes:333.29862557267813
Integrating with 100 nodes:333.3248307995783. . .
Decent approximation with only 5 nodes. We get pretty close to the answer once we move up to 100 nodes.
Increase complexity by 1 degree:
. . .
\int_{x_i}^{x_{i+1}} f(x) dx \approx \frac{h}{2}[f(x_i) + f(x_{i+1})]
. . .
We can aggregate this up to
\int_{a}^{b} f(x) dx \approx \sum_{i=1}^n w_i f(x_i)
where w_1=w_n = h/2 and w_i = h otherwise
Trapezoid rule is O(h^2) and first-order exact: it can integrate any linear function exactly!
Increase complexity by 1 degree:
Let n be odd then approximate the function across a pair of subintervals by a quadratic interpolation passing through
. . .
\int_{x_{i-1}}^{x_{i+1}} f(x) dx \approx \frac{h}{3}[f(x_{i-1}) + 4f(x_{i}) + f(x_{i+1})]
\int_{x_{i-1}}^{x_{i+1}} f(x) dx \approx \frac{h}{3}[f(x_{i-1}) + 4f(x_{i}) + f(x_{i+1})]
We can aggregate this by summing over every two subintervals
\int_{a}^{b} f(x) dx \approx \sum_{i=1}^n w_i f(x_i)
where w_1=w_n = h/3, otherwise and w_i = 4h/3 if i is even and w_i = 2h/3 if i is odd
Simpson’s rule is O(h^4) and third-order exact: it can integrate any cubic function exactly
. . .
That’s weird! Why do we gain 2 orders of accuracy when increasing one order of approximation complexity?
. . .
. . .
. . .
The derivation of the order of error terms is similar to what we did for forward differences: we use Taylor expansions
But there are some extra steps because we have to bring the integral in. We’re going to skip these derivations
BUT it is important to know that these derivations assume smooth functions. That might not be the case in economic models because corner solutions might exist!
We can generalize the unidimensional Newton-Cotes quadrature with tensor products
. . .
Suppose we want to integrate over a rectangle
\{(x_1, x_2) | a_1 \leq x_1 \leq b_1, a_2 \leq x_2 \leq b_2\} \in \mathbb{R}^2
. . .
We compute separately the n_1 nodes and weights (x_{1i}, w_{1i}) and n_2 for (x_{2j}, w_{2j}). The tensor product will give us a grid with n = n_1 n_2 points
. . .
\int_{a_1}^{b_1} \int_{a_2}^{b_2} f(x_1, x_2) dx_2 dx_1 \approx \sum_{i=1}^{n_1} \sum_{j=1}^{n_2} w_{1i} w_{2j} f(x_{1i}, x_{2j})
. . .
We can use the same idea for 3+ dimensions. But the number of points increases very quickly: curse of dimensionality
How did we pick the x_i quadrature nodes for Newton-Cotes rules?
. . .
Evenly spaced, but no particular reason for doing so…
. . .
Gaussian quadrature selects these nodes more efficiently and relies on weight functions w(x)
Gaussian rules try to exactly integrate some finite dimensional collection of functions (i.e. polynomials up to some degree)
. . .
For a given order of approximation n, the weights w_1,...,w_n and nodes x_1,...,x_n are chosen to satisfy 2n moment matching conditions:
. . .
\int_I s^kw(s)ds = \sum_{i=1}^n w_i x^k_i, \;\text{for}\; k=0,...,2n-1
where I is the interval over which we are integrating and w(x) is a given weight function
The moment matching conditions pin down w_is and x_is so we can approximate an integral by a weighted sum of the function at the prescribed nodes
. . .
\int_I f(s) w(s)ds \approx \sum_{i=1}^n w_i f(x_i)
. . .
Gaussian rules are 2n-1 order exact: we can exactly compute the integral of any polynomial order 2n-1
Gaussian quadrature effectively discretizes some distribution p(x) into mass points (nodes) and probabilities (weights) for some other discrete distribution \bar{p}(x)
. . .
. . .
But what do we pick for the weighting function w(s)?
We can start out with a simple w(s) = 1: this gives us Gauss-Legendre quadrature
This can approximate the integral of any function arbitrarily well by increasing n
Sometimes we want to compute exponentially discounted sums like
\int_I f(s) e^{-s} ds
The weighting function e^{-s} is Gauss-Laguerre quadrature
. . .
Sometimes we want to take expectations of functions of normally distributed variables
\int_I f(s) e^{-s^2} ds
. . .
OK, Gaussian quadrature sounds nice. But how do we actually find the nodes and weights?
. . .
It’s a non-trivial task: we have to solve 2n nonlinear equations to get n pairs (x_i, w_n)
We’ll use a package for that: QuantEcon.jl
FastGaussQuadrature.jl does a better job with high-order integrals or weights with singularitiesQuadGK has methods to improve precision further with adaptive subdomainsQuantEcon.jl: Gaussian quadratureLet’s go back to our previous example and integrate x^2 from 0 to 10 again (the answer is 333.333…)
using QuantEcon;
# Our generic function to integrate with Gaussian quadrature
function integrate_gauss(f, a, b, N)
# This function get nodes and weights for Gauss-Legendre quadrature
x, w = qnwlege(N, a, b)
# Calculate expectation
expectation = sum(w .* f(x))
end;QuantEcon.jl: Gaussian quadraturef(x) = x.^2;
println("Integrating with 1 node:$(integrate_gauss(f, 0, 10, 1))")
println("Integrating with 2 nodes:$(integrate_gauss(f, 0, 10, 2))")Integrating with 1 node:250.0
Integrating with 2 nodes:333.3333333333334. . .
All we needed was 2 nodes!
QuantEcon.jl: Other quadraturesQuantEcon.jl has a lot more to offer1
qnwtrap) and Simpson’s (qnwsimp) ruleqnwnorm) or Gaussian with any continuous distribution (qnwdist)quadrect)QuantEcon.jlLet’s integrate f(x) = \frac{1}{x^2 + 1} from 0 to 1
Tip: Type ? in the REPL and the name of the function to read its documentation and learn how to call it
quadrect to integrate using Monte Carlof to accept an array (. syntax)qnwtrap to integrate using the Trapezoid rule quadratureqnwlege to integrate using Gaussian quadrature. . .
Do we know what is the analytic solution to this integral?
QuantEcon.jlf(x) = 1.0./(x.^2 .+ 1.0);
# 1. Use `quadrect` to integrate using Monte Carlo with 1000 nodes
quadrect(f, 1000, 0.0, 1.0, "R")0.7896227636186307# 2. Use `qnwtrap` to integrate using the Trapezoid rule quadrature with 7 nodes
x, w = qnwtrap(7, 0.0, 1.0);
# Perform quadrature approximation
sum(w .* f(x))0.7842407666178157# 3. Use `qnwlege` to integrate using Gaussian quadrature with 7 nodes
x, w = qnwlege(7, 0.0, 1.0);
# Perform quadrature approximation
sum(w .* f(x))0.785398164063432. . .
By the way, \int^1_0 1/(x^2 + 1) = \arctan(1) = \pi/4
π/40.7853981633974483Computing expected utility under quality uncertainty
We’ll compare Trapezoid, Gauss-Legendre, and Gauss-Hermite quadrature on a problem that arises naturally in economics: taking expectations over normally distributed shocks
A consumer has Cobb-Douglas utility over two goods x_1, x_2 with an uncertain quality shifter \theta on good 2:
u(x_1, x_2; \theta) = x_1^\alpha \cdot (\theta \, x_2)^{1-\alpha}
. . .
The quality \theta is log-normally distributed: \ln\theta \sim N(\mu, \sigma^2)
. . .
For fixed consumption bundle (x_1, x_2), the consumer’s expected utility is:
E[u] = x_1^\alpha \, x_2^{1-\alpha} \int_{-\infty}^{\infty} e^{(1-\alpha)z} \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(z-\mu)^2}{2\sigma^2}} dz
where z = \ln\theta. Can you compute this integral by hand?
This integral has a known closed-form solution. Since e^{(1-\alpha)z} with z \sim N(\mu, \sigma^2) is exactly the MGF of the normal evaluated at 1-\alpha, we have:
E[u] = x_1^\alpha \, x_2^{1-\alpha} \exp\left[(1-\alpha)\mu + \frac{(1-\alpha)^2 \sigma^2}{2}\right]
. . .
So we can check exactly how well each numerical method does
The integral we need is E[e^{(1-\alpha)z}] where z = \ln\theta \sim N(\mu, \sigma^2). This happens to be the moment generating function of the normal evaluated at t = (1-\alpha)
. . .
The MGF of z \sim N(\mu, \sigma^2) is M_z(t) = E[e^{tz}] = \exp\left(t\mu + \frac{t^2 \sigma^2}{2}\right)
. . .
Derivation: the product e^{tz} \cdot \phi(z) has combined exponent
tz - \frac{(z-\mu)^2}{2\sigma^2} = -\frac{1}{2\sigma^2}\left[z^2 - 2(\mu + t\sigma^2)z + \mu^2\right]
. . .
Complete the square by adding and subtracting (\mu + t\sigma^2)^2:
= -\frac{(z - (\mu + t\sigma^2))^2}{2\sigma^2} + t\mu + \frac{t^2\sigma^2}{2}
\int_{-\infty}^{\infty} e^{tz} \phi(z) dz = \exp\!\left(t\mu + \frac{t^2\sigma^2}{2}\right) \underbrace{\int_{-\infty}^{\infty} \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{(z-(\mu+t\sigma^2))^2}{2\sigma^2}} dz}_{= 1 \text{ (kernel of } N(\mu+t\sigma^2,\, \sigma^2)\text{)}}
. . .
So M_z(t) = \exp(t\mu + t^2\sigma^2/2). Substituting t = 1-\alpha:
E[u] = x_1^\alpha \, x_2^{1-\alpha} \exp\left[(1-\alpha)\mu + \frac{(1-\alpha)^2 \sigma^2}{2}\right]
Let’s set \alpha = 0.4, x_1 = 2, x_2 = 3, \mu = 0, \sigma = 0.5 and calculate the true value
using QuantEcon, Plots, Distributions, LaTeXStrings;
# Parameters
α = 0.4
x₁, x₂ = 2.0, 3.0
μ, σ = 0.0, 0.5
# The integrand: g(z) * ϕ(z) where g(z) = x₁^α * (exp(z)*x₂)^(1-α) and ϕ is Normal pdf
g(z) = x₁^α * (exp.(z) .* x₂).^(1-α)
# Analytic solution via MGF of normal distribution
EU_true = x₁^α * x₂^(1-α) * exp((1-α)*μ + 0.5*(1-α)^2*σ^2)
println("True expected utility: $(round(EU_true, digits=8))")True expected utility: 2.66825912The trapezoid rule needs a finite interval, so we truncate the normal at \mu \pm 4\sigma and integrate g(z)\phi(z) directly
function eu_trapezoid(n)
a, b = μ - 4σ, μ + 4σ # Truncated domain
x, w = qnwtrap(n, a, b)
# We must include the normal pdf in the integrand
ϕ = pdf.(Normal(μ, σ), x)
return sum(w .* g(x) .* ϕ)
end;
for n in [5, 11, 21, 51]
eu = eu_trapezoid(n)
err = abs(eu - EU_true)
println("Trapezoid n=$n: EU=$(round(eu, digits=8)), error=$(round(err, sigdigits=3))")
endTrapezoid n=5: EU=2.68958165, error=0.0213
Trapezoid n=11: EU=2.66772788, error=0.000531
Trapezoid n=21: EU=2.66788868, error=0.00037
Trapezoid n=51: EU=2.66793885, error=0.00032Same truncation strategy, but now nodes are not evenly spaced: they cluster where they matter most for polynomial exactness
function eu_legendre(n)
a, b = μ - 4σ, μ + 4σ
x, w = qnwlege(n, a, b)
ϕ = pdf.(Normal(μ, σ), x)
return sum(w .* g(x) .* ϕ)
end;
for n in [5, 11, 21, 51]
eu = eu_legendre(n)
err = abs(eu - EU_true)
println("Legendre n=$n: EU=$(round(eu, digits=8)), error=$(round(err, sigdigits=3))")
endLegendre n=5: EU=2.78605187, error=0.118
Legendre n=11: EU=2.66795567, error=0.000303
Legendre n=21: EU=2.6679487, error=0.00031
Legendre n=51: EU=2.6679487, error=0.00031Gauss-Hermite is designed for integrals of the form \int f(x) e^{-x^2} dx
The normal pdf is this form (after a change of variables), so qnwnorm absorbs the density into the weights: no need to multiply by \phi(z)
function eu_hermite(n)
# qnwnorm gives nodes/weights for N(μ, σ²) directly
x, w = qnwnorm(n, μ, σ^2)
return sum(w .* g(x))
end;
for n in [3, 5, 7, 11]
eu = eu_hermite(n)
err = abs(eu - EU_true)
println("Hermite n=$n: EU=$(round(eu, digits=8)), error=$(round(err, sigdigits=3))")
endHermite n=3: EU=2.6682433, error=1.58e-5
Hermite n=5: EU=2.66825912, error=5.09e-10
Hermite n=7: EU=2.66825912, error=1.02e-14
Hermite n=11: EU=2.66825912, error=2.22e-15. . .
Notice Gauss-Hermite needs far fewer nodes. The weight function e^{-x^2} is doing the heavy lifting
ns_trap = 3:2:51 # odd numbers for compatibility
ns_lege = 3:2:51
ns_herm = 3:2:11
errs_trap = [abs(eu_trapezoid(n) - EU_true) for n in ns_trap]
errs_lege = [abs(eu_legendre(n) - EU_true) for n in ns_lege]
errs_herm = [abs(eu_hermite(n) - EU_true) for n in ns_herm]
plot(ns_trap, errs_trap, label="Trapezoid", lw=2, marker=:circle, ms=3, yscale=:log10,
xlabel="Number of nodes", ylabel="Absolute error (log scale)",
title="Convergence: Expected utility of Cobb-Douglas",
legend=:bottomright, size=(700, 400), ylims=(1e-16, 1e0))
plot!(ns_lege, errs_lege, label="Gauss-Legendre", lw=2, marker=:square, ms=3)
plot!(ns_herm, errs_herm, label="Gauss-Hermite", lw=2, marker=:diamond, ms=3)
hline!([eps()], label="Machine epsilon", ls=:dash, color=:gray)The integral has the form \int g(z) \phi(z) dz where \phi is the normal pdf
. . .
. . .
Takeaway: when your integrand has a known density component, choose the quadrature whose weight function matches it
Let’s see where each method places its 5 nodes and how it weights them
# Get nodes and weights for n=5
x_t, w_t = qnwtrap(5, μ-4σ, μ+4σ)
x_l, w_l = qnwlege(5, μ-4σ, μ+4σ)
x_h, w_h = qnwnorm(5, μ, σ^2)
# For Trapezoid and Legendre, effective weight = w * ϕ(x) (contribution to final sum)
eff_t = w_t .* pdf.(Normal(μ, σ), x_t)
eff_l = w_l .* pdf.(Normal(μ, σ), x_l)
# For Hermite, the weights already incorporate the density
eff_h = w_h
p1 = bar(x_t, eff_t, bar_width=0.03, label="Trapezoid", alpha=0.8, color=:steelblue,
title="Effective weights (wᵢ × ϕ(xᵢ) or wᵢ for Hermite)", ylabel="Effective weight",
xlabel=L"z = \ln\theta", legend=:topright, size=(700, 400))
bar!(x_l, eff_l, bar_width=0.03, label="Gauss-Legendre", alpha=0.8, color=:coral)
bar!(x_h, eff_h, bar_width=0.03, label="Gauss-Hermite", alpha=0.8, color=:seagreen)# Plot the integrand g(z)*ϕ(z) and overlay the quadrature nodes
zz = range(μ - 4σ, μ + 4σ, length=200)
integrand = g(zz) .* pdf.(Normal(μ, σ), zz)
p = plot(zz, integrand, lw=2.5, label=L"g(z)\phi(z)", color=:black,
xlabel=L"z = \ln\theta", ylabel="Integrand value",
title="Quadrature nodes overlaid on integrand (n=5)",
size=(700, 400), legend=:topright)
# Trapezoid nodes
scatter!(x_t, g(x_t) .* pdf.(Normal(μ, σ), x_t), ms=8, marker=:circle,
color=:steelblue, label="Trapezoid nodes")
# Legendre nodes
scatter!(x_l, g(x_l) .* pdf.(Normal(μ, σ), x_l), ms=8, marker=:square,
color=:coral, label="Legendre nodes")
# Hermite nodes (plotted at g(z)*ϕ(z) for visual comparison)
scatter!(x_h, g(x_h) .* pdf.(Normal(μ, σ), x_h), ms=8, marker=:diamond,
color=:seagreen, label="Hermite nodes"). . .
Hermite concentrates nodes near the center of the distribution where the integrand has most mass. Trapezoid wastes nodes in the tails where \phi(z) \approx 0
Remember: Gauss-Hermite discretizes the continuous normal into mass points
# Compare the continuous normal pdf with the Hermite discrete approximation for n = 3, 5, 9
ns_plot = [5, 9, 21]
p_disc = plot(layout=(1,3), size=(800, 400),
plot_title="Gauss-Hermite discretization of N($(μ), $(σ))")
for (idx, n) in enumerate(ns_plot)
xn, wn = qnwnorm(n, μ, σ^2)
plot!(p_disc[idx], zz, pdf.(Normal(μ, σ), zz), lw=2, label="Normal pdf", color=:black)
bar!(p_disc[idx], xn, wn ./ 0.05, bar_width=0.05, alpha=0.6,
label="Hermite (n=$n)", color=:seagreen,
xlabel=L"z", title="n = $n", legend=:topleft)
# Scale: divide weights by bar width to get comparable density height
end
p_disc. . .
As n grows, the discrete distribution converges to the continuous one. This is the moment-matching at work
| Method | Best for | Watch out for |
|---|---|---|
| Trapezoid | Bounded domains, non-smooth functions | Slow convergence on smooth problems |
| Gauss-Legendre | General-purpose bounded integrals | Must truncate for unbounded domains |
| Gauss-Hermite | Expectations over normal distributions | Only optimal when density matches |
| Gauss-Laguerre | Integrals with exponential decay | Only optimal for semi-infinite domains |
Others (qnwdist) |
Specialized distributions (Poisson, Chi-Square, Gamma, etc) | Must know the distribution in advance |
qnwdist)Check out the documentation↩︎