Trivariate*normal cdf function?

[Matthew Burnell]

I am wondering if anyone knows whether there exists a stata function to calculate the trivariate normal cdf? I have had a search, and come up short, so I suspect not. I know there is the mvnormal command, but that is based on simulation.
The reason I ask is because I was hoping to be able to use it in an expression in a margins command (to calculate joint triple correlated probabilities), the same way one can replicate the default joint predictions and se's following the biprobit command by inserting the binorm function into the margins expression option.
Thank you
Matthew

[John Mullahy]

Matthew: I'm 99% sure that you're correct and that there's not such a function in Stata (though am cc'ing Stephen Jenkins who I suspect would know for sure). So my hunch is that you may be left having to resort to some Mata-based calculations using mvnormal which can (as I suspect you already know) provide the joint probability calculations and derivatives thereof. Obtaining se's thereof is a different matter, I'd suppose.

[Stephen Jenkins]

I agree with John, also reserving at least 1% for the possibility of being wrong! [Re mv normal densities, also note lnmvnormalden(M,V,X) function in Stata 14 and later (and maybe earlier).]

Matthew is referring to mvnormal which is part of the mvtnorm package by Michael Grayling and Adrian Mander (ssc describe mvtnorm). [NB Matthew: please ensure you tell Forum readers about the provenance of user-written commands -- see the Forum FAQ). For an earlier discussion of how to calculate multivariate normal probabilities, with code and examples, see Cappellari/Jenkins in Stata Journal 2006, 6(2) -- free download from SJ website.

[Matthew Burnell]

Thank you John and Stephen for your comments. I shall try and think of a work around. The actual probabilities I can of course calculate with mvnormal. As for standard errors, I wonder how far out the SEs for P(A)&P(B&C) compared to P(A)&P(B)&P(C) would be? If similar, the former I could calculate by inserting the binorm function within one of the arguments of another binorm function in the expression option, and then highlighting they are very "approximate" CIs?!

[John Mullahy]

Matthew: I'm not sure what I'd expect from the comparison you suggest, although it would be interesting to know this. I suppose my instinct at this point would be—if your data structure allowed it—to bootstrap the whole process (probability estimation and margins, combined) and compute SEs from those results.

[Matthew Burnell]

Hi John , yes you're right, some bootstrapping will give me an idea of the discrepancy involved. I can't realistically bootstrap the whole process with my particular model as it takes a very long time to fit (bivariate random effects probit model and the predictions are population averaged i.e integrating out the random effect), though I suppose I can easily fit a univariate probit model with no REs and hence simple predictions too, and then see.

[Tiago Pereira]

Do you intend to calculate the Pr(X < t), where

X = max( |Z1|, |Z2|,|Z1|), Z are standard normals, t is a critical value, and Zi ~ mvt(0, sigma), and sigma is a correlation matrix with diagonal 1?

[Matthew Burnell]

Hi Tiago, no not the cdf for the max of Z1, Z2, Z3, but the cdf of the joint probability P(Z1<t1, Z2<t2, Z3<t3) where the Zi's are MVN with a non-diagonal correlation matrix i.e correlations existing.
However, as I mention at the top, I can do the calculation with the mvnormal command, but I was hoping for a proper stata function (like binorm for bivariate normals) that can be used like any other (simple) function in command such margins, predictnl or nlcom.

[Matthew Burnell]

Hi John just in case you are interested (or indeed anyone else reading the thread), I did some bootstrapping on the simpler bivariate probit model (not univariate as I accidentally wrote above), which still took a while because calculating the prediction for all using mvnormal takes some time. Obviously it is odd I talk about trivariate joint prediction with a bivariate model, but I don't think that matters for the purpose of this experiment: I just need three predictions (even the same one) and then supply a correlation between them - which doesn't have to be real, I can just use the error correlation of the bivariate terms as long as I refer to it using the proper stata estimation name i.e. exp(2*_b[athrho:_cons])-1)/(exp(2*_b[athrho:_cons])+1) ), in order to calculate SEs via margins delta methods (and same for the linear predictor part too).

Anyway the real dataset had 926 samples, and overall probability for all 3 (the same) predictions was 0.774, meaning a joint prediction of 0.464 if they were independent. I instead pretended they were each correlated with each other with rho=0.372, the actual bivariate error correlation:
So for the fully margins-method i.e. P(A,P(B,C)) - that is, the 'binorm within binorm' version, the "joint" probability was 0.569 with a SE of .0218296.
Using a fully-bootstrapped procedure (1000 replications) the proper joint calculation of P(A, B, C) using mvnromal was 0.574 with a SE of .0233272. However, the SE of the margins-method when also bootstrapped was .0233017.
So, in this particular example at least, the actual point estimates are very similar (though I wouldn't suggest replacing one for the other), but the SEs when estimated with the same method (bootstrapping) gives near identical results. I imagine the discrepancy will grow as you start to further agglomerate the approximation e.g. P(P(A,B),P(C,D)) for P(A,B,C,D) etc.