First page of Chapter 13 (User-defined functions), point 2 of the enumeration, "2. return-type, which states the type of value returned by the function, if any. This must be one of void (if the function does not return anything), scalar, series, matrix, list or string." But this also works for 'bundle', please add it to the list.

On Thu, 17 Jul 2014, David van Herick wrote: 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. > [...] It would be very nice to port Genz's code for the trivariate normal to libgretl. However, we've already got GHK for that and more, since you can use GHK for evaluating integrals of the multivariate normal in any (well, sort of) number of dimensions. Sure, it's a simulation-based method; however, there's a reason why it's become the "industry standard" in microeconometrics: it's really really REALLY good. As for the trivariate probit model, it's one of the examples I've been using to test our GHK implementation, and I could send you an example sample script if you want, or post it to the list if other people are interested. For your project, I'd suggest you take a look at the documentation for the GHK hansl function; even if you want to use C rather than hansl for your project, the corresponding functions you'll find in the C API are remarkably similar. ------------------------------------------------------- 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 ------------------------------------------------------- _______________________________________________ Gretl-devel mailing list Gretl-devel(a)lists.wfu.edu http://lists.wfu.edu/mailman/listinfo/gretl-devel

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

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