Sequential Conditioning and GHK Simulation

Why Sequential Conditioning Helps

Some econometric likelihoods contain probabilities that several latent variables jointly satisfy inequalities. When the latent vector is multivariate normal, the target is a normal probability over a rectangle or another region described by sequential bounds. Direct integration can become difficult as the dimension grows. The Geweke-Hajivassiliou-Keane (GHK) simulator addresses this by conditioning on coordinates in sequence, drawing each coordinate from its appropriate truncated normal distribution, and attaching a probability weight to the resulting path.

The central idea is to replace a rare joint event with draws that are constructed to satisfy its restrictions. The simulator does not simply count how often unrestricted normal draws land inside the region. Instead, it records the conditional probability of satisfying each restriction along a simulated path. Averaging the resulting path weights estimates the desired joint probability.

Definition: Let \(Y\sim N(\mu,\Sigma)\), where \(Y\) is \(d\)-dimensional and \(\Sigma\) is positive definite. For lower and upper bounds \(l\) and \(u\), with \(l_i\leq u_i\), the multivariate normal rectangle probability is \(P(l_i\leq Y_i\leq u_i,\ i=1,\ldots,d)\). The GHK simulator estimates this probability by sequentially generating draws that satisfy the bounds and averaging their conditional-probability weights.

Sequential Conditioning and the GHK Weight

Take a Cholesky factorization \(\Sigma=LL^\mathsf{T}\), where \(L\) is lower triangular with positive diagonal elements. Write \(Y=\mu+LZ\), where the elements of \(Z\) are independent standard normal variables. When coordinates \(1,\ldots,i-1\) have already been drawn, the restriction on coordinate \(i\) becomes a conditional interval for \(Z_i\):

a_i=\frac{l_i-\mu_i-\sum_{j<i}L_{ij}z_j}{L_{ii}}, \qquad b_i=\frac{u_i-\mu_i-\sum_{j<i}L_{ij}z_j}{L_{ii}}.

Because \(Z_i\) is standard normal and independent of the earlier \(Z_j\)'s, the conditional probability of satisfying this interval is \(q_i=\Phi(b_i)-\Phi(a_i)\), where \(\Phi\) is the standard normal cumulative distribution function. Bounds at infinity are allowed, with \(\Phi(-\infty)=0\) and \(\Phi(\infty)=1\). For a uniform draw \(U_i\sim\text{Uniform}(0,1)\), generate

z_i=\Phi^{-1}\!\left(\Phi(a_i)+U_iq_i\right).

This inverse-CDF step draws from a standard normal distribution restricted to \([a_i,b_i]\). The path weight is the product of the conditional probabilities along that path:

W=\prod_{i=1}^{d}q_i.

The same sequence of steps is applied for each simulation replication. The simulated probability is the average of the resulting weights, not the average of indicators for whether unrestricted normal draws satisfy all bounds.

Why Averaging Path Weights Works

For any fixed values of the preceding coordinates, the truncated-normal draw at step \(i\) has density \(\phi(z_i)/q_i\) on \([a_i,b_i]\), where \(\phi\) is the standard normal density. Multiplying that density by the step weight \(q_i\) recovers the unrestricted normal density \(\phi(z_i)\) over the allowed interval. Applying this argument at each step shows that the weight corrects for drawing from the restricted distributions.

Factor the covariance matrix.
Write \(\Sigma=LL^\mathsf{T}\) and represent the normal vector as \(Y=\mu+LZ\).

Find each conditional interval.
Given earlier draws, calculate \(a_i\) and \(b_i\), then the conditional probability \(q_i\).

Draw within the interval and weight the path.
Use the inverse CDF to generate \(z_i\), and multiply the \(q_i\)'s to obtain \(W\).

Average independent path weights.
For \(M\) draws, estimate the probability by \(\hat P_{\mathrm{GHK}}=M^{-1}\sum_{m=1}^{M}W_m\).

More formally, the joint density factors into successive conditional densities. Integrating over the rectangle can therefore be written as an expectation over the sequential truncated draws: the integrand is the product \(W\). If all uniforms across coordinates and replications are mutually independent and uniform, then each \(W_m\) has expectation equal to the target rectangle probability. Thus the ordinary GHK estimate is unbiased for that probability, apart from numerical error in evaluating the normal CDF and its inverse.

P(l_i\leq Y_i\leq u_i,\ i=1,\ldots,d)=E(W), \qquad \hat P_{\mathrm{GHK}}=\frac{1}{M}\sum_{m=1}^{M}W_m.

Worked Examples

Worked Example: An Independent Rectangle Probability

Suppose \(Y_1\) and \(Y_2\) are independent standard normal variables, and the event is \(0\leq Y_1\leq1\) and \(-1\leq Y_2\leq1\). This is a useful check because the Cholesky factor is the identity matrix, so neither conditional bound depends on an earlier draw. The probability is the product of the two one-dimensional interval probabilities.

For the first coordinate, \(q_1=\Phi(1)-\Phi(0)=0.841345-0.500000=0.341345\). For the second, \(q_2=\Phi(1)-\Phi(-1)=0.841345-0.158655=0.682690\). Consequently every GHK path has the same weight:

$$ W=q_1q_2=(0.341345)(0.682690)=0.2330. $$

The exact rectangle probability is also \(0.2330\), to four significant figures. Here simulation introduces no variation in the weights. That is a special feature of independence and fixed marginal bounds, not a general property of GHK.

Worked Example: A Correlated Bivariate Path

Let \(Y\sim N(0,\Sigma)\), with unit marginal variances and correlation \(0.5\). Consider the event \(Y_1>0\) and \(0<Y_2<1\). A Cholesky factor is \(L=\begin{pmatrix}1&0\\0.5&0.866025\end{pmatrix}\), so \(Y_1=Z_1\) and \(Y_2=0.5Z_1+0.866025Z_2\). Use the hypothetical uniform draws \(U_1=0.4\) and \(U_2=0.6\).

At the first step the allowed interval is \(Z_1>0\), so \(q_1=0.5\). The inverse-CDF draw is \(z_1=\Phi^{-1}(0.5+0.4(0.5))=\Phi^{-1}(0.7000)\approx0.5244\).

Given this draw, the restriction on the second coordinate is \(0<0.5(0.5244)+0.866025Z_2<1\). Therefore \(a_2=-0.5(0.5244)/0.866025\approx-0.3028\) and \(b_2=(1-0.5(0.5244))/0.866025\approx0.85194\). Using these unrounded bounds gives \(q_2=\Phi(b_2)-\Phi(a_2)\approx0.42184\). The second draw is

z_2=\Phi^{-1}\!\left(\Phi(a_2)+0.6q_2\right) \approx\Phi^{-1}(0.38103+0.6(0.42184)) =\Phi^{-1}(0.63413)\approx0.3428.

The corresponding value is \(Y_2\approx0.5(0.5244)+0.866025(0.3429)=0.5592\), which lies between zero and one. The path weight is \(W=(0.5000)(0.42184)=0.2109\). This one path contributes a weight to the estimate; it is not itself an estimate of the full bivariate probability.

Worked Example: A Three-Dimensional Path and Its Weight

Continue with a zero-mean normal vector whose covariance has unit variances and pairwise correlations \(0.5\). Its Cholesky factor is \(L=\begin{pmatrix}1&0&0\\0.5&0.866025&0\\0.5&0.288675&0.816497\end{pmatrix}\). Consider the event \(Y_1>0,\ Y_2>0,\ Y_3>0.0507\), and use \(U_1=0.4,\ U_2=0.409,\ U_3=0.779\). These are hypothetical simulation draws.

As before, \(q_1=0.5\) and \(z_1\approx0.5244\). The second-coordinate lower bound is \(a_2=-0.5z_1/0.866025\approx-0.3028\), with no finite upper bound. Thus \(q_2=1-\Phi(a_2)\approx0.6190\). The inverse-CDF draw is approximately \(z_2=0.3430\), since \(\Phi(a_2)+0.409q_2\approx0.6342\). The value \(Y_2=0.5z_1+0.866025z_2\approx0.5592\) satisfies its positive restriction.

The conditional mean contribution to \(Y_3\) from the first two draws is \(0.5z_1+0.288675z_2\approx0.3612\). Hence its standardized lower bound is \(a_3=(0.0507-0.3612)/0.816497\approx-0.3803\), and \(q_3=1-\Phi(a_3)\approx0.6481\). Applying the inverse CDF with \(U_3=0.779\) gives \(z_3\approx1.067\). The resulting component is \(Y_3\approx0.3612+0.816497(1.067)=1.232\), which exceeds \(0.0507\). The full path weight is

$$ W=(0.5000)(0.6190)(0.6482)=0.2006. $$

The conditional draws satisfy all three constraints, while the weight records the probability mass associated with those restrictions. In particular, a draw's satisfying the event does not make its contribution equal to one.

Assessing Numerical Accuracy

The variation across path weights determines the Monte Carlo error of an ordinary GHK estimate. Given \(M\) independent paths, let \(s_W\) be their sample standard deviation. The estimated standard error is \(s_W/\sqrt{M}\). This measures simulation variability for the chosen covariance factorization, ordering, and draw design; it is distinct from sampling uncertainty in the econometric data.

Worked Example: Estimating Simulation Error

Suppose a hypothetical run with \(M=20{,}000\) independent GHK paths gives average weight \(\hat P_{\mathrm{GHK}}=0.2268\) and sample standard deviation \(s_W=0.0990\). Then

\widehat{SE}(\hat P_{\mathrm{GHK}}) =\frac{s_W}{\sqrt{M}} =\frac{0.0990}{\sqrt{20{,}000}} =\frac{0.0990}{141.421} \approx0.000700.

A rough simulation-error range using 1.96 estimated standard errors is \(0.2268\pm1.96(0.000700)\), or approximately \([0.2254,0.2282]\). This is a diagnostic range for Monte Carlo error under the stated simulation, not a confidence interval for an unknown econometric parameter. Repeated independent runs can help reveal instability that one estimated standard error may not make apparent.

The coordinate ordering can affect precision. The exact probability does not depend on the order in which the coordinates are processed, but the conditional bounds and therefore the variation in path weights can. An ordering that avoids very small conditional probabilities may yield more stable weights. In practice, one can compare the estimated standard errors from plausible orderings or compare estimates across independent runs. For very small tail probabilities, direct subtraction of two nearly equal CDF values can lose numerical precision; stable tail-probability calculations are preferable in that situation.

Common Mistakes

Counting only whether a path satisfies the event. The truncated draws are constructed to satisfy the sequential constraints. The probability estimate is the average of the path weights \(W\), not an average of all-one indicators.
Omitting a conditional probability from the weight. Every step contributes its own \(q_i\). The path weight is their product, including the first-coordinate probability.
Using bounds that ignore earlier draws. For \(i>1\), the bounds for \(Z_i\) depend on the previously generated \(z_j\)'s through the Cholesky factor. Reusing unconditional bounds changes the target calculation.
Assuming a different ordering changes the probability. It can change simulation efficiency and the distribution of weights, but a correctly implemented GHK estimator targets the same probability.
Interpreting one path weight as a precise probability estimate. A path is one random contribution. Precision is assessed from the variability of many weights, not from whether an individual weight appears plausible.

Key takeaway: GHK factors a multivariate normal rectangle probability into sequential conditional probabilities. It draws each normal coordinate from its current truncated interval, multiplies the conditional probabilities into a path weight, and averages those weights. Correct conditional bounds and explicit checks of simulation error are essential.

Check Your Understanding

Use the definitions and examples above to answer these questions.

Why does the GHK draw for coordinate \(i\) depend on the previously generated coordinates?
For independent standard normal coordinates with fixed rectangle bounds, explain why each path can have the same weight.
In a two-coordinate calculation, what does the product \(q_1q_2\) represent for a single simulated path?
If the sample standard deviation of GHK path weights stays fixed while the number of independent draws is multiplied by four, how does the estimated standard error change?
Can changing the coordinate ordering change the target probability? What aspect of the simulation can it change?