Simulation of Multivariate Normal Orthant Probabilities: Theoretical and Computational Results

An extensive literature in econometrics and in numerical analysis has considered the problem of evaluating the multiple integral P(B; µ, Ω) = Integralabn(v - µ, Ω)dv = EV1(VcB), where V is a m-dimensional normal vector with mean µ, covariance matrix Ω, and density n(v - µ, Ω) and 1(VcB) is an indicator for the event B = {V|a < V < b}. A leading case of such an integral is the negative orthant probability, where B = {v|v < 0}. The problem is computationally difficult except in very special cases. The multinomial probit (MNP) model used in econometrics and biometrics has cell probabilities that are negative orthant probabilities, with µ and depending on unknown parameters (and, in general, on covariates). Estimation of this model requires, for each trial parameter vector and each observation in a sample, evaluation of P(µ;B) and of its derivatives with respect to µ and Ω. This paper surveys Monte Carlo techniques that have been developed for approximations of P(µ;Ω) and its linear and logarithmic derivatives that limit computation while possessing properties that facilitate their use in iterative calculations for statistical inference: the Crude Frequency Simulator (CFS), Normal Importance Sampling (NIS), a Kernel-Smoothed Frequency Simulator (KFS), Stern’s Decomposition Simulator (SDS), the Geweke–Hajivassiliou–Keane Simulator (GHK), a Parabolic Cylinder Function Simulator (PCF), Deák’s Chi-squared Simulator (DCS), an Acceptance/Rejection Simulator (ARS), the Gibbs Sampler Simulator (GSS), a Sequentially Unbiased Simulator (SUS), and an Approximately Unbiased Simulator (AUS). We also discuss Gauss and FORTRAN implementations of these algorithms and present our computational experience with them. We find that GHK is overall the most reliable method.