Estimating a Latent-Variable Model by Simulation

Why Simulate a Latent-Variable Likelihood?

A latent variable can represent persistent but unobserved differences across individuals. For example, repeated binary choices by the same person may be correlated because of a stable preference that is not measured. When the latent effect is integrated out, the likelihood contribution can involve an integral without a convenient closed form. Simulated maximum likelihood (SML) approximates that integral numerically and then maximizes the resulting approximate likelihood.

The goal here is to work through that procedure for a binary panel model with a normally distributed random intercept. The example is deliberately compact so the likelihood calculations can be checked by hand; a real application would generally use more individuals, more simulation draws, and careful optimization diagnostics.

Model and Likelihood

Definition: For individuals $i=1,\ldots,n$ observed at periods $t=1,\ldots,T_i$, let $y_{it}\in\{0,1\}$ be the outcome and $x_{it}$ a vector of observed covariates. The random-intercept probit model is $$ P(y_{it}=1\mid x_{it},b_i)=\Phi(x_{it}'\beta+b_i), $$ where $\Phi$ is the standard normal cumulative distribution function, $\beta$ is a coefficient vector, and $b_i\sim N(0,\sigma_b^2)$ is an individual-specific latent effect. Conditional on $b_i$ and the covariates, outcomes are independent across periods. The random effects are independent across individuals and independent of the covariates.

The probit error scale is normalized by the standard normal link. Consequently, $\sigma_b$ describes the standard deviation of the latent individual effect relative to that normalized scale; it is not a freely estimated residual scale. Write the parameter vector as $\theta=(\beta',\sigma_b)'$, with $\sigma_b\geq 0$.

For a fixed value of the latent effect, define the conditional probability of the observed outcome as $$ q_{it}(b_i;\beta)=\Phi(x_{it}'\beta+b_i)^{y_{it}}\left[1-\Phi(x_{it}'\beta+b_i)\right]^{1-y_{it}}. $$ Conditional independence means the probability of the individual’s entire outcome history is the product of these period-specific probabilities. For $\sigma_b>0$, integrating over the random-effect distribution gives the individual likelihood contribution

L_i(\beta,\sigma_b)= \int_{-\infty}^{\infty} \left\{\prod_{t=1}^{T_i}q_{it}(b;\beta)\right\} \frac{1}{\sigma_b}\phi\!\left(\frac{b}{\sigma_b}\right)\,db,

where $\phi$ is the standard normal density. At the boundary $\sigma_b=0$, the random effect is degenerate at zero, not a normal density with a zero denominator. Define the likelihood there directly:

L_i(\beta,0)= \prod_{t=1}^{T_i} \Phi(x_{it}'\beta)^{y_{it}} \left[1-\Phi(x_{it}'\beta)\right]^{1-y_{it}}.

Assuming individuals are independent, the sample likelihood is $\prod_i L_i(\theta)$, and the log likelihood is $\ell(\theta)=\sum_i\log L_i(\theta)$. The dependence among repeated outcomes for the same person is represented by integrating one shared $b_i$ across all periods. Integrating separately in every period would instead remove that shared source of dependence.

Constructing the Simulated Likelihood

Let $z_{ir}$, for $r=1,\ldots,R$, be standard normal draws, so that $b_{ir}=\sigma_b z_{ir}$ has the required random-effect distribution. For $\sigma_b>0$, approximate the integral for individual $i$ by averaging the conditional history probabilities over the draws:

\widehat L_i(\theta)= \frac{1}{R}\sum_{r=1}^{R} \prod_{t=1}^{T_i} q_{it}(\sigma_b z_{ir};\beta), \qquad \widehat\ell_R(\theta)=\sum_{i=1}^{n}\log\widehat L_i(\theta).

The same expression has a well-defined limit at $\sigma_b=0$: every draw then gives $b_{ir}=0$, and $\widehat L_i(\beta,0)=L_i(\beta,0)$. The simulation draws may be random or generated by a deterministic design such as a suitable quasi-Monte Carlo rule. For optimization, keeping the draw design fixed as parameters change usually makes the objective smoother than drawing new numbers at every evaluation. This is the common-random-numbers approach.

The order of operations matters. Each draw first gives the probability of the whole observed history; those probabilities are averaged; only then is the logarithm taken. The average of log probabilities is a different quantity and does not approximate the log of the integral.

Worked Example: One Individual’s Two-Draw Contribution

Consider one hypothetical individual observed twice, with outcomes $y_i=(1,0)$. Suppose the two linear indices excluding the random effect are $-0.2$ and $0.5$, and evaluate the likelihood at $\sigma_b=0.6$. Use two standard normal draws $z_{i1}=-1$ and $z_{i2}=1$, so the corresponding random effects are $-0.6$ and $0.6$.

For the first draw, the two full indices are $-0.8$ and $-0.1$. The probability of observing a success and then a failure is $$ \Phi(-0.8)\{1-\Phi(-0.1)\} =\Phi(-0.8)\Phi(0.1) \approx 0.211855\times0.539828 \approx 0.1144. $$

For the second draw, the full indices are $0.4$ and $1.1$. Its history probability is $$ \Phi(0.4)\{1-\Phi(1.1)\} =\Phi(0.4)\Phi(-1.1) \approx 0.655422\times0.135666 \approx 0.08892. $$

Using the unrounded draw-specific probabilities, the two-draw simulated likelihood contribution is approximately $$ \widehat L_i=\frac{0.114365+0.0889185}{2}\approx 0.1016, \qquad \log\widehat L_i\approx -2.286. $$ This is a likelihood contribution at the specified parameter values and draw design, not an estimate of the parameters by itself.

From Contributions to an Objective

For a sample, compute one simulated contribution per individual and add their log contributions. Do not multiply the individual draws across people or average the likelihood contributions across people: the independent sampling units are the individuals, so their likelihood contributions multiply and their log contributions add. Within each individual, the same latent-effect draw applies to every period.

Worked Example: A Two-Individual Simulated Log Likelihood

Continue with two simulation draws $z=-1,1$ and $\sigma_b=0.6$. Individual 1 has the outcomes and indices from the previous calculation, so $\widehat L_1\approx0.101642$. For a second hypothetical individual, let $y_i=(0,1)$, with indices excluding the random effect $0.3$ and $-0.4$.

With $z=-1$, the full indices are $-0.3$ and $-1.0$. The history probability is $$ \{1-\Phi(-0.3)\}\Phi(-1.0) =\Phi(0.3)\Phi(-1.0) \approx 0.617911\times0.158655 \approx 0.09803. $$ With $z=1$, the indices are $0.9$ and $0.2$, giving $$ \{1-\Phi(0.9)\}\Phi(0.2) =\Phi(-0.9)\Phi(0.2) \approx 0.184060\times0.579260 \approx 0.1066. $$ Thus $\widehat L_2\approx(0.09804+0.1066)/2\approx0.1023$. The two-individual simulated log likelihood is $$ \widehat\ell_R =\log(0.101642)+\log(0.1023) \approx -2.2863-2.2798 \approx -4.566. $$ The average of the two draw-specific log probabilities for individual 1 would not give $\log(0.101642)$; the required operation is the logarithm of the average probability.

Optimization and a Short Justification

SML estimates the parameters by maximizing $\widehat\ell_R(\theta)$ over the admissible parameter space. In practice, an optimizer needs parameter constraints, starting values, and a fixed simulation design during each optimization run. One can optimize $\sigma_b$ subject to $\sigma_b\geq0$, or use a transformation such as $\sigma_b=\exp(\gamma)$ for interior solutions; the latter cannot attain the boundary value zero. Several starting values and checks of the final objective and gradient are useful because a numerical stopping message alone does not establish that a relevant maximum was found.

The integral approximation has a direct justification. For fixed $\theta$, define $g_i(z;\theta)=\prod_t q_{it}(\sigma_b z;\beta)$. If the draws are independent standard normal variables, then

E[g_i(Z;\theta)] =\int g_i(z;\theta)\phi(z)\,dz =L_i(\theta), \qquad E\!\left[\frac{1}{R}\sum_{r=1}^{R}g_i(z_{ir};\theta)\right] =L_i(\theta).

Thus the simulated likelihood contribution is an unbiased Monte Carlo approximation to the likelihood contribution at a fixed parameter. But taking logs changes that property. Since the logarithm is concave, Jensen’s inequality gives $$ E[\log \widehat L_i(\theta)]\leq \log E[\widehat L_i(\theta)] =\log L_i(\theta). $$ Therefore the simulated log likelihood is generally not an unbiased estimate of the exact log likelihood. Increasing the number of draws reduces integration error under standard conditions, but finite-draw optimization can still be affected by simulation error.

A practical numerical detail is to avoid underflow when probabilities are very small. Compute each draw’s log history probability, use a log-sum-exp calculation to form the logarithm of their average, and retain adequate numerical precision. A low draw count may make computation quick but can distort the objective or its derivatives. Draw-count and design choices should be checked for stability, especially by verifying that estimates and reported uncertainty do not change materially when the simulation is made more accurate.

Reporting Estimates and Simulation-Aware Uncertainty

An SML report should identify the model, the simulation method and draw count, the optimization procedure, and the uncertainty calculation. Sampling uncertainty reflects which individuals were observed. Simulation uncertainty reflects approximation of the likelihood by a finite simulation design. These sources are conceptually distinct; reporting only a conventional standard error can understate uncertainty when simulation error is not negligible.

Under regularity conditions and sufficiently accurate simulation, an estimator’s sampling covariance can be estimated from the curvature of the log likelihood or from a score-based sandwich estimator. With finite simulation, simulation error can be assessed by repeating the optimization with independent simulation designs while holding the observed data fixed, or by a suitable score-based variance calculation. If sampling and simulation components are estimated as independent variance components, their variances—not their standard errors—are added. In applications, use a method appropriate to the simulation design and estimation procedure rather than assuming this independence automatically.

The boundary $\sigma_b=0$ also deserves attention: a variance component at zero is on the edge of the parameter space, so ordinary interior-point normal approximations for that parameter need not apply. The example below reports an interior estimate and uses approximate intervals only as a transparent illustration. A real analysis should justify its inferential method, particularly when the estimate is close to zero.

Worked Example: Combining Sampling and Simulation Uncertainty

Suppose an optimizer applied to a hypothetical panel data set reports $\widehat\beta_x=0.700$ and $\widehat\sigma_b=0.550$. A sampling-uncertainty calculation gives standard errors $0.120$ and $0.160$, respectively. Independent simulation-design replications estimate simulation standard deviations of $0.030$ for $\widehat\beta_x$ and $0.040$ for $\widehat\sigma_b$. For illustration, suppose the two components are treated as independent.

The combined standard errors are the square roots of the sums of component variances: $$ SE(\widehat\beta_x)=\sqrt{0.120^2+0.030^2} =\sqrt{0.0153}\approx0.124, $$ $$ SE(\widehat\sigma_b)=\sqrt{0.160^2+0.040^2} =\sqrt{0.0272}\approx0.165. $$ Using a normal approximation for this interior estimate, the illustrative 95% interval for $\beta_x$ is $$ 0.700\pm1.96(0.124)=[0.457,0.943], $$ and for $\sigma_b$ it is $$ 0.550\pm1.96(0.165)=[0.227,0.873]. $$ The first interval concerns the change in the latent probit index associated with a one-unit increase in the covariate, holding other covariates fixed; it is not a percentage-point change in outcome probability. The second interval describes the latent-effect standard deviation on the normalized probit scale. These intervals are approximate, and the stated method requires special care if the estimate is near the zero boundary.

Common Pitfalls and Practical Checks

Averaging log probabilities. This substitutes the mean of draw-specific log probabilities for the log of the simulated integral. Average the conditional history probabilities first, then take their logarithm.
Using separate latent effects by period. This fails to represent a shared individual effect and its induced dependence across repeated outcomes. Use one draw of $b_i$ for the individual’s full history.
Ignoring the boundary. The density formula with $1/\sigma_b$ applies only when $\sigma_b>0$. At zero, use the likelihood with $b_i=0$ directly.
Treating simulation error as sampling error. Repeated simulation designs can reveal integration sensitivity, but they do not replace uncertainty due to sampling individuals. State how each component was estimated.
Interpreting a probit coefficient as a probability change. A coefficient shifts the latent index. The probability change also depends on the covariate values and the probit link.

A credible SML analysis connects the model assumptions to the likelihood, the likelihood to a transparent simulation calculation, and the numerical maximum to a defensible uncertainty statement. Check convergence, the boundary behavior, sensitivity to starting values, and stability under a more accurate simulation design before interpreting estimates.

Check Your Understanding

Answer these questions using the definitions and calculations in this tutorial.

Why does the likelihood for a person with repeated outcomes integrate one random effect across all periods rather than a separate random effect in each period?
At $\sigma_b=0$, what likelihood contribution replaces the normal-density integral, and why?
Explain why the simulated likelihood averages draw-specific history probabilities before taking a logarithm.
If the estimated sampling standard error is $0.10$ and the estimated independent simulation standard deviation is $0.05$, what is the combined standard error under a variance-components calculation?
What does a coefficient in the random-intercept probit model describe, and why is it not directly a percentage-point change in outcome probability?