Consistency and Asymptotic Normality of Simulated Maximum Likelihood

Why Simulation Changes Likelihood Asymptotics

Maximum likelihood is often defined using observation-specific probabilities or densities that require an integral over unobserved variables. Simulated maximum likelihood (SML) replaces those quantities with simulated approximations and maximizes the resulting criterion. This replacement has two distinct consequences: a finite number of draws can change the criterion’s population maximizer, and simulation error can contribute to the estimator’s sampling variation. The consequences depend on how the simulation draws are generated and how their number changes with the sample size.

We study a setting with independent observations and simulation draws generated independently across observations. This makes the simulation contribution explicit; designs that reuse common draws across observations can have different dependence and asymptotic variance. The focus here is on conditions under which the simulated criterion converges to the genuine likelihood criterion and on the first-order distribution when the draw count increases.

Definition: Let $Z_i$, $i=1,\ldots,n$, be independent observations with likelihood contribution $p_i(\theta)$, where $\theta\in\Theta\subset\mathbb{R}^k$. Suppose $p_i(\theta)=E[m_i(\theta,U_i)\mid Z_i]$, where $m_i$ is a simulated contribution and the expectation is over the simulator. With $R_n$ independent draws $U_{ir}$ for observation $i$, define $$ \hat p_{i,R_n}(\theta)=\frac{1}{R_n}\sum_{r=1}^{R_n}m_i(\theta,U_{ir}),\qquad Q_{n,R_n}(\theta)=\frac{1}{n}\sum_{i=1}^n\log\hat p_{i,R_n}(\theta). $$ An SML estimator $\hat\theta_{n,R_n}$ maximizes this simulated criterion, or comes within a stated negligible optimization error of its supremum.

The target criterion is $Q(\theta)=E[\log p_i(\theta)]$, and its unique maximizer is denoted $\theta_0$. In general, $E[\log\hat p_{i,R}(\theta)]$ is not $E[\log p_i(\theta)]$ for finite $R$: the logarithm is nonlinear. Consequently, consistency for $\theta_0$ is not automatic just because each simulated probability is an unbiased estimate of $p_i(\theta)$.

Consistency: Identification and Uniform Approximation

A standard consistency argument requires a compact parameter space, a uniquely identified population maximum, and uniform convergence of the simulated criterion to $Q$. In particular, for every fixed $\varepsilon>0$, identification requires a positive criterion gap:

Q(\theta_0)-\sup_{\theta:\,\|\theta-\theta_0\|\geq\varepsilon}Q(\theta)>0.

Uniform convergence means that $\sup_{\theta\in\Theta}|Q_{n,R_n}(\theta)-Q(\theta)|$ tends to zero in probability. It can be obtained under suitable moment, smoothness, positivity, and uniform law of large numbers conditions when $R_n$ grows, so that simulation approximation error vanishes uniformly. The precise conditions depend on the simulator and the model. If the draw count is fixed, the criterion typically converges instead to a finite-draw population criterion $Q_R$; its maximizer need not be $\theta_0$.

Theorem: Suppose $\Theta$ is compact, $Q$ is continuous with unique maximizer $\theta_0$, and $\sup_{\theta\in\Theta}|Q_{n,R_n}(\theta)-Q(\theta)|\xrightarrow{p}0$. If the estimator’s criterion value is within $o_p(1)$ of the supremum of $Q_{n,R_n}$, then $\hat\theta_{n,R_n}\xrightarrow{p}\theta_0$.

The key inequality behind the result is the identification gap. Outside any fixed $\varepsilon$-neighbourhood of $\theta_0$, the true criterion is strictly worse by a positive amount. Once uniform approximation error and optimization error are smaller than that gap, no approximate maximizer can remain outside the neighbourhood.

Worked Example: What a Uniform Error Bound Does—and Does Not—Show

For illustration, let $\Theta=[0,1]$ and $Q(\theta)=-(\theta-0.5)^2$, whose unique maximizer is $\theta_0=0.5$. Suppose a simulated criterion satisfies $\sup_{\theta\in\Theta}|Q_{n,R_n}(\theta)-Q(\theta)|\leq0.02$. At $\theta_0$, its value is at least $-0.02$. At any point with $|\theta-\theta_0|>0.2$, the true criterion is less than $-0.04$, so its simulated value is less than $-0.02$. Thus no point strictly farther than $0.2$ can strictly outperform $\theta_0$.

At the boundary $\theta=0.3$, however, $Q(0.3)=-0.04$. The allowed error could give $Q_{n,R_n}(0.3)=-0.02$, tying a simulated value of $-0.02$ at $\theta_0$. Therefore this bound does not guarantee that every exact maximizer lies strictly within distance $0.2$. A strict criterion gap, together with approximation error small enough to preserve it, would rule out such a tie. For consistency, the relevant argument uses a fixed positive gap outside each chosen $\varepsilon$-neighbourhood and lets the uniform error tend to zero.

How Simulation Enters the Score

For the asymptotic distribution, assume $p_i(\theta_0)>0$, that the contributions are sufficiently smooth near $\theta_0$, and that the required derivatives and moments exist. Also assume the observations are independent and identically distributed and that, conditional on the data, simulation draws are independent across both observations and draw indices. Define the exact score $s_i=\nabla_\theta\log p_i(\theta_0)$. Under standard likelihood regularity conditions, $E[s_i]=0$ and its covariance is $\Omega=E[s_is_i']$.

The first-order simulation fluctuation in one observation’s score has a useful form. Write $m_i=m_i(\theta_0,U_{ir})$, $m_{i,\theta}=\nabla_\theta m_i(\theta_0,U_{ir})$, and $p_i=p_i(\theta_0)$, $p_{i,\theta}=\nabla_\theta p_i(\theta_0)$. A first-order expansion of the ratio $\nabla\hat p_{i,R}/\hat p_{i,R}$ gives the draw-level term

\psi_{ir}=\frac{m_{i,\theta}}{p_i} -\frac{p_{i,\theta}}{p_i^2}m_i,\qquad E[\psi_{ir}\mid Z_i]=0.

The conditional mean is zero because $E[m_i\mid Z_i]=p_i$ and $E[m_{i,\theta}\mid Z_i]=p_{i,\theta}$. Averaging these terms over $R_n$ draws produces simulation noise. But the logarithm also creates a smaller, generally nonzero finite-draw bias in the expected simulated score, typically of order $1/R_n$. Let $b$ denote its leading coefficient under an expansion $E[\nabla Q_{n,R_n}(\theta_0)]=b/R_n+o(1/R_n)$. The score expansion, under suitable uniform derivative and moment conditions, is

\nabla Q_{n,R_n}(\theta_0) = \frac{1}{n}\sum_{i=1}^n s_i +\frac{1}{nR_n}\sum_{i=1}^n\sum_{r=1}^{R_n}\psi_{ir} +\frac{b}{R_n} +O_p\!\left(\frac{1}{\sqrt n R_n}\right)+o(R_n^{-1}).

When $R_n\to\infty$ and $\sqrt n/R_n\to0$, the root-$n$-scaled remainder vanishes: $\sqrt n\,O_p(1/(\sqrt n R_n))=O_p(1/R_n)=o_p(1)$, and $\sqrt n\,o(R_n^{-1})=o(\sqrt n/R_n)=o(1)$.

Here the second term has covariance of order $1/(nR_n)$ before multiplication by $\sqrt n$. Thus, with independent simulation draws by observation, its contribution to the covariance of the root-$n$ score is of order $1/R_n$. The finite-draw bias has a different implication: after multiplication by $\sqrt n$, its size is governed by $\sqrt n/R_n$. These two effects must not be conflated.

Asymptotic Normality and the Integrated Hessian

Assume that $Q$ has a negative-definite Hessian at $\theta_0$, written $\nabla^2 Q(\theta_0)=-A$, where $A$ is positive definite. Also assume the simulated Hessian converges uniformly in a neighbourhood of $\theta_0$ to the corresponding population Hessian, and that the estimator is consistent and an exact interior maximizer. Then its score is zero.

For a vector parameter, the appropriate expansion integrates the Hessian along the segment joining the true parameter and the estimator. It does not rely on a scalar mean-value theorem with one common intermediate point. Define

\bar H_{n,R_n}= \int_0^1\nabla^2 Q_{n,R_n}\!\left(\theta_0+t(\hat\theta_{n,R_n}-\theta_0)\right)\,dt.

The fundamental theorem of calculus applied along that segment gives the exact first-order score relation (when the estimator’s score is zero):

0=\nabla Q_{n,R_n}(\theta_0) +\bar H_{n,R_n}(\hat\theta_{n,R_n}-\theta_0).

Consistency and uniform Hessian convergence imply $\bar H_{n,R_n}\xrightarrow{p}-A$. Therefore the estimator’s leading error is $A^{-1}\nabla Q_{n,R_n}(\theta_0)$. This integral expansion is valid for a vector parameter and makes clear which Hessian limit supports the rate and distribution calculations.

Theorem: Under the stated regularity conditions, suppose $R_n\to\infty$ and $\sqrt n/R_n\to0$. Then the leading finite-draw score bias is negligible at the root-$n$ scale, and $$ \sqrt n(\hat\theta_{n,R_n}-\theta_0)\xrightarrow{d} N\!\left(0,A^{-1}\Omega A^{-1}\right). $$ Under independent simulation draws across observations, the root-$n$ simulation-noise covariance is of order $1/R_n$, so it vanishes as $R_n\to\infty$. The displayed limit is the same first-order covariance as exact maximum likelihood.

The draw-count condition in this result controls the typical $1/R_n$ bias caused by taking a logarithm of a simulated likelihood contribution. It is a sufficient rate condition under the stated expansion, not a universal requirement for every SML construction: special structure or a different simulation design can change the bias order or dependence. If $R_n$ grows too slowly for $\sqrt n/R_n$ to vanish, a centered root-$n$ limit at $\theta_0$ is not justified by this theorem.

Worked Example: Computing a First-Order Variance

Suppose a scalar parameter $\theta$ has curvature $A=4$, exact-score variance $\Omega=9$, and draw-level simulation variance $\Psi=16$, where $\Psi=E[\psi_{ir}^2]$. Consider $R=25$ independent draws per observation and focus on the leading covariance of the root-$n$ score. The data component contributes $9$; the simulation component contributes $16/25$, since its root-$n$ variance is $\Psi/R$:

\operatorname{Var}\{\sqrt n\,\nabla Q_{n,R}(\theta_0)\} \approx \Omega+\frac{\Psi}{R} =9+\frac{16}{25}=9.64.

The scalar integrated Hessian tends to $-4$, so the corresponding leading estimator variance is

\operatorname{Var}\{\sqrt n(\hat\theta_{n,R}-\theta_0)\} \approx \frac{9.64}{4^2} =\frac{9.64}{16}=0.6025.

For $n=10{,}000$, this approximation gives estimator standard error $\sqrt{0.6025/10{,}000}\approx0.007762$, in the units of $\theta$. This calculation isolates variance contributions and omits any finite-$R$ bias; a centered approximation also requires that bias to be negligible at the relevant scale.

Worked Example: Checking Whether the Bias Is Negligible

Suppose the leading simulated-score bias is $b/R_n$, with $b=0.3$, and compare two draw-count sequences. If $R_n=n^{0.6}$, then

\sqrt n\,\frac{b}{R_n}=0.3n^{0.5-0.6}=0.3n^{-0.1}\longrightarrow0.

Thus this bias vanishes on the root-$n$ scale, as required by the theorem. In contrast, if $R_n=n^{0.4}$, then $\sqrt n\,b/R_n=0.3n^{0.1}$, which does not tend to zero. The usual centered root-$n$ normal approximation is then not supported by this bias calculation. This does not by itself determine the estimator’s alternative rate or limit; those require analyzing the bias and other terms at their appropriate scales.

Common Pitfalls and Practical Meaning

Assuming unbiased simulated probabilities guarantee an unbiased criterion. In general, $E[\hat p_{i,R}]=p_i$ does not imply $E[\log\hat p_{i,R}]=\log p_i$. The finite-draw criterion can have a different maximizer.
Using pointwise convergence where uniform convergence is needed. Pointwise approximation alone does not rule out maximizers moving into parameter regions where approximation is poor. Consistency relies on uniform control together with an identification gap.
Ignoring ties at a boundary of a stated error bound. A non-strict bound can show that a point does not strictly outperform another while still allowing a tie. A strict gap is needed to exclude the point from the maximizer set.
Treating all simulation designs as independent across observations. Reusing draws can induce dependence in simulation errors across observations. The variance order and limiting distribution must then reflect that design rather than use the independent-draw formula automatically.
Dropping the finite-draw bias from a root-$n$ argument. Under an $O(1/R_n)$ score-bias expansion, it is $\sqrt n/R_n$, not merely $1/R_n$, that determines whether the bias disappears at the root-$n$ scale.

These results explain why SML can be both consistent and asymptotically normal even though its likelihood contributions are simulated: the approximation must become sufficiently accurate, and the objective must remain well behaved in the parameter. They also show why uncertainty calculations should distinguish sampling variation from simulation variation and bias, while respecting how the draws were generated.

Check Your Understanding

Use the definitions and conditions above to answer the following questions.

Why does an unbiased simulated likelihood contribution not generally make the simulated log-likelihood criterion unbiased?
What role does the positive criterion gap outside an $\varepsilon$-neighbourhood play in the consistency argument?
With independent simulation draws by observation, what is the order of the simulation-noise covariance in the root-$n$ score?
Why does an $O(1/R_n)$ score bias require checking $\sqrt n/R_n$ for a centered root-$n$ limit?
Why is an integrated Hessian expansion preferable to asserting a single intermediate parameter value for a vector-valued score?