On Fri, 18 Jul 2014, David van Herick wrote:

Thanks.  I will look at the documentation on the GHK hansl function.  Would
you please send (or post) the example script for the trivariate probit
model?

Here you are. This uses numerical derivatives. It could also be done with analytical derivatives, but the script would be a bit longer.

<hansl>
function series trivp_ll(matrix param, series y1, series y2, series y3,
                         list X1, list X2, list X3, matrix *U)
#    printf "%8.3f", param'

    set warnings off

    k1 = nelem(X1)
    k2 = nelem(X2)
    k3 = nelem(X3)

    cf1 = param[1:k1]
    cf2 = param[k1+1:k1+k2]
    cf3 = param[k1+k2+1:k1+k2+k3]

    c = param[k1+k2+k3+1:k1+k2+k3+3]
    if maxc(abs(c)) > 1
        series ret = NA
        return ret
    endif

    S = zeros(3,3)
    S[2,1] = c[1]
    S[3,1] = c[2]
    S[3,2] = c[3]

    c = sumr(S.^2)
    if maxc(c) > 1
        series ret = NA
        return ret
    endif

    S[diag] = 1 - c

    loop foreach j 1 2 3 -q
        series ndx = -lincomb(X$j, cf$j)
        series t$j = y$j ? $huge : ndx
        series b$j = y$j ? ndx : -$huge
    endloop

    matrix A = {b1, b2, b3}
    matrix B = {t1, t2, t3}
    ret = ghk(S, A, B, U)

    return ln(ret)

end function

# --- main ----------------------------------------------

set echo off
set messages off
set skip_missing off
set seed 732237

# generate some artificial data

nulldata 1200

k1 = 1
k2 = 2
k3 = 3

C12 = 0.3
C13 = 0.5
C23 = 0.6

Sigma = {1; C12; C13; 1; C23; 1}
Sigma = unvech(Sigma)

print Sigma

list X = const
loop i=1..k1 -q
    x_$i = uniform()
    X += x_$i
endloop
b1 = mnormal(k1+1, 1)

list Z = const
loop i=1..k2 -q
    z_$i = uniform()
    Z += z_$i
endloop
b2 = mnormal(k2+1, 1)

list W = const
loop i=1..k3 -q
    w_$i = uniform()
    W += w_$i
endloop
b3 = mnormal(k3+1, 1)

matrix E = mnormal($nobs, 3) * cholesky(Sigma)'

series y1star = lincomb(X, b1)+E[,1]
y1 = (y1star>0)
series y2star = lincomb(Z, b2)+E[,2]
y2 = (y2star>0)
series y3star = lincomb(W, b3)+E[,3]
y3 = (y3star>0)

# store prova.csv y1 y2 y3 x* z* w* --csv

# --- proceed with estimation ------------------

rep = 100
U = muniform(3, rep)

# initialise the params by univariate probit models
probit y1 X -q
b1 = $coeff
probit y2 Z -q
b2 = $coeff
probit y3 W -q
b3 = $coeff

matrix theta = b1 | b2 | b3 | zeros(3,1)

# go with MLE
set stopwatch
mle ll = trivp_ll(theta, y1, y2, y3, X, Z, W, &U)
    params theta
end mle -v
t0 = $stopwatch
printf "Time = %g\n", t0
</hansl>

I know Limdep also uses a GHK simulator for its MPROBIT algorithm.
However, the problem is that even a simple trivariate model can take quite
a while to estimate with Limdep's default of 100 pts (on the order or 10-20
minutes for the *N*=750 test model I'm running).  Is gretl's GHK simulator
any faster, or are slow speeds just the nature of the GHK simulation?

Of course, that depends on the number of covariates you're using. Our GHK machine is quite good, if I may say so myself, both in speed and accuracy. Please let us know your impressions.

I think most of Genz's n>3 algorithms are simulation based also, but he has functions that estimate the special trivariate case much more quickly.  I think if those were implemented, the special case of the trivariate probit would run much faster than using GHK.  However, I'm sure someone must have thought of this before, so I'm a little wary that I'm missing something.

True. As I said, it'd be very nice to add Genz's algorithm for the trivariate normal to our rangen() weaponry. We're going to release a new version very soon, so we have all the time to work on this after release.

It's not GHK related, but is there a way run to run a simple (two-level)
FIML nested logit model in gretl?

Not natively, but a Hansl script shouldn't be hard to write. If I'm not mistaken, Claudia Pigini could already have something ready, but I may be wrong. Claudia?


-------------------------------------------------------
  Riccardo (Jack) Lucchetti
  Dipartimento di Scienze Economiche e Sociali (DiSES)

  Università Politecnica delle Marche
  (formerly known as Università di Ancona)

  r.lucchetti@univpm.it
  http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------