Mata's mvnormal(...) suite of functions

[John Mullahy]

On various occasions I've used many of the mvnormal(...) functions to simulate multivariate normal orthant probabilities and related quantities.

A recent simulation exercise convinced me that if your task can afford the additional computation time required then instead of mvnormal(...) using mvnormalqp(...) with a large number of quadrature points is recommended.

(In a nutshell I was quite confident but not 100% certain about a result involving probability inequalities and decided to do "proof-by-Stata." Most of the time, but not always, mvnormal(...) returned results consistent with my expectations. mvnormalqp(...) with the number of quadrature points set at 5,000 returned the "correct" results in each instance.)

Depending on the nature of your task the additional computation time may or may not be problematic. Here's a simple illustration.

mata

reps=1000

m=4
corrmat=J(m,m,.25)+.75*I(m)
mvec=J(1,m,0)

timer_clear()

timer_on(1)

for (j=1;j<=reps;j++) {
 p=mvnormal(mvec,vech(corrmat)')
}

timer_off(1)

timer_on(2)

for (j=1;j<=reps;j++) {
 p=mvnormalqp(mvec,vech(corrmat)',5000)
}

timer_off(2)

timer()

end

which returned

: timer()

---------------------------------------------------------------------------------------------
timer report
  1.       .099 /        1 =      .099
  2.        6.4 /        1 =     6.404
---------------------------------------------------------------------------------------------

Of course there may be a happy middle ground that delivers sufficient accuracy with fewer than 5,000 quadrature points.

[Fang Wang (StataCorp)]

By default, mvnormal() uses 128 quadrature points, which may not be enough in some cases.

Here I tried your example with p = mvnormal() with default points and p5000 = mvnormalqp() with 5000 points, where the relative difference of the results is about 1e-9. If using p1000 = mvnormalqp() with 1000 points, the relative difference with p5000 would be around 1e-13 and the time used is about one third or one fourth compared to the time of computing of p5000 on my machine. If using p500 = mvnormalqp() with 500 points, the relative difference with p5000 would be around 1e-11 and the time used is slightly more than the time used to compute p. So it would depend on how accurate you need in order to balance between accuracy and time.

In addition, the running time and accuracy also depend on the dimension of the problem and upper-limits and correlation matrix, so there may not be a uniformly good number of quadrature points and you may want to try different numbers and see how it works. Let us know if there's any more questions.

[John Mullahy]

Thanks very much for this helpful information, Fang.