Importance sampling method
The GHK algorithm (Geweke, Hajivassiliou and Keane) is an importance sampling method for simulating choice probabilities in the multivariate probit model. These simulated probabilities can be used to recover parameter estimates from the maximized likelihood equation using any one of the usual well known maximization methods (Newton's method, BFGS, etc.). Train has well documented steps for implementing this algorithm for a multinomial probit model. What follows here will apply to the binary multivariate probit model.
Consider the case where one is attempting to evaluate the choice probability of formula_1 where formula_2 and where we can take formula_3 as choices and formula_4 as individuals or observations, formula_5 is the mean and formula_6 is the covariance matrix of the model. The probability of observing choice formula_7 is 
 formula_8
Where formula_9 and,
 formula_10
Unless formula_11 is small (less than or equal to 2) there is no closed form solution for the integrals defined above (some work has been done with formula_12). The alternative to evaluating these integrals closed form or by quadrature methods is to use simulation. GHK is a simulation method to simulate the probability above using importance sampling methods.
Evaluating formula_13 is simplified by recognizing that the latent data model formula_14 can be rewritten using a Cholesky factorization, formula_15. This gives formula_16 where the formula_17 terms are distributed formula_18.
Using this factorization and the fact that the formula_17 are distributed independently one can simulate draws from a truncated multivariate normal distribution using draws from a univariate random normal.
For example, if the region of truncation formula_20 has lower and upper limits equal to formula_21 (including a,b = formula_22) then the task becomes 
 formula_23 
Note: formula_16, substituting:
 formula_25
Rearranging above, 
 formula_26
Now all one needs to do is iteratively draw from the truncated univariate normal distribution with the given bounds above. This can be done by the inverse CDF method and noting the truncated normal distribution is given by, 
formula_27
Where formula_28 will be a number between 0 and 1 because the above is a CDF. This suggests to generate random draws from the truncated distribution one has to solve for formula_29 giving, 
 formula_30
where formula_31 and formula_32 and formula_33 is the standard normal CDF. With such draws one can reconstruct the formula_34 by its simplified equation using the Cholesky factorization. These draws will be conditional on the draws coming before and using properties of normals the product of the conditional PDFs will be the joint distribution of the formula_34, 
 formula_36
Where formula_37 is the multivariate normal distribution.
Because formula_38 conditional on formula_39 is restricted to the set formula_40 by the setup using the Cholesky factorization then we know that formula_37 is a truncated multivariate normal. The distribution function of a truncated normal is, 
 formula_42
Therefore, formula_38 has distribution, 
formula_44
where formula_45 is the standard normal pdf for choice formula_3.
Because formula_47 the above standardization makes each term mean 0 variance 1.
Let the denominator formula_48 and the numerator formula_49 where formula_50 is the multivariate normal PDF.
Going back to the original goal, to evaluate the 
formula_51
Using importance sampling we can evaluate this integral, 
formula_52
This is well approximated by formula_53.
References.
 ;templatestyles src="Reflist/styles.css" / ;