Monte Carlo Integration for Econometric Estimation

Why Integrate by Simulation?

In econometric models with latent variables, an observation's likelihood can require integrating over an unobserved quantity. The integral may have no convenient closed form, or it may be difficult to evaluate numerically with deterministic methods. Monte Carlo integration replaces the integral with an average of function values at random draws. The approximation is simple to construct, and its variability can be quantified when the relevant variance is finite.

The key is to draw from the distribution that appears in the integral. The average then targets the desired integral directly; it is not an average over arbitrary values of the latent variable. The standard error describes how much the estimate would vary across repeated sets of simulation draws, with the data and model parameters held fixed.

From a Likelihood Integral to a Monte Carlo Average

Let \(y\) be an observed outcome, \(U\) a latent variable, and \(\theta\) a fixed parameter value. Suppose \(f(u\mid\theta)\) is the density of the latent variable and \(p(y\mid u,\theta)\) is the conditional likelihood contribution for the observed outcome given \(U=u\). Integrating out \(U\) gives the marginal likelihood contribution

Draw \(U_1,\ldots,U_R\) independently from \(f(u\mid\theta)\), evaluate the conditional likelihood at each draw, and average:

This is a direct application of the expectation formula. In particular, the draws must follow the specified latent-variable distribution at the parameter value being evaluated. If they instead follow a different distribution, this average generally estimates a different expectation.

To define the expectation and form the average, it is enough that \(h(U,\theta)\) be integrable, meaning \(E[|h(U,\theta)|\mid\theta]<\infty\). With independent draws, the law of large numbers then gives convergence of the average to the expectation as \(R\) increases. A finite second moment, \(E[h(U,\theta)^2\mid\theta]<\infty\), is an additional condition for the variance formula and conventional standard-error calculation below. Finite variance is not required merely to compute the simulation average.

Standard Error from Independent Draws

Write \(H_r=h(U_r,\theta)\). If the draws are independent and identically distributed and \(H_r\) has finite variance \(\sigma_h^2\), the variance of their average follows directly from independence:

The middle equality uses the fact that the covariance between distinct independent draws is zero. Since \(\sigma_h^2\) is usually unknown, estimate it with the sample variance of the simulated function values:

This is an estimated standard error across repeated sets of simulation draws, conditional on the target distribution and fixed model inputs. It is not the standard deviation of the individual function values: that is estimated by \(s_h\). Averaging \(R\) independent values reduces the variance by a factor of \(R\) and the standard error by a factor of \(\sqrt{R}\).

Simulated Probabilities

A probability can be written as the expectation of an indicator. For an event \(A\), define \(H_r=\mathbf{1}\{U_r\in A\}\), which equals 1 if draw \(r\) is in the event and 0 otherwise. The mean of these indicators is the simulated probability estimate. Each indicator has variance \(p(1-p)\), where \(p=P(U\in A)\), so independent draws give a standard error that can be estimated by substituting the estimated probability for \(p\).

Estimating Error by Repeated Simulation

A second way to assess Monte Carlo variability is to repeat the entire simulation calculation independently. For each replication \(m=1,\ldots,M\), use \(R\) fresh independent draws and compute an estimate \(\hat I_R^{(m)}\). The sample standard deviation across the \(M\) estimates estimates the standard deviation of one run with \(R\) draws. This approach is useful when the function values are not saved or when the calculation is treated as a complete simulation routine.

Derivation and Practical Checks

The standard-error formula follows from a short variance calculation. For independent \(H_r\) with common variance \(\sigma_h^2\), scaling a random variable by \(1/R\) scales its variance by \(1/R^2\). Variances of independent terms add, so the sum of \(R\) terms has variance \(R\sigma_h^2\). Combining these facts yields

Replacing the unknown \(\sigma_h\) with the sample standard deviation \(s_h\) gives the estimated standard error \(s_h/\sqrt{R}\). This calculation depends on independent draws and finite variance. If draws are dependent, the covariance terms need not be zero, so this formula does not in general give the correct standard error. If the function has infinite variance, the average can still be computed when the expectation exists, but this conventional variance formula and standard-error estimate are not justified.

For a likelihood integral, check that each simulated latent value comes from the latent distribution specified by the model at the parameter value under evaluation. For a probability, check that the event indicator corresponds exactly to the event being estimated. In either case, report the number of independent draws, the estimate, and the Monte Carlo standard error in the units of the estimated quantity.

Common Mistake

A frequent error is to say that finite variance is required just to calculate a Monte Carlo average. It is not. Given simulated values, their arithmetic mean can be calculated, and an integrable function has a well-defined expectation that the average can estimate under the law of large numbers. Finite variance is needed for the variance formula and the conventional standard-error calculation presented here. A separate mistake is to report \(s_h\) itself as the standard error of the average; the standard error is \(s_h/\sqrt{R}\), not \(s_h\).