GHK Simulator

Definition

The GHK simulator — named for Geweke, Hajivassiliou, and Keane — is a smooth, unbiased Monte Carlo estimator of the probability that a multivariate-normal vector falls in a rectangle, i.e. of orthant/rectangle probabilities $\Pr(a < Z < b)$ with $Z \sim \mathcal{N}(0,\Sigma)$ . It is the workhorse for simulating the choice probabilities and likelihood contributions of high-dimensional multinomial probit and other correlated limited-dependent-variable models, where the required Gaussian integral has no closed form once the dimension exceeds three or four.

Key Ideas

Recursive triangularization. Factor $\Sigma = LL'$ by Cholesky and change variables so the correlated rectangle probability becomes a product of one-dimensional conditional probabilities: $\Pr(a<Z<b) = \prod_{j} \Pr(\ell_{aj} < \eta_j < \ell_{bj} \mid \eta_1,\ldots,\eta_{j-1})$ , each factor a univariate standard-normal interval whose bounds depend recursively on the previously drawn components.
Sequential truncated-normal draws. The simulator walks the recursion, at step $j$ drawing $\eta_j$ from the standard normal truncated to the current interval and accumulating the analytic univariate probability $\Phi(\cdot)-\Phi(\cdot)$ of that interval. The product of these accumulated factors is one draw of the probability; averaging over $R$ such draws gives the estimate.
Smooth and bounded. Because every factor is an analytic normal CDF difference rather than a 0/1 indicator, the estimator is continuous and differentiable in the model parameters $\theta$ and lies strictly in $(0,1)$ — unlike the crude frequency (acceptance–rejection) simulator, which is a step function and can return exactly zero.
Importance-sampling interpretation. The sequential truncation is exactly importance sampling from the conditional distribution, which is why GHK has the lowest variance among the standard simulators in essentially every benchmark.
One of three MSS simulators. In the method of simulated scores taxonomy of Hajivassiliou-McFadden (1998), GHK is the "recursive-triangularization" device, sitting between the acceptance–rejection simulator and the Gibbs-resampling (data augmentation) simulator.

How It Works

Order the alternatives and Cholesky-factor the covariance $\Sigma=LL'$ (the ordering affects finite-sample variance but not consistency).
For $j=1,\ldots,m$ : compute the truncation bounds for $\eta_j$ implied by the rectangle and the already-drawn $\eta_{1:j-1}$ ; record the univariate probability mass $q_j$ of that interval; draw $\eta_j$ from the truncated standard normal.
Multiply the $q_j$ to obtain one simulated probability $\hat p^{(r)}=\prod_j q_j$ .
Average $\hat p = R^{-1}\sum_r \hat p^{(r)}$ ; plug into a simulated maximum likelihood (MSL), method of simulated moments, or Bayesian likelihood evaluation. Consistency of the estimator requires the number of draws $R$ to grow with the sample (faster than $\sqrt N$ for the score-based version).

Why It Matters

Makes correlated discrete choice estimable. Full-covariance multinomial probit — infeasible by quadrature beyond a handful of alternatives — becomes routine once choice probabilities are evaluated by GHK, the standard engine in both classical MSL and Bayesian MNP.
Threads through Bayesian hierarchical models. GHK evaluates the integrated likelihood $f(y_i\mid\beta,D)$ in longitudinal probit models (e.g. the Chib-Carlin blocked sampler), dramatically cutting autocorrelation by marginalising random effects.
A reusable primitive. The same smooth simulator underlies mixed-logit / random-coefficient discrete choice, multivariate probit panels, and simulated Bayes factors — anywhere a Gaussian rectangle probability must be differentiated in $\theta$ .

Minimax exponential tilting (Botev 2017)

Botev (2017) is the modern successor for the same problem — the probability that a multivariate normal falls in a rectangle, and iid sampling from that truncated distribution. It keeps Genz's separation-of-variables (SOV) importance-sampling proposal but adds an exponential tilt whose parameters solve a minimax (saddle-point) problem: choose the tilt that minimizes the worst-case variance of the estimator. The optimization is efficient because it is log-concave, and the tilted proposal has high acceptance in over 100 dimensions (about the cost of one Gibbs cycle) and — crucially — accurate relative error even for rare-event (tiny tail) probabilities, exactly where plain GHK/SOV degrade. It also delivers exact iid draws from the truncated normal, not just a probability estimate.

Open Questions

Ordering and dimension. Simulation variance depends on the Cholesky ordering of alternatives and grows with dimension; adaptive or antithetic orderings help but there is no universally optimal rule. Minimax exponential tilting (Botev 2017) largely tames the dimension/rare-event problem by optimizing the tilt against worst-case variance.
Bias in the log. GHK is unbiased for the probability but the log-likelihood it feeds is a nonlinear transform, so MSL carries a simulation bias of order $1/R$ that must be controlled by scaling $R$ with $N$ .
Competition from quasi-Monte Carlo. Sparse-grid and randomised QMC integration can beat GHK in moderate dimensions; when each is preferable in practice is problem-specific.