Hello,

This is a "big" question, and I apologize in advance for the long message. I'm new to this list, so I'm not sure if this is the right place to ask about this.

I'm attempting to extend some of gretl's discrete choice models in my own programming project (c, mingw, win32). One of the models I've been trying to get working is trivariate probit (and hopefully eventually multivariate probit) - I've adapted Alan Genz's matlab script for c (http://www.math.wsu.edu/faculty/genz/software/matlab/tvnl.m) for calculating trivariate normal cdfs (which matches his results), and I've been using the code for gretl's bivariate model as a basis to attempt to extend to three dimensions. In the bivariate case, there's a place in the code (biprobit.c) where the first partials of the bivariate cdf are computed (at least I'm pretty sure that's what d1 and d2 are):

    ca = cosh(bp->arho);
    sa = sinh(bp->arho);

    u_ba = (ca*b - ssa*a);
    u_ab = (ca*a - ssa*b);

    d1 = exp(-0.5*a*a) * normal_cdf(u_ba) / (P * SQRT_2_PI);
    d2 = exp(-0.5*b*b) * normal_cdf(u_ab) / (P * SQRT_2_PI);

After taking care of all the sign switches, I imagine the code for the trivariate case is something along the lines of (but obviously not correct since it's not working):

       ca12 = cosh(tp->arho12);
       sa12 = sinh(tp->arho12);
       ca13 = cosh(tp->arho13);
       sa13 = sinh(tp->arho13);
       ca23 = cosh(tp->arho23);
       sa23 = sinh(tp->arho23);

        u_ba = (ca12*b - ssa12*a);
        u_ab = (ca12*a - ssa12*b);
        u_ca = (ca13*c - ssa13*a);
        u_ac = (ca13*a - ssa13*c);
        u_cb = (ca23*c - ssa23*b);
        u_bc = (ca23*b - ssa23*c);

        d1 = 2. * exp(-0.5*a*a) * normal_cdf(u_ba) * normal_cdf(u_ca) * bvnorm_cdf(r23, u_ba, u_ca) / (P * SQRT_2_PI);
        d2 = 2. * exp(-0.5*b*b) * normal_cdf(u_ab) * normal_cdf(u_cb) * bvnorm_cdf(r13, u_ab, u_cb) / (P * SQRT_2_PI);
        d3 = 2. * exp(-0.5*c*c) * normal_cdf(u_ac) * normal_cdf(u_bc) * bvnorm_cdf(r12, u_ac, u_bc) / (P * SQRT_2_PI);

At this point, I'm just trying to follow patterns, though, and my knowledge of derivatives of multivariate normal cdfs is somewhat limited.

Are there resources that give a more step-by-step derivation of the derivatives, other than general calculus principle resources? And, given what I've written above, does anyone know if I need a bivaraite PDF somewhere in there?

Any help would be appreciated. My hope is that a discussion about this might lead to this becoming a feature of gretl. I'm just trying to get BFGS to work first, since there's less to go wrong. If I provide the correct parameter estimates, the LL matches Limdep, so I think everything is okay up through the calculation of the triviariate CDF, and it *may* just be a matter of getting the right derivatives.

Thanks,
David