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 -------------------------------------------------------

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 -------------------------------------------------------