st: Replicating Poi's quaids command

["Cruz Angel Echevarria" <cruz.echevarria@ehu.es>]

To Whom It May Concern,

I would really appreciate some help on this. Running the following commands
to estimate a QUAIDS model

use food.dta, clear;
des;
quaids w1-w4, anot(10) prices(p1-p4) expenditure(expfd) nolog;
estat expenditure e*;
summarize e_1-e_4;

I obtain the following output for expenditure elasticities

    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
         e_1 |      4039    .1742147     1.21937  -60.39818   .7226771
         e_2 |      4029    .7332107    .2953817  -8.625433   .9495344
         e_3 |      3996    2.008486     1.26706   1.133542   36.11684
         e_4 |      4047    2.292034    2.938157   1.211376   169.9647

However, running my own code to estimate QUAIDS model, I obtain the same
parameter estimates that Boid's procedure does. I interpret this as that the
nlsurquaids program I have written is OK. However, when I try to compute the
expenditure elasticities, I obtain different results. This makes me think
that I am not translating the theoretical formulae properly into Stata code.

My code follows:

cscript 
#delimit; 
set more off;
version 13;
drop _all;
clear matrix;
capture log close;
log using output.log, replace;
set memory 500m, permanently; 
use food.dta, clear;
rename lnexp lnm;
des;
scalar a_0 = 10;

nlsur quaids @ w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnm,
  ifgnls neq(3) 
  param(a1 a2 a3 b1 b2 b3 g11 g12 g13 g22 g23 g33 l1 l2 l3)
  title(     Example QUAIDS)
  title2(    4 goods)
  nolog; 
 
matrix B = e(b)';
matrix list B;	  

/* a4 */	  
nlcom a4: 1 - (_b[/a1] + _b[/a2] + _b[/a3]);

/* b4 */
nlcom b4: 0 - (_b[/b1] + _b[/b2] + _b[/b3]);

/* l4 */
nlcom l4: 0 - (_b[/l1] + _b[/l2] + _b[/l3]);

/* g14 */
nlcom g14: 0 - (_b[/g11] + _b[/g12] + _b[/g13]);

/* g21 */
nlcom g21: _b[/g12];

/* g24 */
nlcom g24: 0 - (_b[/g12] + _b[/g22] + _b[/g23]);

/* g31 */
nlcom g31: _b[/g13];

/* g32 */
nlcom g32: _b[/g23];

/* g34 */
nlcom g34: 0 - (_b[/g13] + _b[/g23] + _b[/g33]);

/* g41 */
nlcom g41: 0 - (_b[/g11] + _b[/g12] + _b[/g13]);

/* g42 */
nlcom g42: 0 - (_b[/g12] + _b[/g22] + _b[/g23]);

/* g43 */
nlcom g43: 0 - (_b[/g13] + _b[/g23] + _b[/g33]);

/* g44 */
nlcom g44: 0 - (0 - (_b[/g11] + _b[/g12] + _b[/g13])) 
             - (0 - (_b[/g12] + _b[/g22] + _b[/g23])) 
			 - (0 - (_b[/g13] + _b[/g23] + _b[/g33]));


/* Scalars for parameter estimates */

scalar a1 = _b[/a1];
scalar a2 = _b[/a2];
scalar a3 = _b[/a3];
scalar a4 = 1 - a1 - a2 - a3;
scalar b1 = _b[/b1];
scalar b2 = _b[/b2];
scalar b3 = _b[/b3];
scalar b4 = 1 - b1 - b2 - b3;
scalar g11 = _b[/g11];
scalar g12 = _b[/g12];
scalar g13 = _b[/g13];
scalar g14 = 0 - (g11 + g12 + g13);
scalar g21 = g12;
scalar g22 = _b[/g22];
scalar g23 = _b[/g23];
scalar g24 = 0 - (g12 + g22 + g23);
scalar g31 = g13;
scalar g32 = g23;
scalar g33 = _b[/g33];
scalar g34 = 0 - (g13 + g23 + g33);
scalar g41 = g14;
scalar g42 = g24;
scalar g43 = g34;
scalar g44 = 0 - (g14 + g24 + g34);
scalar l1 = _b[/l1];
scalar l2 = _b[/l2];
scalar l3 = _b[/l3];
scalar l4 = 0 - (l1 + l2 + l3);

display a1;
display a2;
display a3;
display a4;
display b1;
display b2;
display b3;
display b4;
display l1;
display l2;
display l3;
display l4;
display g11;
display g12;
display g13;
display g14;
display g21;
display g22;
display g23;
display g24;
display g31;
display g32;
display g33;
display g34;
display g41;
display g42;
display g43;
display g44;

/* This produces the same estimates tan those obtained 
by running the aforementioned quaids command */
/* ===== */
/* FOR ELASTICITIES                      */
/* ===== */

/* Preliminaries */

gen lnb = b1*lnp1 + b2*lnp2 + b3*lnp3 + b4*lnp4;

gen bindex = exp(lnb);

gen lna = a_0 + a1*lnp1 + a2*lnp2 + a3*lnp3 + a4*lnp4 +
         (1/2)*(
   g11*lnp1*lnp1 + g12*lnp1*lnp2 + g13*lnp1*lnp3 + g14*lnp1*lnp4
 + g21*lnp2*lnp1 + g22*lnp2*lnp2 + g23*lnp2*lnp3 + g24*lnp2*lnp4
 + g31*lnp3*lnp1 + g32*lnp3*lnp2 + g33*lnp3*lnp3 + g34*lnp3*lnp4
 + g41*lnp4*lnp1 + g42*lnp4*lnp2 + g43*lnp4*lnp3 + g44*lnp4*lnp4
             );
 
gen aindex = exp(lna);


/* INCOME ELASTICITIES */
gen mu1 = b1 + ((l1*2)/bindex)*(lnm - lna);
gen mu2 = b2 + ((l2*2)/bindex)*(lnm - lna);
gen mu3 = b3 + ((l3*2)/bindex)*(lnm - lna);
gen mu4 = b4 + ((l4*2)/bindex)*(lnm - lna);
gen em1 = (mu1/w1) + 1; 
gen em2 = (mu2/w2) + 1; 	
gen em3 = (mu3/w3) + 1; 	
gen em4 = (mu4/w4) + 1; 

sum em*;

And this is the output that I obtain


    Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
         em1 |      4039    1.556595    1.060537  -1.413979   37.93764
         em2 |      4029    .9256395    .1425985  -2.658406   2.997787
         em3 |      3996    .4810829    1.376963  -46.75089   4.259475
         em4 |      4047    5.482535    8.378185  -.6769157   437.3061


which clearly differs from the one above shown.


I would really appreciate any help on this. 
Regards, 
Cruz

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/