Re: st: Error in Poi's QUAIDS user-command

[Jorge Eduardo Pérez Pérez <jorge_perez@brown.edu>]

Which version of -quaids- are you using? Mine is version 1.1.0 and the
code reads:

if (quadratics == "") {
bofp = exp(lnp*beta')
lnpindex = lnp*alpha' :+ st_numscalar("e(anot)")
for(i=1; i<=rows(lnpindex); ++i) {
lnpindex[i] = lnpindex[i] + 0.5*lnp[i,.]*gamma*lnp[i,.]'
}
}


--------------------------------------------
Jorge Eduardo Pérez Pérez
Graduate Student
Department of Economics
Brown University


On Fri, Apr 18, 2014 at 2:43 PM, Robin Winkler <robinwinkler@gmail.com> wrote:
> Hi all,
>
> I believe I have found a serious error in the way in which the QUAIDS
> user command written by Brian Poi computes elasticities.
>
> In the Mata code underlying the ado file - the relevant bit is pasted
> below - the price index *lnpindex* used to compute the elasticities
> does not include alpha 0, which is entered into the command suite as
> *anot*. Hence, I have had the problem that  my price index was
> underestimated, and the magnitude of the computed expenditure
> elasticities tend to be way overstated.
>
> Should Mr Poi read this, I would be grateful for a comment.
>
> Robin
>
>
>
> *! version 1.0.0  29dec2011
> /*
> _quaids__utils.mata
> Mata routines called by
> nlsur__quaids.ado
> quaids_estat.ado
> quaids_p.ado
> */
>
> mata
>
> mata set matastrict on
>
> void _quaids__expelas(string scalar touses,
>      string scalar quadratics,
>      string scalar atmeans,
>      string scalar lnps,
>      string scalar lnexps,
>      real   scalar ndemo,
>      string scalar demos,
>      string scalar expelass)
> {
> real scalar i
> real vector alpha, beta, lambda, rho
> real vector bofp, cofp, lnpindex, lnexp, mbar
> real matrix gamma, eta
> real matrix shares, expelas, lnp, demo
> alpha  = st_matrix("e(alpha)")
> beta   = st_matrix("e(beta)")
> gamma  = st_matrix("e(gamma)")
> if (quadratics == "") {
> lambda = st_matrix("e(lambda)")
> }
> else {
> lambda = J(1, cols(beta), 0)
> }
> if (ndemo > 0) {
> eta = st_matrix("e(eta)")
> rho = st_matrix("e(rho)")
> if (atmeans == "") {
> st_view(demo=., ., demos, touses)
> }
> else {
> demo = mean(st_data(., demos, touses))
> }
> }
>
> if (atmeans == "") {
> st_view(shares=.,  ., st_global("e(lhs)"),  touses)
> st_view(expelas=., ., expelass, touses)
> st_view(lnp=.,     ., lnps, touses)
> st_view(lnexp=.,   ., lnexps, touses)
> }
> else {
> shares  = mean(st_data(., st_global("e(lhs)"), touses))
> lnp     = mean(st_data(., lnps, touses))
> lnexp   = mean(st_data(., lnexps, touses))
> expelas = J(1, cols(shares), .)
> }
> if (quadratics == "") {
> bofp = exp(lnp*beta')
> lnpindex = lnp*alpha'
> for(i=1; i<=rows(lnpindex); ++i) {
> lnpindex[i] = lnpindex[i] + lnp[i,.]*gamma*lnp[i,.]'
> }
> }
> else {
> bofp = J(rows(lnp), 1, 1) // 1, so now div0 problem
> lnpindex = J(rows(lnp), 1, 0)
> }
> if (ndemo > 0) {
> cofp = J(rows(lnp), 1, 0)
> for(i=1; i<=rows(lnp); ++i) {
> cofp[i] = lnp[i,.]*(eta'*demo[i,.]')
> }
> cofp = exp(cofp)
> mbar = 1 :+ demo*rho'
> for(i=1; i<=rows(expelas); ++i) {
> expelas[i,.] = 1 :+ 1:/shares[i,.]:*
> (beta + demo[i,.]*eta + 2*lambda/bofp[i]/cofp[i]*
> (lnexp[i]-ln(mbar[i])-lnpindex[i]))
> }
> }
> else {
> for(i=1; i<=rows(expelas); ++i) {
> expelas[i,.] = 1 :+ 1:/shares[i,.]:*
> (beta + 2*lambda/bofp[i]*(lnexp[i]-lnpindex[i]))
> }
> }
> if (atmeans != "") {
> st_matrix(expelass, expelas)
> *
> *   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/

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