Re: [Gretl-devel] Trivariate Probit
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(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------