No, in fact what you get from my script is the bottom part of the
Cholesky factor. If you want to have the corresponding values for the
correlation matrix, you should adapt the following hansl snippet:...
Just a follow-up. After correcting to get the correlation matrix instead of the bottom part of the Cholesky factor, gretl's trivariate probit script appears to match Limdep with a LOT of uniform draws (and slightly fewer :) halton draws). However, to replicate Limdep's results (or pretty close, anyway), it took a little longer than I said before, but peanuts compared to Limdep.
Gretl took 91 seconds (and my 8-core CPU was shaking at 96.5% during that time - nice job on efficiency).
You should thank Allin for that, he did a very fine job at parallelising GHK.
Thanks Allin - the speed gain from parallelising is quite nice.
I got an "analytical scores" trivariate probit model to work. However, the results only match Limdep to 3-4 decimal places using Genz's tvnl algorithm (and OPG for the covariance matrix), and that's for a relatively small sample (1200 from the example script). Yes, I know that's not terrible considering the nature of multivariate normals CDFs, but it clearly won't do. I am now starting to understand the reason for the GHK simulator, but...I'm determined to find an even faster (and hopefully as accurate, of course) solution, unless I'm just hitting my head against the wall, like the many who have thought about this before me.
In theory, is there some way to use the tnvl algorithm as a starting point for the trivariate cdf since it's fast (or something like Genz's mvndst.f for higher order multivariate normal CDF), then use the GHK simulator to improve the accuracy of said CDF, which could then be used directly in the analytical scores calculation?