Integration

Paul Schrimpf

2026-02-25

\[ \def\Er{{\mathrm{E}}} \def\En{{\mathbb{E}_n}} \def\cov{{\mathrm{Cov}}} \def\var{{\mathrm{Var}}} \def\R{{\mathbb{R}}} \def\arg{{\mathrm{arg}}} \newcommand\norm[1]{\left\lVert#1\right\rVert} \def\rank{{\mathrm{rank}}} \newcommand{\inpr}{ \overset{p^*_{\scriptscriptstyle n}}{\longrightarrow}} \def\inprob{{\,{\buildrel p \over \rightarrow}\,}} \def\indist{\,{\buildrel d \over \rightarrow}\,} \DeclareMathOperator*{\plim}{plim} \DeclareMathOperator*{\argmax}{argmax} \DeclareMathOperator*{\argmin}{argmin} \]

Example: Random Coefficients Demand

Example: Random Coefficients Demand

• Berry, Levinsohn, and Pakes (1995)
• Consumers choose product: \[ j = \argmax_{j \in \{0, ..., J\}} x_{jt}' (\bar{\beta} + \Sigma \nu_i) + \xi_{jt} + \epsilon_{ijt} \]
  • \(\nu_i \sim N(0,I_k)\), \(\epsilon_{ijt} \sim\) Type I Extreme Value
  • Unobserved demand shock \(\xi_{jt}\)

Example: Random Coefficients Demand

• Aggregate demand: \[ s_{jt} = \int \frac{e^{x_{jt}'(\bar{\beta} + \Sigma \nu) + \xi_{jt}}} {\sum_{k = 0}^J e^{x_{kt}'(\bar{\beta} + \Sigma \nu) + \xi_{kt}} } dF\nu \]

Example: Random Coefficients Demand

• Instruments \(Z_{jt}\) with \(E[\xi_{jt} Z_{jt}] = 0\)
• \(g(s_{jt},x_{jt}, Z_{jt}, \bar{\beta},\Sigma) = \left(\delta_{jt}(s_{\cdot t}, x_{\cdot t},\beta,\Sigma) - x_{jt}'\bar{\beta}\right) Z_{jt}\)
• where \(\delta_{jt}\) solves \[ s_{jt} = \int \frac{e^{\delta_{jt} + x_{jt}'\Sigma \nu}} {\sum_{k = 0}^J e^{\delta_{kt} + x_{kt}'\Sigma \nu}} dF\nu \]

Code

Shares: version 1

function share(δ, Σ, dFν, x, ∫)
    function shareν(ν)
        s = exp.(δ .+ x * Σ * ν)
        return (s ./ sum(s))
    end
    return (∫(shareν, dFν))
end

share (generic function with 1 method)

• ∫(shareν, dFν) is not a function that exists, will create it next
• clear?
• correct?
• robust?

Shares: version 2

@doc raw"""
    share(δ, Σ, dFν, x, ∫)

Computes shares in random coefficient logit with mean tastes `δ`, observed characteristics `x`, unobserved taste distribution `dFν`, and taste covariances `Σ`.

# Arguments

- `δ` vector of length `J`
- `Σ` `K` by `K` matrix
- `dFν` distribution of length `K` vector
- `x` `J` by `K` array
- `∫` function that takes a function and a distribution and returns the integral of the function with respect to the distribution


# Returns

- vector of length `J` consisting of $s_1$, ..., $s_J$
"""
function share(δ, Σ, dFν, x, ∫)
    J, K = size(x)
    (length(δ) == J) || error("length(δ)=$(length(δ)) != size(x,1)=$J")
    K, K == size(Σ) || error("size(x,1)=$J != size(Σ)=$(size(Σ))")
    function shareν(ν)
        s = δ .+ x * Σ * ν
        s .-= maximum(s)
        s .= exp.(s)
        s ./= sum(s)
        return (s)
    end
    return (∫(shareν, dFν))
end

Main.Notebook.share

• adds documentation
• some error checking on inputs
• protects against overflow in exp
• Are these changes needed & desirable?

Integration: Monte-Carlo

using Distributions, Statistics
∫mc(f, dx; ndraw=100) = mean(f(rand(dx)) for i in 1:ndraw)

∫mc (generic function with 1 method)

module EC

mutable struct EvaluationCounter
    f::Function
    n::Int
end

EvaluationCounter(f) = EvaluationCounter(f, 0)

function (ec::EvaluationCounter)(x)
    ec.n += 1
    return ec.f(x)
end

import Base: show
function show(io::IO, ec::EvaluationCounter)
    print(io, "evaluated $(ec.f) $(ec.n) times")
end

function reset!(ec::EvaluationCounter)
    ec.n = 0
end

end

using Distributions
dx = Normal(0, 1)

f = EC.EvaluationCounter(x -> x^2)
trueint = 1.0

EC.reset!(f)
intmc = ∫mc(f, dx, ndraw=100)
@show intmc, trueint, intmc - trueint
@show f
EC.reset!(f)

(intmc, trueint, intmc - trueint) = (1.0611575919429865, 1.0, 0.06115759194298653)
f = evaluated #6 100 times

0

Integration: Quasi-Monte-Carlo

• Use “low discrepency sequences” to reduce variance of Monte-Carlo integrals
• Error is \(O\left(\frac{\log(n)^d}{n}\right)\) for integrating \(d\) dimensional function with \(n\) draws
• See Owen (2023) for details

using Sobol
import Base.Iterators: take
function ∫s(f, dx::AbstractMvNormal; ndraw=100)
    marginals = [Normal(dx.μ[i], sqrt(dx.Σ[i, i])) for i in 1:length(dx)]
    invcdf(x) = quantile.(marginals, x)
    ss = skip(SobolSeq(length(dx)), ndraw)
    mean(f(invcdf(x)) for x in take(ss, ndraw))
end

function ∫s(f, dx::Normal; ndraw=100)
    invcdf(x) = quantile(dx, x)
    ss = skip(SobolSeq(length(dx)), ndraw)
    mean(f(invcdf(x[1])) for x in take(ss, ndraw))
end

∫s (generic function with 2 methods)

f = x -> x^2
S = 1_000
@show mean(abs(∫s(f, dx, ndraw=100) - trueint) for s in 1:S)
@show mean(abs(∫mc(f, dx, ndraw=100) - trueint) for s in 1:S)

mean((abs(∫s(f, dx, ndraw = 100) - trueint) for s = 1:S)) = 0.02875751885407857
mean((abs(∫mc(f, dx, ndraw = 100) - trueint) for s = 1:S)) = 0.1173042241441584

0.1173042241441584

Integration: Gaussian Quadrature

\[ \int f(x) dF(x) \approx \sum_{i=1}^n w_i f(x_i) \] - \(x_i\) and \(w_i\) chosen to make approximation exact for polynomials of degree \(\leq 2n-1\) - Error is \(O(n^{-r/d} )\) for \(r\) times differentiable \(f: \R^d \to \R\)

Integration: Gaussian Quadrature

using FastGaussQuadrature, LinearAlgebra
import Base.Iterators: product, repeated
function ∫q(f, dx::MvNormal; ndraw=100)
    n = Int(ceil(ndraw^(1 / length(dx))))
    x, w = gausshermite(n)
    L = cholesky(dx.Σ).L
    sum(f(√2 * L * vcat(xs...) + dx.μ) * prod(ws)
        for (xs, ws) ∈ zip(product(repeated(x, length(dx))...),
        product(repeated(w, length(dx))...))
    ) / (π^(length(dx) / 2))
end

∫q (generic function with 1 method)

Integration: Sparse Grid Quadrature

• Naive Gaussian quadrature in \(d\) dimensions takes product of 1-dimensional quadrature rules – curse of dimensionality
• Sparse grids use subset of 1-dimensional quadrature points
• Error \(O(n^{-r}(\log n)^{(d-1)(r+1)})\)

Integration: Sparse Grid Quadrature

using SparseGrids
function ∫sgq(f, dx::MvNormal; order=5)
    X, W = sparsegrid(length(dx), order, gausshermite, sym=true)
    L = cholesky(dx.Σ).L
    d = length(dx)
    rootdetΣ = abs(det(L))
    wnorm = (π)^(d/2) #* rootdetΣ
    sum(f(√2 * L * x + dx.μ) * w for (x, w) ∈ zip(X, W)) / wnorm
end

∫sgq (generic function with 1 method)

Integration: Adaptive Cubature

• Adaptive integration methods compute \[ \int f(x) dF(x) \approx \sum_{i=1}^n w_i f(x_i) = I_n(f) \] for increasing \(n\) until \(|I_n(f) - I_m(f)|\) is small
• Provides explicit error estimate
• Be careful when nested inside other interative methods

Integration: Adaptive Cubature

using HCubature
function ∫cuba(f, dx; rtol=1e-4)
    D = length(dx)
    x(t) = t ./ (1 .- t .^ 2)
    Dx(t) = prod((1 .+ t .^ 2) ./ (1 .- t .^ 2) .^ 2)
    hcubature(t -> f(x(t)) * pdf(dx, x(t)) * Dx(t), -ones(D), ones(D), rtol=rtol)[1]
end

∫cuba (generic function with 1 method)

Choosing an Integration Method

flowchart TD
    Dim[Dimension of Integral, d?]
    Dim -->|d <= 2| Quadrature[Gaussian Quadrature<br>Adaptive Cubature]
    Dim -->|2 < d < 10| Moderate[Moderate Dimension]
    Dim -->|d >= 10| High[High Dimension]

    Moderate --> Smooth{Function Smooth and<br>Well-Behaved?}
    Smooth -->|Yes| Sparse[Sparse Grids / Smolyak]
    Smooth -->|No| QMC1[Quasi-Monte Carlo]

    High --> Dist{Complex Target<br>Distribution?}
    Dist -->|Yes| MCMC[MCMC / Importance Sampling]
    Dist -->|No| QMC2[Quasi-Monte Carlo<br>Monte Carlo]

Testing

using Test
@testset "Integration" begin
    testfunctions = (
        (x->1,"x->1"),
        (x->sum(x),"x->sum(x)"),
        (x->norm(x)^2, "x->norm(x)^2"),
        (x->exp(norm(x))/(1+exp(norm(x))), "x->exp(norm(x))/(1+exp(norm(x)))")
    )

    C = rand(4,4)
    tol = 1e-2
    distributions = (
        MvNormal(zeros(2),I(2)),
        MvNormal(rand(4), C'*C + I(4))
    )
    baselinemethod=(f,dx)->∫mc(f,dx,ndraw=ceil(Int,10/tol^2))
    methods = (((f,dx)->∫s(f,dx,ndraw=10_000), "QMC"),
               ((f,dx)->∫q(f,dx,ndraw=80), "gauss-hermite") ,
               ((f,dx)->∫sgq(f,dx,order=5), "sparse grid"),
               (∫cuba, "cuba") )
    for ∫ ∈ methods
        for (fs,dx) ∈ Iterators.product(testfunctions, distributions)
            f=fs[1]
            target = baselinemethod(f,dx)
            result = ∫[1](f,dx)
            atol = 2*max(abs(target)*tol,tol)
            #@test target ≈ result atol=atol
            if norm(target-result)>atol
                x = rand(dx)
                println("failed on f=$(fs[2]), dx=$(dx), ∫=$(∫[2])")
                println(target," ", result)
            end
        end
    end
end

failed on f=x->exp(norm(x))/(1+exp(norm(x))), dx=DiagNormal(
dim: 2
μ: [0.0, 0.0]
Σ: [1.0 0.0; 0.0 1.0]
)
, ∫=sparse grid
0.7589297236273503 0.6050491899328833
failed on f=x->sum(x), dx=FullNormal(
dim: 4
μ: [0.3253255867557334, 0.22854123078291455, 0.792631687270502, 0.4018903696552496]
Σ: [2.8965965547240584 2.286883364057918 2.1230930571767 1.232741682274224; 2.286883364057918 4.053132264887253 2.7880861184516266 1.6986855638063563; 2.1230930571767 2.7880861184516266 3.5831543559483365 1.5638761846196862; 1.232741682274224 1.6986855638063563 1.5638761846196862 2.036434752417636]
)
, ∫=sparse grid
1.711156567835575 1.7483888744643947
failed on f=x->exp(norm(x))/(1+exp(norm(x))), dx=FullNormal(
dim: 4
μ: [0.3253255867557334, 0.22854123078291455, 0.792631687270502, 0.4018903696552496]
Σ: [2.8965965547240584 2.286883364057918 2.1230930571767 1.232741682274224; 2.286883364057918 4.053132264887253 2.7880861184516266 1.6986855638063563; 2.1230930571767 2.7880861184516266 3.5831543559483365 1.5638761846196862; 1.232741682274224 1.6986855638063563 1.5638761846196862 2.036434752417636]
)
, ∫=sparse grid
0.9262812346551518 0.9477629126956173
Test Summary: | Total  Time
Integration   |     0  9.0s

Test.DefaultTestSet("Integration", Any[], 0, false, false, true, 1.772051036770856e9, 1.772051045811771e9, false, "/home/paul/ECON622/qmd/integration.qmd", Random.Xoshiro(0x12cae4d216e78aea, 0xe01a1dabd6be4d08, 0x5d6dee1cf982bc93, 0xf52cac47c43e46c0, 0x3a779a5cd8641021))

References

Berry, Steven, James Levinsohn, and Ariel Pakes. 1995. “Automobile Prices in Market Equilibrium.” Econometrica 63 (4): pp. 841–890. http://www.jstor.org/stable/2171802.