An accurate and efficient numerical approximation of the multivariate normal (MVN) distribution function is necessary for obtaining maximum likelihood estimates for models involving the MVN distribution. Numerical integration through simulation (Monte Carlo) or number-theoretic (quasi–Monte Carlo) techniques is one way to accomplish this task. One popular simulation technique is the Geweke–Hajivassiliou–Keane MVN simulator. This paper reviews this technique and introduces a Mata function that implements it. It also computes analytical first-order derivatives of the simulated probability with respect to the variables and the variance–covariance parameters.