Numerical integration
University of Illinois Urbana-Champaign
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\)
Quite 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…)
Integrating with 5 nodes:328.125
Integrating with 50 nodes:333.29862557267813
Integrating with 100 nodes:333.3248307995783Decent 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_i\)s and \(x_i\)s 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…)
QuantEcon.jl: Gaussian quadratureIntegrating with 1 node:250.0
Integrating with 2 nodes:333.3333333333334All 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 quadratureDo we know what is the analytic solution to this integral?
QuantEcon.jl0.78674618361748420.78424076661781570.785398164063432Computing 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
Legendre 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)\)
Hermite 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-15Notice 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_discAs \(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)