Simulation-Based Estimation

Definition

Simulation-based estimation covers methods that replace analytically intractable moment conditions or likelihood contributions with unbiased simulation approximations. The generic problem is evaluating $Eh(U)=\int h(u)f(u)\,du$ where the integral is high-dimensional — e.g., $(J-1)$ dimensions in multinomial probit or $T$ dimensions in dynamic programming with horizon $T$ . Given $R$ draws $u^r\sim f$ , the sample mean $(1/R)\sum_r h(u^r)$ is unbiased with variance $=\mathrm{Var}[h(U)]/R\to 0$ as $R\to\infty$ (Stern 1997).

Key Ideas

Four estimation families:

Method of Simulated Moments (MSM) — replace true moment $h(y,X|\theta)$ with simulated counterpart $\hat{h}$ ; find $\theta$ solving $0=Q'\hat{h}$ ; consistent at any fixed $R$ as $N\to\infty$ (LLN washes out simulation error); asymptotic covariance inflated by $[1+1/R]$ ; GLS-optimal instruments $Q=D^{-1}(y_i)\partial h/\partial\theta$ make MSM asymptotically equivalent to ML.
Maximum Simulated Likelihood (MSL) — simulate $\hat{L}_i=(1/R)\sum_r \hat{L}_i^r$ and maximize $\sum_i\log\hat{L}_i$ ; technically requires $R\to\infty$ as $N\to\infty$ for consistency (Jensen inequality bias: $E\log\hat{P}_i\neq\log EP_i=\log P_i$ ); Börsch-Supan–Hajivassiliou (1993) show bias negligible in practice with smooth simulators; more efficient than MSM.
Method of Simulated Scores (MSS) — simulate the score statistic $\partial\log L_i/\partial\theta$ directly; achieves efficiency of ML plus consistency of MSM; theoretically promising but the hardest to implement; least developed as of 1997.
MCMC — Gibbs sampling, Metropolis-Hastings, and related algorithms form a fourth estimation family; characterized by Stern (1997) as "very expensive relative to MSM and MSL," a 1997 assessment that was overtaken by KSC (1998) mixture approximations, FFBS smoothers, and GPU computing.

Simulator types (ranked by quality):

Frequency (direct) simulator: $\hat{P}_i=(1/R)\sum_r \mathbf{1}[u^r<X_i\beta]$ ; unbiased, $\mathrm{Var}=P_i(1-P_i)/R$ ; but a step function in $\theta$ — discontinuous and non-differentiable, unsuitable for derivative-based optimizers; the naive baseline.
GHK simulator (Geweke 1991; Hajivassiliou 1990; Keane 1994): for MNP, alternates between computing $\hat{p}_k=\Phi[(\text{upper truncation})/\sigma^2]$ analytically and simulating $u_{jk}^*$ from truncated normals; products $\prod_k \hat{p}_k$ accumulate the simulator; strictly bounded $(0,1)$ ; continuous and differentiable in $\theta$ ; importance-sampling interpretation (oversamples joint density correctly via conditioning); best Monte Carlo performance in all known benchmarks.
Decomposition simulator (Stern 1992): decompose $u^*=Z_1+Z_2$ with $Z_1\sim N[0,\lambda]$ , $Z_2\sim N[0,\Omega_j^*-\lambda]$ independent; then $P_j=\int\prod_k\Phi\!\left(\frac{V_{jk}^*-z_{2k}}{\lambda_k}\right)g(z_2)\,dz_2$ simulatable as $(1/R)\sum_r\prod_k\Phi\!\left(\frac{V_{jk}^*-z_{2k}^r}{\lambda_k}\right)$ ; shares GHK's smoothness and boundedness.
Importance sampling: rewrite $Eh(U)=\int [h(u)f(u)/g(u)]\,g(u)\,du$ and draw from well-chosen $g$ ; requirements: same support as $f$ , easy to draw from $g$ , $h(u)f(u)/g(u)$ bounded and smooth; optimal $g\propto |h(u)f(u)|$ achieves zero variance.
Antithetic acceleration: average $h(F^{-1}(z^r))$ and $h(F^{-1}(1-z^r))$ ; for symmetric $F$ and monotone $h$ , covariance term is negative so variance is reduced; reduces simulation loss in $\mathrm{Var}(\hat{\theta})$ to order $N^{-1}$ , making standard-error adjustment unnecessary.

How It Works

MSL roadmap for MNP (Stern 1997, §2.2):

Choose identifiable parameterization of $\Omega$ and initial $\theta_0$ .
Fix $2NJR$ random numbers in separate instrument/estimation files — never change during optimization.
Simulate instruments $Q_{ij}=(1/R)\sum_r(\tilde{P}_{ij}^r)^{-1}(\partial\tilde{P}_{ij}^r/\partial\theta)$ using GHK.
Given instruments and estimation file, find $\theta$ maximizing $\sum_i\sum_j y_{ij}\log\hat{P}_{ij}$ ; use derivative-based optimizer (GHK is smooth).
Update instruments using consistent $\hat{\theta}$ ; reoptimize once for (approximately) efficient estimates.

Why It Matters

Before simulation, empirical economists were restricted to models with analytically tractable integrals — iid Extreme Value errors (logit), or at most Gauss-Hermite with scalar heterogeneity. Simulation unlocks:

MNP with free correlation structure $\Omega$ (no IIA; AR(1) + random effects in Börsch-Supan et al. 1992 improve pseudo- $R^2$ from 40% to 60%).
Dynamic programming with T-dimensional planning horizon (Pakes 1986 patent renewal: $T$ -dim integral for reservation return sequence).
Market entry games with unobserved heterogeneity (Berry 1992: $EN_i$ easy to simulate even when $N_i$ likelihood is not).
Demand systems with interacting consumer/brand unobservables (BLP 1995: $\xi_j$ treated as fixed effect to avoid endogeneity).

Open Questions