st: Error in Poi's QUAIDS user-command

[Robin Winkler <robinwinkler@gmail.com>]

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/