Re: st: Replicating Poi's quaids command

["Brian P. Poi" <bpoi@stata.com>]

On 09/04/2013 08:26 AM, Cruz Angel Echevarria wrote:

> 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;	
>
> /* 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);
>
> /* 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

Cruz,

I found two issues with your code, and one error in mine.  First,

> scalar b4 = 1 - b1 - b2 - b3;

should be

  scalar b4 = 0 - b1 - b2 - b3;

since the betas should sum to zero, not one.

Second, you have

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

This -nlsur- command is going to use the nlsurquaids.ado file that was shown in my 2008 Stata Journal article (assuming it was already installed).  However, that file has a_0 hardcoded as equal to 5.  So to get comparable results using -quaids- and your program, you must specify anot(5) and set your scalar a_0 = 5.

With those two changes to your code, we get much closer results.

The last error is in my code, and it only affects the computation of the expenditure elasticities.  There I forgot to include the factor of 0.5 that appears in the price index term.

I will get that fixed and in the "Software Update" section of the next issue of the Stata Journal.  That update will also have a new feature several users requested.  When you use the -estat- postestimation commands to get elasticities either at the variable means or else for a single observation in the dataset, those commands will now also return matrices containing the standard errors of the estimated elasticities.

   -- Brian Poi
      bpoi@stata.com

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