[tareq]

Statalisters,

After ample attempts at constructing a nlsur quaids function evaluator program to estimate a QUAIDS demand model with demographic scaling and censoring and even more browsing through the statalist, I have hit a brick wall, particularly with the demographic scaling. I would like to apply the demographic scaling approach applied by Poi (in Stata Journal (2012) 12, No. 3) but am not able to estimate the parameter rho (ρ) needed to calculate the mbar function which is necessary to scale total expenditure.

The main question is how can one arrive at rho within the framework of nlsur? Having studied the ado files for Stata's -quaids command I can see it produces a 1xn matrix (where n = number of demographic variables included) for rho, yet my mata skills are not strong enough to figure out how this matrix is populated. It should be noted that this is different from the n x g matrix for the eta's (where g is the number of goods). Unfortunately all of the examples included in this string of posts estimate the eta's yet ignore rho and none have actually taken the steps necessary to estimate the cofp and mbar functions needed for the scaling.

I would appreciate any help or suggestions you may provide.

Many thanks.

[Debapriya Chakraborty]

Hi Statalisters,

I am modifying Nigussie's program quaidsNNP for the purpose of my 3 goods model but I keep encountering the following error message:
nlsurquaidsdc1 returned 103
verify that nlsurquaidsdc1 is a function evaluator program

Model Description: 3 electric vehicles charging period
y-variable: expenditure share in period 'i 'of day 'd' where i=1,2,3 and d= date of treatment
x-variables: lnp1, lnp2, lnp3, lnexp, pdf1, pdf2,pdf3, cdf1, cdf2, cdf3. I do not have any demographic information. In the first stage, I estimated a mixed probit model with random effects (meprobit) to obtain the pdf and CDF estimates.

Now, I am working on the nlsur quaids model to obtain the own- and cross-price elasticities. For this purpose, I want to modify the nlsurquaidsNNP program. But I keep getting the error message. This is the first time I am modifying a STATA program. I tried debugging with trace on but I could not identify the issue. I was hoping to get some help from you all Example.xlsx .

Please find attached my code for the nlsur program and a snapshot of the data .

*! version 1.0.1 08may2008
program nlsurquaidsdc1

version 14
syntax varlist(min=8 max=8) if, at(name)
tokenize `varlist'
args w1 w2 lnp1 lnp2 lnp3 lnexp pdf1 pdf2 cdf1 cdf2
// With three periods, there are 8 parameters that can be
// estimated, after eliminating one of the periods and
// imposing adding up, symmetry, and homogeneity
// constraints, in the QUAIDS model
// Here, we extract those parameters from the `at'
// vector, and impose constraints as we go along

tempname a1 a2 a3
scalar `a1' = `at'[1,1]
scalar `a2' = `at'[1,2]
scalar `a3' = 1 - `a1' - `a2'
tempname b1 b2 b3
scalar `b1' = `at'[1,3]
scalar `b2' = `at'[1,4]
scalar `b3' = -`b1' - `b2'
tempname g11 g12 g13
tempname g21 g22 g23
tempname g31 g32 g33
scalar `g11' = `at'[1,5]
scalar `g12' = `at'[1,6]
scalar `g13' = -`g11' - `g12'
scalar `g21' = `g12'
scalar `g22' = `at'[1,7]
scalar `g23' = -`g21' - `g22'
scalar `g31' = `g13'
scalar `g32' = `g23'
scalar `g33' = -`g31' - `g32'
tempname l1 l2 l3
scalar `l1' = `at'[1,8]
scalar `l2' = `at'[1,9]
scalar `l3' = -`l1' - `l2'
tempname d1 d2
scalar `d1' = `at'[1,10]
scalar `d2' = `at'[1,11]
// Okay, now that we have all the parameters, we can
// calculate the expenditure shares.
quietly {
// First get the price index
// I set a_0 = 5
tempvar lnpindex
gen double `lnpindex' = 5 + `a1'*`lnp1' + `a2'*`lnp2' ///
+ `a3'*`lnp3'
forvalues i = 1/3 {
forvalues j = 1/3 {
replace `lnpindex' = `lnpindex' + ///
0.5*`g`i'`j''*`lnp`i''*`lnp`j''
}
}
// The b(p) term in the QUAIDS model:
tempvar bofp
gen double `bofp' = 0
forvalues i = 1/3 {
replace `bofp' = `bofp' + `lnp`i''*`b`i''
}
replace `bofp' = exp(`bofp')
// Finally, the expenditure shares for 2 of the 3
// periods (the third is dropped to avoid singularity)
replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' + ///
`g13'*`lnp3' + ///
`b1'*(`lnexp' - `lnpindex') + ///
`l1'/`bofp'*(`lnexp' - `lnpindex')^2)*`cdf1'+ ///
`d1'*`pdf1'
replace `w2' = (`a2' + `g21'*`lnp1' + `g22'*`lnp2' + ///
`g23'*`lnp3' + ///
`b2'*(`lnexp' - `lnpindex') + ///
`l2'/`bofp'*(`lnexp' - `lnpindex')^2)*`cdf2' + ///
`d2'*`pdf2'
}
end
I am aware of the QUAIDS model in STATA but I need to do the 2 step correction to account for zero expenditure. There are no missing values in the variables included in the model.
Hoping to get some feedback on how to debug the program. Thanks in advance.

[Temitope Oluwole]

Hi Statalisters

I am trying to estimate a 4 goods censored QUAIDS model and I am using the Shonkwiler and Yen's method (1999) to correct for the zero expenditure reported by respondents.

Also, I am using the “nlsur” command to estimate this censored QUAIDS model. However, I get the following error message: [d1] not found

Any assitance or suggestions will be appreciated.

Thank you very much in advance.

Oluwole Temitope.


Here are the Stata codes:

cap program drop nlsurquaidsTPP
program nlsurquaidsTPP

version 14.0

syntax varlist(min=16 max=16) if, at(name)

tokenize `varlist'
args w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnexp z1 z2 pdf1 pdf2 pdf3 cdf1 cdf2 cdf3

// With four goods, there are 15 parameters that can be
// estimated, after eliminating one of the goods and
// imposing adding up, symmetry, and homogeneity
// constraints, in the QUAIDS model
// Here, we extract those parameters from the `at'
// vector, and impose constraints as we go along

tempname a1 a2 a3 a4
scalar `a1' = `at'[1,1]
scalar `a2' = `at'[1,2]
scalar `a3' = `at'[1,3]
scalar `a4' = 1 - `a1' - `a2' - `a3'

tempname b1 b2 b3 b4
scalar `b1' = `at'[1,4]
scalar `b2' = `at'[1,5]
scalar `b3' = `at'[1,6]
scalar `b4' = -`b1' - `b2' - `b3'

tempname g11 g12 g13 g14
tempname g21 g22 g23 g24
tempname g31 g32 g33 g34
tempname g41 g42 g43 g44
scalar `g11' = `at'[1,7]
scalar `g12' = `at'[1,8]
scalar `g13' = `at'[1,9]
scalar `g14' = -`g11' - `g12' - `g13'

scalar `g21' = `g12'
scalar `g22' = `at'[1,10]
scalar `g23' = `at'[1,11]
scalar `g24' = -`g21' - `g22' - `g23'

scalar `g31' = `g13'
scalar `g32' = `g23'
scalar `g33' = `at'[1,12]
scalar `g34' = -`g31' - `g32' - `g33'

scalar `g41' = `g14'
scalar `g42' = `g24'
scalar `g43' = `g34'
scalar `g44' = -`g41' - `g42' - `g43'

tempname l1 l2 l3 l4
scalar `l1' = `at'[1,13]
scalar `l2' = `at'[1,14]
scalar `l3' = `at'[1,15]
scalar `l4' = -`l1' - `l2' - `l3'

// constant and household demographics
tempname r11 r12
tempname r21 r22
tempname r31 r32

scalar `r11' = `at'[1,16]
scalar `r12' = `at'[1,17]
scalar `r21' = `at'[1,18]
scalar `r22' = `at'[1,19]
scalar `r31' = `at'[1,20]
scalar `r32' = `at'[1,21]

// pdf
tempname d1 d2 d3
scalar `d1' = `at'[1,22]
scalar `d2' = `at'[1,23]
scalar `d3' = `at'[1,24]

// Okay, now that we have all the parameters, we can
// calculate the expenditure shares.
quietly {
// First get the price index
// I set a_0 = 5
tempvar lnpindex
gen double `lnpindex' = 5 + `a1'*`lnp1' + `a2'*`lnp2' ///
+ `a3'*`lnp3' + `a4'*`lnp4'
forvalues i = 1/4 {
forvalues j = 1/4 {
replace `lnpindex' = `lnpindex' + ///
0.5*`g`i'`j''*`lnp`i''*`lnp`j''
}
}
// The b(p) term in the QUAIDS model:
tempvar bofp
gen double `bofp' = 0
forvalues i = 1/4 {
replace `bofp' = `bofp' + `lnp`i''*`b`i''
}
replace `bofp' = exp(`bofp')
// Finally, the expenditure shares for 3 of the 4
// goods (the fourth is dropped to avoid singularity)
replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' + ///
`g13'*`lnp3' + `g14'*`lnp4' + ///
`b1'*(`lnexp' - `lnpindex') + ///
`l1'/`bofp'*(`lnexp' - `lnpindex')^2 + ///
`r11'*`z1' +`r12'*`z2')*`cdf1'+ ///
`d1'*`pdf1'

replace `w2' = (`a2' + `g21'*`lnp1' + `g22'*`lnp2' + ///
`g23'*`lnp3' + `g24'*`lnp4' + ///
`b2'*(`lnexp' - `lnpindex') + ///
`l2'/`bofp'*(`lnexp' - `lnpindex')^2 + ///
`r21'*`z1' +`r22'*`z2')*`cdf2' + ///
`d2'*`pdf2'

replace `w3' = (`a3' + `g31'*`lnp1' + `g32'*`lnp2' + ///
`g33'*`lnp3' + `g34'*`lnp4' + ///
`b3'*(`lnexp' - `lnpindex') + ///
`l3'/`bofp'*(`lnexp' - `lnpindex')^2 + ///
`r31'*`z1' +`r32'*`z2')*`cdf3' + ///
`d3'*`pdf3'
}

end



****************************************


nlsur quaidsTPP @ w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnexp ///
z1 z2 pdf1 pdf2 pdf3 cdf1 cdf2 cdf3, ifgnls nequations(3) ///
param(a1 a2 a3 ///
g11 g21 g31 ///
g22 g32 ///
g33 ///
b1 b2 b3 ///
l1 l2 l3 ///
r1 r2 ///
m11 m12 m13 ///
m21 m22 m23 ///
d1 d2 d3 ) nolog

est store quaidsTPP2

*set trace on
*set tracedepth 4

*Share, price and total expenditure means
quietly {
foreach x of varlist w* lnp* lnexp {
sum `x'
scalar `x'mean=r(mean)
}
* Price indexes
glo asum "_b[a1]*lnp1mean"
forv i=2(1)4 {
glo asum "${asum} + _b[a`i']*lnp`i'mean"
}
glo gsum ""
forv i=1(1)4 {
forv j=1(1)4 {
glo gsum "${gsum} + 0.5*_b[g`i'`j']*lnp`i'mean*lnp`j'mean"
}
}
glo ap "10.0 + ${asum} ${gsum}"
glo bp "_b[b1]*lnp1mean"
forv i=2(1)4 {
glo bp "${bp} + _b[b`i']*lnp`i'mean"
}
glo bp "(exp(${bp}))"
* Mus
forv i=1(1)4 {
glo mu`i' "_b[b`i'] + 2*_b[l`i']/${bp}*(lnexpmean-(${ap}))"
}
forv j=1(1)4 {
glo gsum2`j' ""
forv k=1(1)4 {
glo gsum2`j' "${gsum2`j'} + _b[g`j'`k']*lnp`k'mean"
}
}
}
nlcom (a1:_b[/a1])(a2:_b[/a2])(a3:_b[/a3])(a4:1-_b[/a1]-_b[/a2]-_b[/a3]) ///
(b1:_b[/b1])(b2:_b[/b2])(b3:_b[/b3])(b4:-_b[/b1]-_b[/b2]-_b[/b3]) ///
(g11:_b[/g11])(g12:_b[/g21])(g13:_b[/g31]) ///
(g21:_b[/g21])(g22:_b[/g22])(g23:_b[/g32]) ///
(g31:_b[/g31])(g32:_b[/g32])(g33:_b[/g33]) ///
(g14:-_b[/g11]-_b[/g21]-_b[/g31]) ///
(g24:-_b[/g21]-_b[/g22]-_b[/g32]) ///
(g34:-_b[/g31]-_b[/g32]-_b[/g33]) ///
(g41:-_b[/g11]-_b[/g21]-_b[/g31]) ///
(g42:-_b[/g21]-_b[/g22]-_b[/g32]) ///
(g43:-_b[/g31]-_b[/g32]-_b[/g33]) ///
(g44:-(-_b[/g11]-_b[/g21]-_b[/g31])-(-_b[/g21]-_b[/g22]-_b[/g32])-(-_b[/g31]-_b[/g32]-_b[/g33])) ///
(l1:_b[/l1])(l2:_b[/l2])(l3:_b[/l3])(l4:-_b[/l1]-_b[/l2]-_b[/l3]) ///
(d1:_b[/d1]) (d2:_b[/d2]) (d3:_b[/d3]), post noheader


* est store quaidsTPP2


forv i=1(1)4 {
forv j=1(1)4 {
glo delta=cond(`i'==`j',1,0)
glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j'] ${gsum2`j'})-_b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2"
cap nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}), post noheader
if _rc {
qui nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f: (1e+2)*(${mu`i'`j'})), post noheader
qui nlcom (elasexp`i': _b[elasexp`i']) (mu`i'`j':_b[mu`i'`j'f]/(1e+2)), post noheader
}
* Uncompensated price elasticity
nlcom (elasexp`i': _b[elasexp`i']*_b[d`i']) (elu`i'`j':(_b[mu`i'`j']/w`i'mean - ${delta})*_b[d`i']) , post noheader
* Compensated price elasticity
nlcom (elc`i'`j': _b[elu`i'`j'] + _b[elasexp`i']*w`j'mean), noheader
qui est restore quaidsTPP2
}
}

This is the place I got the error message:

* est store quaidsTPP2
.
.
. forv i=1(1)4 {
  2. forv j=1(1)4 {
  3. glo delta=cond(`i'==`j',1,0)
  4. glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j'] ${gsum2`j'})-_b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2"
  5. cap nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}), post noheader
  6. if _rc {
  7. qui nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f: (1e+2)*(${mu`i'`j'})), post noheader
  8. qui nlcom (elasexp`i': _b[elasexp`i']) (mu`i'`j':_b[mu`i'`j'f]/(1e+2)), post noheader
  9. }
 10. * Uncompensated price elasticity
. nlcom (elasexp`i': _b[elasexp`i']*_b[d`i']) (elu`i'`j':(_b[mu`i'`j']/w`i'mean - ${delta})*_b[d`i']) , post noheader
 11. * Compensated price elasticity
. nlcom (elc`i'`j': _b[elu`i'`j'] + _b[elasexp`i']*w`j'mean), noheader
 12. qui est restore quaidsTPP2
 13. }
 14. }
[d1] not found
r(111);

end of do-file

r(111);

[Jorge Eduardo Perez Perez]

You need to restore the original quaids estimates before you calculate the price elasticity. This is because your previous nlcom calculations -post- their results, so d1 from the original estimates is not posted when you try to calculate the price elasticity.

[Temitope Oluwole]

Thank you for your assistance, Jorge. But, how do I restore the original quaids and where do I put the code for restore.

[Srabashi Ray]

Hello everyone! Thanks for sharing codes to solve the censoring problem in QAIDS.
I tried to run the 4 good code as in post #16 posted on 17th Sept 2015 but got the following error :

[b1] not found
r(111);

while the last 'forv' loop was running.

Also, do these codes have calculations for Income elasticity?

Thank you for any help!!
Just to clarify: I did not make any changes to the code I ran from post #16. Just copied it from here and ran it on Stata.

[Srabashi Ray]

Hi ... I was able to run codes in post#16 by allowing line 'est store quaidsNNP2' to run just before the elasticity estimations in the end. Is that the right way to make it run?

However, when I compare these results with the direct QAIDS estimation using Poi's 2012 QUAIDS command with the demographics option for nkids and rural, I do not get the same results.

The QAIDS code I use to test this is - quaids w1-w4, anot(10) prices(p1-p4) expenditure(expfd) nolog

A closer look at the codes in post #16, it seems like they do not define w1, w2,... correctly when demographics is used.
The part where w1, w2, w3... is defined - it seems to me that they are defined as the usual functional form for QAIDS as in Equation 10 in the Banks et. al 1997 paper. However if demographics are included (and I think x1 and x2 demographic variables here), the w1, w2, w3... equations should be re-written as in Poi's 2012 paper on pp. 436. Is that right? It could be great if I could hear someone else's thoughts on this before I start playing around with these codes!

Hope to hear from someone!
Regards,
Srabashi

[Moses Muwonge]

Hi Srabashi,

Please i have come across you recent post but this is a similar problem am experiencing in estimation of a 6 good Quaids model but so far i have got a solution . Some help on how you have managed to deal with this sort of problem.

Thank you in advance
Regards
Moses

[Ahalya Ramanathan]

Dear David,
I have been trying to run the codes for quaids model including demographics and censoring with 13 goods and 10 demographic variables. Although I have run the model, but I am unable to use these codes to calculate elasticities as it shows error. Kindly help me with the codes for this model.

[Ahalya Ramanathan]

Hello,
I have run the following code for a quaids model with demographics (x1, x2,....) and censoring comprising of 13 goods. Although the model is running perfectly, I am not able to write the codes for calculating elasticity. Kindly help.

clear
use "C:\Users\Ramanathan\Documents\Demand system estimation\Food_components_data_rural"

cap program drop nlsurquaids
program nlsurquaids

version 12.1

syntax varlist(min=63 max=63) if, at(name)
tokenize `varlist'
args w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 w11 w12 w13 lnp1 lnp2 lnp3 lnp4 lnp5 lnp6 lnp7 lnp8 lnp9 lnp10 lnp11 lnp12 lnp13 lnexp x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 pdf1 pdf2 pdf3 pdf4 pdf5 pdf6 pdf7 pdf8 pdf9 pdf10 pdf11 pdf12 pdf13 cdf1 cdf2 cdf3 cdf4 cdf5 cdf6 cdf7 cdf8 cdf9 cdf10 cdf11 cdf12 cdf13
// With 13 goods, there are 114 parameters that can be estimated

*
**parameters are best viewed as a single vector
*
* first 13 alphas
*
tempname a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13
scalar `a1' = `at'[1,1]
scalar `a2' = `at'[1,2]
scalar `a3' = `at'[1,3]
scalar `a4' = `at'[1,4]
scalar `a5' = `at'[1,5]
scalar `a6' = `at'[1,6]
scalar `a7' = `at'[1,7]
scalar `a8' = `at'[1,8]
scalar `a9' = `at'[1,9]
scalar `a10' = `at'[1,10]
scalar `a11' = `at'[1,11]
scalar `a12' = `at'[1,12]
scalar `a13' = `at'[1,13]

*
* first 13 betas
*
tempname b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13
scalar `b1' = `at'[1,14]
scalar `b2' = `at'[1,15]
scalar `b3' = `at'[1,16]
scalar `b4' = `at'[1,17]
scalar `b5' = `at'[1,18]
scalar `b6' = `at'[1,19]
scalar `b7' = `at'[1,20]
scalar `b8' = `at'[1,21]
scalar `b9' = `at'[1,22]
scalar `b10' = `at'[1,23]
scalar `b11' = `at'[1,24]
scalar `b12' = `at'[1,25]
scalar `b13' = `at'[1,26]

*
* vector gamma star k rows and columns of gamma
*

tempname g11 g12 g13 g14 g15 g16 g17 g18 g19 g110 g111 g112 g113
tempname g21 g22 g23 g24 g25 g26 g27 g28 g29 g210 g211 g212 g213
tempname g31 g32 g33 g34 g35 g36 g37 g38 g39 g310 g311 g312 g313
tempname g41 g42 g43 g44 g45 g46 g47 g48 g49 g410 g411 g412 g413
tempname g51 g52 g53 g54 g55 g56 g57 g58 g59 g510 g511 g512 g513
tempname g61 g62 g63 g64 g65 g66 g67 g68 g69 g610 g611 g612 g613
tempname g71 g72 g73 g74 g75 g76 g77 g78 g79 g710 g711 g712 g713
tempname g81 g82 g83 g84 g85 g86 g87 g88 g89 g810 g811 g812 g813
tempname g91 g92 g93 g94 g95 g96 g97 g98 g99 g910 g911 g912 g913
tempname g101 g102 g103 g104 g105 g106 g107 g108 g109 g1010 g1011 g1012 g1013
tempname g111 g112 g113 g114 g115 g116 g117 g118 g119 g1110 g1111 g1112 g1113
tempname g121 g122 g123 g124 g125 g126 g127 g128 g129 g1210 g1211 g1212 g1213
tempname g131 g132 g133 g134 g135 g136 g137 g138 g139 g1310 g1311 g1312 g1313


scalar `g11' = `at'[1,27]
scalar `g12' = `at'[1,28]
scalar `g13' = `at'[1,29]
scalar `g14' = `at'[1,30]
scalar `g15' = `at'[1,31]
scalar `g16' = `at'[1,32]
scalar `g17' = `at'[1,33]
scalar `g18' = `at'[1,34]
scalar `g19' = `at'[1,35]
scalar `g110' = `at'[1,36]
scalar `g111' = `at'[1,37]
scalar `g112' = `at'[1,38]
scalar `g113' = `at'[1,39]


scalar `g21' = `g12'
scalar `g22' = `at'[1,40]
scalar `g23' = `at'[1,41]
scalar `g24' = `at'[1,42]
scalar `g25' = `at'[1,43]
scalar `g26' = `at'[1,44]
scalar `g27' = `at'[1,45]
scalar `g28' = `at'[1,46]
scalar `g29' = `at'[1,47]
scalar `g210' = `at'[1,48]
scalar `g211' = `at'[1,49]
scalar `g212' = `at'[1,50]
scalar `g213' = `at'[1,51]



scalar `g31' = `g13'
scalar `g32' = `g23'
scalar `g33' = `at'[1,52]
scalar `g34' = `at'[1,53]
scalar `g35' = `at'[1,54]
scalar `g36' = `at'[1,55]
scalar `g37' = `at'[1,56]
scalar `g38' = `at'[1,57]
scalar `g39' = `at'[1,58]
scalar `g310' = `at'[1,59]
scalar `g311' = `at'[1,60]
scalar `g312' = `at'[1,61]
scalar `g313' = `at'[1,62]



scalar `g41' = `g14'
scalar `g42' = `g24'
scalar `g43' = `g34'
scalar `g44' = `at'[1,63]
scalar `g45' = `at'[1,64]
scalar `g46' = `at'[1,65]
scalar `g47' = `at'[1,66]
scalar `g48' = `at'[1,67]
scalar `g49' = `at'[1,68]
scalar `g410' = `at'[1,69]
scalar `g411' = `at'[1,70]
scalar `g412' = `at'[1,71]
scalar `g413' = `at'[1,72]


scalar `g51' = `g15'
scalar `g52' = `g25'
scalar `g53' = `g35'
scalar `g54' = `g45'
scalar `g55' = `at'[1,73]
scalar `g56' = `at'[1,74]
scalar `g57' = `at'[1,75]
scalar `g58' = `at'[1,76]
scalar `g59' = `at'[1,77]
scalar `g510' = `at'[1,78]
scalar `g511' = `at'[1,79]
scalar `g512' = `at'[1,80]
scalar `g513' = `at'[1,81]


scalar `g61' = `g16'
scalar `g62' = `g26'
scalar `g63' = `g36'
scalar `g64' = `g46'
scalar `g65' = `g56'
scalar `g66' = `at'[1,82]
scalar `g67' = `at'[1,83]
scalar `g68' = `at'[1,84]
scalar `g69' = `at'[1,85]
scalar `g610' = `at'[1,86]
scalar `g611' = `at'[1,87]
scalar `g612' = `at'[1,88]
scalar `g613' = `at'[1,89]


scalar `g71' = `g17'
scalar `g72' = `g27'
scalar `g73' = `g37'
scalar `g74' = `g47'
scalar `g75' = `g57'
scalar `g76' = `g67'
scalar `g77' = `at'[1,90]
scalar `g78' = `at'[1,91]
scalar `g79' = `at'[1,92]
scalar `g710' = `at'[1,93]
scalar `g711' = `at'[1,94]
scalar `g712' = `at'[1,95]
scalar `g713' = `at'[1,96]

scalar `g81' = `g18'
scalar `g82' = `g28'
scalar `g83' = `g38'
scalar `g84' = `g48'
scalar `g85' = `g58'
scalar `g86' = `g68'
scalar `g87' = `g78'
scalar `g88' = `at'[1,97]
scalar `g89' = `at'[1,98]
scalar `g810' = `at'[1,99]
scalar `g811' = `at'[1,100]
scalar `g812' = `at'[1,101]
scalar `g813' = `at'[1,102]


scalar `g91' = `g19'
scalar `g92' = `g29'
scalar `g93' = `g39'
scalar `g94' = `g49'
scalar `g95' = `g59'
scalar `g96' = `g69'
scalar `g97' = `g79'
scalar `g98' = `g89'
scalar `g99' = `at'[1,103]
scalar `g910' = `at'[1,104]
scalar `g911' = `at'[1,105]
scalar `g912' = `at'[1,106]
scalar `g913' = `at'[1,107]



scalar `g101' = `g110'
scalar `g102' = `g210'
scalar `g103' = `g310'
scalar `g104' = `g410'
scalar `g105' = `g510'
scalar `g106' = `g610'
scalar `g107' = `g710'
scalar `g108' = `g810'
scalar `g109' = `g910'
scalar `g1010' = `at'[1,108]
scalar `g1011' = `at'[1,109]
scalar `g1012' = `at'[1,110]
scalar `g1013' = `at'[1,111]


scalar `g111' = `g111'
scalar `g112' = `g211'
scalar `g113' = `g311'
scalar `g114' = `g411'
scalar `g115' = `g511'
scalar `g116' = `g611'
scalar `g117' = `g711'
scalar `g118' = `g811'
scalar `g119' = `g911'
scalar `g1110' = `g1011'
scalar `g1111' = `at'[1,112]
scalar `g1112' = `at'[1,113]
scalar `g1113' = `at'[1,114]

scalar `g121' = `g112'
scalar `g122' = `g212'
scalar `g123' = `g312'
scalar `g124' = `g412'
scalar `g125' = `g512'
scalar `g126' = `g612'
scalar `g127' = `g712'
scalar `g128' = `g812'
scalar `g129' = `g912'
scalar `g1210' = `g1012'
scalar `g1211' = `g1112'
scalar `g1212' = `at'[1,115]
scalar `g1213' = `at'[1,116]

scalar `g131' = `g113'
scalar `g132' = `g213'
scalar `g133' = `g313'
scalar `g134' = `g413'
scalar `g135' = `g513'
scalar `g136' = `g613'
scalar `g137' = `g713'
scalar `g138' = `g813'
scalar `g139' = `g913'
scalar `g1310' = `g1013'
scalar `g1311' = `g1113'
scalar `g1312' = `g1213'
scalar `g1313' = `at'[1,117]

* first 13 lambdas
*
tempname l1 l2 l3 l4 l5 l6 l7 l8 l9 l10 l11 l12 l13
scalar `l1' = `at'[1,118]
scalar `l2' = `at'[1,119]
scalar `l3' = `at'[1,120]
scalar `l4' = `at'[1,121]
scalar `l5' = `at'[1,122]
scalar `l6' = `at'[1,123]
scalar `l7' = `at'[1,124]
scalar `l8' = `at'[1,125]
scalar `l9' = `at'[1,126]
scalar `l10' = `at'[1,127]
scalar `l11' = `at'[1,128]
scalar `l12' = `at'[1,129]
scalar `l13' = `at'[1,130]


**add household demographics variables
*
tempname r11 r12 r13 r14 r15 r16 r17 r18 r19 r110
tempname r21 r22 r23 r24 r25 r26 r27 r28 r29 r210
tempname r31 r32 r33 r34 r35 r36 r37 r38 r39 r310
tempname r41 r42 r43 r44 r45 r46 r47 r48 r49 r410
tempname r51 r52 r53 r54 r55 r56 r57 r58 r59 r510
tempname r61 r62 r63 r64 r65 r66 r67 r68 r69 r610
tempname r71 r72 r73 r74 r75 r76 r77 r78 r79 r710
tempname r81 r82 r83 r84 r85 r86 r87 r88 r89 r810
tempname r91 r92 r93 r94 r95 r96 r97 r98 r99 r910
tempname r101 r102 r103 r104 r105 r106 r107 r108 r109 r1010
tempname r111 r112 r113 r114 r115 r116 r117 r118 r119 r1110
tempname r121 r122 r123 r124 r125 r126 r127 r128 r129 r1210
tempname r131 r132 r133 r134 r135 r136 r137 r138 r139 r1310


scalar `r11' = `at'[1,131]
scalar `r12' = `at'[1,132]
scalar `r13' = `at'[1,133]
scalar `r14' = `at'[1,134]
scalar `r15' = `at'[1,135]
scalar `r16' = `at'[1,136]
scalar `r17' = `at'[1,137]
scalar `r18' = `at'[1,138]
scalar `r19' = `at'[1,139]
scalar `r110' = `at'[1,140]


scalar `r21' = `at'[1,141]
scalar `r22' = `at'[1,142]
scalar `r23' = `at'[1,143]
scalar `r24' = `at'[1,144]
scalar `r25' = `at'[1,145]
scalar `r26' = `at'[1,146]
scalar `r27' = `at'[1,147]
scalar `r28' = `at'[1,148]
scalar `r29' = `at'[1,149]
scalar `r210' = `at'[1,150]

scalar `r31' = `at'[1,151]
scalar `r32' = `at'[1,152]
scalar `r33' = `at'[1,153]
scalar `r34' = `at'[1,154]
scalar `r35' = `at'[1,155]
scalar `r36' = `at'[1,156]
scalar `r37' = `at'[1,157]
scalar `r38' = `at'[1,158]
scalar `r39' = `at'[1,159]
scalar `r310' = `at'[1,160]


scalar `r41' = `at'[1,161]
scalar `r42' = `at'[1,162]
scalar `r43' = `at'[1,163]
scalar `r44' = `at'[1,164]
scalar `r45' = `at'[1,165]
scalar `r46' = `at'[1,166]
scalar `r47' = `at'[1,167]
scalar `r48' = `at'[1,168]
scalar `r49' = `at'[1,169]
scalar `r410' = `at'[1,170]

scalar `r51' = `at'[1,171]
scalar `r52' = `at'[1,172]
scalar `r53' = `at'[1,173]
scalar `r54' = `at'[1,174]
scalar `r55' = `at'[1,175]
scalar `r56' = `at'[1,176]
scalar `r57' = `at'[1,177]
scalar `r58' = `at'[1,178]
scalar `r59' = `at'[1,179]
scalar `r510' = `at'[1,180]

scalar `r61' = `at'[1,181]
scalar `r62' = `at'[1,182]
scalar `r63' = `at'[1,183]
scalar `r64' = `at'[1,184]
scalar `r65' = `at'[1,185]
scalar `r66' = `at'[1,186]
scalar `r67' = `at'[1,187]
scalar `r68' = `at'[1,188]
scalar `r69' = `at'[1,189]
scalar `r610' = `at'[1,190]

scalar `r71' = `at'[1,191]
scalar `r72' = `at'[1,192]
scalar `r73' = `at'[1,193]
scalar `r74' = `at'[1,194]
scalar `r75' = `at'[1,195]
scalar `r76' = `at'[1,196]
scalar `r77' = `at'[1,197]
scalar `r78' = `at'[1,198]
scalar `r79' = `at'[1,199]
scalar `r710' = `at'[1,200]

scalar `r81' = `at'[1,201]
scalar `r82' = `at'[1,202]
scalar `r83' = `at'[1,203]
scalar `r84' = `at'[1,204]
scalar `r85' = `at'[1,205]
scalar `r86' = `at'[1,206]
scalar `r87' = `at'[1,207]
scalar `r88' = `at'[1,208]
scalar `r89' = `at'[1,209]
scalar `r810' = `at'[1,210]

scalar `r91' = `at'[1,211]
scalar `r92' = `at'[1,212]
scalar `r93' = `at'[1,213]
scalar `r94' = `at'[1,214]
scalar `r95' = `at'[1,215]
scalar `r96' = `at'[1,216]
scalar `r97' = `at'[1,217]
scalar `r98' = `at'[1,218]
scalar `r99' = `at'[1,219]
scalar `r910' = `at'[1,220]

scalar `r101' = `at'[1,221]
scalar `r102' = `at'[1,222]
scalar `r103' = `at'[1,223]
scalar `r104' = `at'[1,224]
scalar `r105' = `at'[1,225]
scalar `r106' = `at'[1,226]
scalar `r107' = `at'[1,227]
scalar `r108' = `at'[1,228]
scalar `r109' = `at'[1,229]
scalar `r1010' = `at'[1,230]


*r11, r12, r13 estimated with loops:
loc start=231
forv i=1(1)10 {
scalar `r11`i''=`at'[1,`start']
loc start=`start'+1
}
*
*
loc start=241
forv i=1(1)10 {
scalar `r12`i''=`at'[1,`start']
loc start=`start'+1
}
*
*
loc start=251
forv i=1(1)10 {
scalar `r13`i''=`at'[1,`start']
loc start=`start'+1
}
*
*
*
**pdf
*
tempname d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13
scalar `d1' = `at'[1,261]
scalar `d2' = `at'[1,262]
scalar `d3' = `at'[1,263]
scalar `d4' = `at'[1,264]
scalar `d5' = `at'[1,265]
scalar `d6' = `at'[1,266]
scalar `d7' = `at'[1,267]
scalar `d8' = `at'[1,268]
scalar `d9' = `at'[1,269]
scalar `d10' = `at'[1,270]
scalar `d11' = `at'[1,271]
scalar `d12' = `at'[1,272]
scalar `d13' = `at'[1,273]



// Okay, now that we have all the parameters, we can
// calculate the expenditure shares.
quietly {
// First get the price index
// I set a_0 = 5
tempvar lnpindex
gen double `lnpindex' = 5 + `a1'*`lnp1' + `a2'*`lnp2'+ `a3'*`lnp3' + `a4'*`lnp4'+ `a5'*`lnp5'+ `a6'*`lnp6'+ `a7'*`lnp7'+ `a8'*`lnp8'+ `a9'*`lnp9'+ `a10'*`lnp10'+ `a11'*`lnp11'+ `a12'*`lnp12'+ `a13'*`lnp13'
forvalues i = 1/13 {
forvalues j = 1/13 {
replace `lnpindex' = `lnpindex' + 0.5*`g`i'`j''*`lnp`i''*`lnp`j''
}
}
// The b(p) term in the QUAIDS model:
tempvar bofp
gen double `bofp' = 0
forvalues i = 1/13 {
replace `bofp' = `bofp' + `lnp`i''*`b`i''
}
replace `bofp' = exp(`bofp')
// Finally, the expenditure shares for 17 of the 18
// nutrients (the equation 18 is dropped to avoid singularity)
/* Add extra space between + and `d1' */
replace `w1' = (`a1' + `g11'*`lnp1' + `g12'*`lnp2' +`g13'*`lnp3' + `g14'*`lnp4' + `g15'*`lnp5'+ `g16'*`lnp6' + `g17'*`lnp7' + `g18'*`lnp8' + `g19'*`lnp9' + `g110'*`lnp10' + `g111'*`lnp11'+ `g112'*`lnp12' + `g113'*`lnp13'+ `b1'*(`lnexp' - `lnpindex') + `l1'/`bofp'*(`lnexp' - `lnpindex')^2 +`r11'*`x1' +`r12'*`x2' + `r13'*`x3'+ `r14'*`x4' +`r15'*`x5' +`r16'*`x6' +`r17'*`x7' + `r18'*`x8' + `r19'*`x9' + `r110'*`x10') * `cdf1' + `d1'*`pdf1'


replace `w2' = (`a2' + `g21'*`lnp1' + `g22'*`lnp2' +`g23'*`lnp3' + `g24'*`lnp4' + `g25'*`lnp5' + `g26'*`lnp6' + `g27'*`lnp7' + `g28'*`lnp8' + `g29'*`lnp9' + `g210'*`lnp10' + `g211'*`lnp11' + `g212'*`lnp12' + `g213'*`lnp13' + `b2'*(`lnexp' - `lnpindex') + `l2'/`bofp'*(`lnexp' - `lnpindex')^2 +`r21'*`x1' +`r22'*`x2' + `r23'*`x3'+ `r24'*`x4' +`r25'*`x5' +`r26'*`x6' +`r27'*`x7' + `r28'*`x8' + `r29'*`x9' + `r210'*`x10')*`cdf2' +`d2'*`pdf2'


replace `w3' = (`a3' + `g31'*`lnp1' + `g32'*`lnp2' +`g33'*`lnp3' + `g34'*`lnp4' + `g35'*`lnp5' + `g36'*`lnp6' + `g37'*`lnp7' + `g38'*`lnp8' + `g39'*`lnp9' + `g310'*`lnp10' + `g311'*`lnp11' + `g312'*`lnp12' + `g313'*`lnp13' + `b3'*(`lnexp' - `lnpindex') + `l3'/`bofp'*(`lnexp' - `lnpindex')^2 +`r31'*`x1' +`r32'*`x2' + `r33'*`x3'+ `r34'*`x4' +`r35'*`x5' +`r36'*`x6' +`r37'*`x7' + `r38'*`x8' + `r39'*`x9' + `r310'*`x10')*`cdf3' +`d3'*`pdf3'



replace `w4' = (`a4' + `g41'*`lnp1' + `g42'*`lnp2' +`g43'*`lnp3' + `g44'*`lnp4' + `g45'*`lnp5' + `g46'*`lnp6' + `g47'*`lnp7' + `g48'*`lnp8' + `g49'*`lnp9' + `g410'*`lnp10' + `g411'*`lnp11' + `g412'*`lnp12' + `g413'*`lnp13' + `b4'*(`lnexp' - `lnpindex') + `l4'/`bofp'*(`lnexp' - `lnpindex')^2 +`r41'*`x1' +`r42'*`x2' + `r43'*`x3'+ `r44'*`x4' +`r45'*`x5' +`r46'*`x6' +`r47'*`x7' + `r48'*`x8' + `r49'*`x9' + `r410'*`x10')*`cdf4' +`d4'*`pdf4'



replace `w5' = (`a5' + `g51'*`lnp1' + `g52'*`lnp2' +`g53'*`lnp3' + `g54'*`lnp4' + `g55'*`lnp5' + `g56'*`lnp6' + `g57'*`lnp7' + `g58'*`lnp8' + `g59'*`lnp9' + `g510'*`lnp10' + `g511'*`lnp11' + `g512'*`lnp12' + `g513'*`lnp13' + `b5'*(`lnexp' - `lnpindex') + `l5'/`bofp'*(`lnexp' - `lnpindex')^2 +`r51'*`x1' +`r52'*`x2' + `r53'*`x3'+ `r54'*`x4' +`r55'*`x5' +`r56'*`x6' +`r57'*`x7' + `r58'*`x8' + `r59'*`x9' + `r510'*`x10')*`cdf5' +`d5'*`pdf5'



replace `w6' = (`a6' + `g61'*`lnp1' + `g62'*`lnp2' +`g63'*`lnp3' + `g64'*`lnp4' + `g65'*`lnp5' + `g66'*`lnp6' + `g67'*`lnp7' + `g68'*`lnp8' + `g69'*`lnp9' + `g610'*`lnp10' + `g611'*`lnp11' + `g612'*`lnp12' + `g613'*`lnp13' + `b6'*(`lnexp' - `lnpindex') + `l6'/`bofp'*(`lnexp' - `lnpindex')^2 +`r61'*`x1' +`r62'*`x2' + `r63'*`x3'+ `r64'*`x4' +`r65'*`x5' +`r66'*`x6' +`r67'*`x7' + `r68'*`x8' + `r69'*`x9' + `r610'*`x10')*`cdf6' +`d6'*`pdf6'




replace `w7' = (`a7' + `g71'*`lnp1' + `g72'*`lnp2' +`g73'*`lnp3' + `g74'*`lnp4' + `g75'*`lnp5' + `g76'*`lnp6' + `g77'*`lnp7' + `g78'*`lnp8' + `g79'*`lnp9' + `g710'*`lnp10' + `g711'*`lnp11' + `g712'*`lnp12' + `g713'*`lnp13' + `b7'*(`lnexp' - `lnpindex') + `l7'/`bofp'*(`lnexp' - `lnpindex')^2 +`r71'*`x1' +`r72'*`x2' + `r73'*`x3'+ `r74'*`x4' +`r75'*`x5' +`r76'*`x6' +`r77'*`x7' + `r78'*`x8' + `r79'*`x9' + `r710'*`x10')*`cdf7' +`d7'*`pdf7'


replace `w8' = (`a8' + `g81'*`lnp1' + `g82'*`lnp2' +`g83'*`lnp3' + `g84'*`lnp4' + `g85'*`lnp5' + `g86'*`lnp6' + `g87'*`lnp7' + `g88'*`lnp8' + `g89'*`lnp9' + `g810'*`lnp10' + `g811'*`lnp11' + `g812'*`lnp12' + `g813'*`lnp13' + `b8'*(`lnexp' - `lnpindex') + `l8'/`bofp'*(`lnexp' - `lnpindex')^2 +`r81'*`x1' +`r82'*`x2' + `r83'*`x3'+ `r84'*`x4' +`r85'*`x5' +`r86'*`x6' +`r87'*`x7' + `r88'*`x8' + `r89'*`x9' + `r810'*`x10')*`cdf8' +`d8'*`pdf8'



replace `w9' = (`a9' + `g91'*`lnp1' + `g92'*`lnp2' +`g93'*`lnp3' + `g94'*`lnp4' + `g95'*`lnp5' + `g96'*`lnp6' + `g97'*`lnp7' + `g98'*`lnp8' + `g99'*`lnp9' + `g910'*`lnp10' + `g911'*`lnp11' + `g912'*`lnp12' + `g913'*`lnp13' + `b9'*(`lnexp' - `lnpindex') + `l9'/`bofp'*(`lnexp' - `lnpindex')^2 +`r91'*`x1' +`r92'*`x2' + `r93'*`x3'+ `r94'*`x4' +`r95'*`x5' +`r96'*`x6' +`r97'*`x7' + `r98'*`x8' + `r99'*`x9' + `r910'*`x10')*`cdf9' +`d9'*`pdf9'



replace `w10' = (`a10' + `g101'*`lnp1' + `g102'*`lnp2' +`g103'*`lnp3' + `g104'*`lnp4' + `g105'*`lnp5' + `g106'*`lnp6' + `g107'*`lnp7' + `g108'*`lnp8' + `g109'*`lnp9' + `g1010'*`lnp10' + `g1011'*`lnp11' + `g1012'*`lnp12' + `g1013'*`lnp13' + `b10'*(`lnexp'-`lnpindex') + `l10'/`bofp'*(`lnexp' - `lnpindex')^2 +`r101'*`x1' +`r102'*`x2' + `r103'*`x3'+ `r104'*`x4' +`r105'*`x5' +`r106'*`x6' +`r107'*`x7' + `r108'*`x8' + `r109'*`x9' + `r1010'*`x10')*`cdf10' +`d10'*`pdf10'



replace `w11' = (`a11' + `g111'*`lnp1' + `g112'*`lnp2' +`g113'*`lnp3' + `g114'*`lnp4' + `g115'*`lnp5' + `g116'*`lnp6' + `g117'*`lnp7' + `g118'*`lnp8' + `g119'*`lnp9' + `g1110'*`lnp10' + `g1111'*`lnp11' + `g1112'*`lnp12' + `g1113'*`lnp13' + `b11'*(`lnexp' -`lnpindex') + `l11'/`bofp'*(`lnexp' - `lnpindex')^2 +`r111'*`x1' +`r112'*`x2' + `r113'*`x3'+ `r114'*`x4' +`r115'*`x5' +`r116'*`x6' +`r117'*`x7' + `r118'*`x8' + `r119'*`x9' + `r1110'*`x10')*`cdf11' +`d11'*`pdf11'


replace `w12' = (`a12' + `g121'*`lnp1' + `g122'*`lnp2' +`g123'*`lnp3' + `g124'*`lnp4' + `g125'*`lnp5' + `g126'*`lnp6' + `g127'*`lnp7' + `g128'*`lnp8' + `g129'*`lnp9' + `g1210'*`lnp10' + `g1211'*`lnp11' + `g1212'*`lnp12' + `g1213'*`lnp13' + `b12'*(`lnexp' - `lnpindex') + `l12'/`bofp'*(`lnexp' - `lnpindex')^2 +`r121'*`x1' +`r122'*`x2' + `r123'*`x3'+ `r124'*`x4' +`r125'*`x5' +`r126'*`x6' +`r127'*`x7' + `r128'*`x8' + `r129'*`x9'+ `r1210'*`x10')*`cdf12' +`d12'*`pdf12'

replace `w13' = (`a13' + `g131'*`lnp1' + `g132'*`lnp2' +`g133'*`lnp3' + `g134'*`lnp4' + `g135'*`lnp5' + `g136'*`lnp6' + `g137'*`lnp7' + `g138'*`lnp8' + `g139'*`lnp9' + `g1310'*`lnp10' + `g1311'*`lnp11' + `g1312'*`lnp12' + `g1313'*`lnp13' + `b13'*(`lnexp' - `lnpindex') + `l13'/`bofp'*(`lnexp' - `lnpindex')^2 +`r131'*`x1' +`r132'*`x2' + `r133'*`x3'+ `r134'*`x4' +`r135'*`x5' +`r136'*`x6' +`r137'*`x7' + `r138'*`x8' + `r139'*`x9'+ `r1310'*`x10')*`cdf13' +`d13'*`pdf13'

}
end

/* cdfs have to be added to command */
glo cdfs ""
forv i=1(1)13 {
glo cdfs "${cdfs} cdf`i'"
}

glo A_NOT =5
*
noi nlsur quaids @ w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 w11 w12 w13 lnp1 lnp2 lnp3 lnp4 lnp5 lnp6 lnp7 lnp8 lnp9 lnp10 lnp11 lnp12 lnp13 lnexp x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 pdf1 pdf2 pdf3 pdf4 pdf5 pdf6 pdf7 pdf8 pdf9 pdf10 pdf11 pdf12 pdf13 ${cdfs}, ifgnls nequations(13) param(a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 a11 a12 a13 b1 b2 b3 b4 b5 b6 b7 b8 b9 b10 b11 b12 b13 g11 g12 g13 g14 g15 g16 g17 g18 g19 g110 g111 g112 g113 g22 g23 g24 g25 g26 g27 g28 g29 g210 g211 g212 g213 g33 g34 g35 g36 g37 g38 g39 g310 g311 g312 g313 g44 g45 g46 g47 g48 g49 g410 g411 g412 g413 g55 g56 g57 g58 g59 g510 g511 g512 g513 g66 g67 g68 g69 g610 g611 g612 g613 g77 g78 g79 g710 g711 g712 g713 g88 g89 g810 g811 g812 g813 g99 g910 g911 g912 g913 g1010 g1011 g1012 g1013 g1111 g1112 g1113 g1212 g1213 g1313 l1 l2 l3 l4 l5 l6 l7 l8 l9 l10 l11 l12 l13 r11 r12 r13 r14 r15 r16 r17 r18 r19 r110 r21 r22 r23 r24 r25 r26 r27 r28 r29 r210 r31 r32 r33 r34 r35 r36 r37 r38 r39 r310 r41 r42 r43 r44 r45 r46 r47 r48 r49 r410 r51 r52 r53 r54 r55 r56 r57 r58 r59 r510 r61 r62 r63 r64 r65 r66 r67 r68 r69 r610 r71 r72 r73 r74 r75 r76 r77 r78 r79 r710 r81 r82 r83 r84 r85 r86 r87 r88 r89 r810 r91 r92 r93 r94 r95 r96 r97 r98 r99 r910 r101 r102 r103 r104 r105 r106 r107 r108 r109 r1010 r111 r112 r113 r114 r115 r116 r117 r118 r119 r1110 r121 r122 r123 r124 r125 r126 r127 r128 r129 r1210 r131 r132 r133 r134 r135 r136 r137 r138 r139 r1310 d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12 d13)

[Ahalya Ramanathan]

Dear Statalisters,

I have written the code for calculating elasticities for a QUAIDS model with demographics and censoring containing 13 goods. Despite several attempts it shows errors. It will be great help if someone count correct these codes. The Quaids regression is running perfectly, there is problem only in the elasticity calculation part.
Thanking you in anticipation,
R. Ahalya


****quaids elasticities
est store quaids

*set trace on
*set tracedepth 4

quietly {
foreach x of varlist w* lnp* lnm {
sum `x'
scalar `x'mean=r(mean)
}
* Price indexes
glo asum "_b[a1]*lnp1mean"
forv i=2(1)13 {
glo asum "${asum} + _b[a`i']*lnp`i'mean"
}
glo gsum ""
forv i=1(1)13{
forv j=1(1)13 {
glo gsum "${gsum} + 0.5*_b[g`i'`j']*lnp`i'mean*lnp`j'mean"
}
}
glo ap "5.0 + ${asum} ${gsum}"
glo bp "_b[b1]*lnp1mean"
forv i=2(1)13 {
glo bp "${bp} + _b[b`i']*lnp`i'mean"
}
glo bp "(exp(${bp}))"
* Mus
forv i=1(1)13 {
glo mu`i' "_b[b`i'] + 2*_b[l`i']/${bp}*(lnmmean-(${ap}))"
}
forv j=1(1)13 {
glo gsum2`j' ""
forv k=1(1)13 {
glo gsum2`j' "${gsum2`j'} + _b[g`j'`k']*lnp`k'mean"
}
}
}
nlcom (a1:_b[/a1])(a2:_b[/a2])(a3:_b[/a3])(a4:_b[/a4])(a5:_b[/a5])(a6:_b[/a6]) (a7:_b[/a7])(a8:_b[/a8])(a9:_b[/a9])(a10:_b[/a10])(a11:_b[/a11])(a12:_b[/a12])(a13: 1- _b[/a1] -_b[/a2]- _b[/a3] -_b[/a4]- _b[/a5]- _b[/a6]- _b[/a7] - _b[/a8]- _b[/a9] - _b[/a10] - _b[a11]- _b[/a12])///
(b1:_b[/b1])(b2:_b[/b2])(b3:_b[/b3])(b4:_b[/b4])(b5:_b[/b5])(b6:-_b[/b6])(b7:_b[/b7])(b8:_b[/b8])(b9:_b[/b9])(b10:_b[/b10])(b11:_b[/b11])(b12:_b12[/b12])(b13: -_b[/b1]-_b[/b2]-_b[/b3]-_b[/b4]-_b[/b5]-_b[/b6]-_b[/b7]- _b[/b8]-_b[/b9]- _b[/b10]- _b[/b11] - _b[/b12]) ///
(g11:_b[/g11])(g12:_b[/g12])(g13:_b[/g13])(g14:_b[/g14])(g15:_b[/g15])(g16:_b[/g16])(b17:_b[/g17])(b18:_b[/g18])(b19:_b[/g19])(b110:_b[/g110])(b111:_b[/g111])(b112:_b[/g112]) ///
(g21:_b[/g12])(g22:_b[/g22])(g23:_b[/g23])(g24:_b[/g24])(g25:_b[/g25])(g26:_b[/g26])(g27:_b[/g27])(g28:_b[/g28])(g29:_b[/g29])(g210:_b[/g210])(g211:_b[/g211])(g212;_b[/g212]) ///
(g31:_b[/g13])(g32:_b[/g23])(g33:_b[/g33])(g34:_b[/g34])(g35:_b[/g35])(g36:_b[/g36])(g37:_b[/g37])(g38:_b[/g38])(g39:_b[/g39])(g310:_b[/g310])(g311:_b[/g311])(g312;_b[/g312]) ///
(g41:_b[/g14])(g42:_b[/g24])(g43:_b[/g34])(g44:_b[/g44])(g45:_b[/g45])(g46:_b[/g46])(g47:_b[/g47])(g48:_b[/g48])(g49:_b[/g49])(g410:_b[/g410])(g411:_b[/g411])(g412;_b[/g412]) ///
(g51:_b[/g15])(g52:_b[/g25])(g53:_b[/g35])(g54:_b[/g45])(g55:_b[/g55])(g56:_b[/g56])(g57:_b[/g57])(g58:_b[/g58])(g59:_b[/g59])(g510:_b[/g510])(g511:_b[/g511])(g512;_b[/g512])///
(g61:_b[/g16])(g62:_b[/g26])(g63:_b[/g36])(g64:_b[/g46])(g65:_b[/g56])(g66:_b[/g66])(b67:_b[/g67])(b68:_b[/g68])(b69:_b[/g69])(b610:_b[/g610])(b611:_b[/g611])(b612:_b[/g612]) ///
(g71:_b[/g17])(g72:_b[/g27])(g73:_b[/g37])(g74:_b[/g47])(g75:_b[/g57])(g76:_b[/g67])(g77:_b[/g77])(g78:_b[/g78])(g79:_b[/g79])(g710:_b[/g710])(g711:_b[/g711])(g712;_b[/g712]) ///
(g81:_b[/g18])(g82:_b[/g28])(g83:_b[/g38])(g84:_b[/g48])(g85:_b[/g58])(g86:_b[/g68])(g87:_b[/g78])(g88:_b[/g88])(g89:_b[/g89])(g810:_b[/g810])(g811:_b[/g811])(g812;_b[/g812]) ///
(g91:_b[/g19])(g92:_b[/g29])(g93:_b[/g39])(g94:_b[/g49])(g95:_b[/g59])(g96:_b[/g69])(g97:_b[/g79])(g98:_b[/g89])(g99:_b[/g99])(g910:_b[/g910])(g911:_b[/g911])(g912;_b[/g912]) ///
(g101:_b[/g110])(g102:_b[/g210])(g103:_b[/g310])(g104:_b[/g410])(g105:_b[/g510])(g106:_b[/g610])(g107:_b[/g710])(g108:_b[/g810])(g109:_b[/g910])(g1010:_b[/g1010])(g1011:_b[/g1011])(g1012;_b[/g1012])///
(g111:_b[/g111])(g112:_b[/g211])(g113:_b[/g311])(g114:_b[/g411])(g115:_b[/g511])(g116:_b[/g611])(g117:_b[/g711])(g118:_b[/g811])(g119:_b[/g911])(g1110:_b[/g1011])(g1111:_b[/g1111])(g1112;_b[/g1112]) ///
(g121:_b[/g112])(g122:_b[/g212])(g123:_b[/g312])(g124:_b[/g412])(g125:_b[/g512])(g126:_b[/g612])(g127:_b[/g712])(g128:_b[/g812])(g129:_b[/g912])(g1210:_b[/g1012])(g1211:_b[/g1112])(g1212;_b[/g1212])///

(g113:-_b[/g11]-_b[/g12]-_b[/g13]-_b[/g14]-_b[/g15]- _b[/g16]-_b[/g17]- _b[/g18]-_b[/g19]- _b[/g110]- _b[/g111]- _b[/g112]) ///
(g213:-_b[/g12]-_b[/g22]-_b[/g23]-_b[/g24]-_b[/g25] - _b[/g26]- _b[/g27]-_b[/g28]- _b[/g29]-_b[/g210]- _b[/g211] - _b[/g212]) ///
(g313:-_b[/g13]-_b[/g23]-_b[/g33]-_b[/g34]-_b[/g35]- _b[/g36]- _b[/g37]-_b[/g38] - _b[/g39] -_b[/g310]- _b[/g311]-_b[/g312]) ///
(g413:-_b[/g14]-_b[/g24]-_b[/g34]-_b[/g44]-_b[/g45]-_b[/g46]- _b[/g47]-_b[/g48] - _b[/g49] -_b[/g410]- _b[/g411]-_b[/g412]) ///
(g513:-_b[/g15]-_b[/g25]-_b[/g35]-_b[/g45]-_b[/g55]- _b[/g56]- _b[/g57]-_b[/g58] - _b[/g59] -_b[/g510]- _b[/g511]-_b[/g512]) ///
(g613:-_b[/g16]-_b[/g26]-_b[/g36]-_b[/g46]-_b[/g56]- _b[/g66]- _b[/g67]-_b[/g68] - _b[/g69] -_b[/g610]- _b[/g611]-_b[/g612]) ///
(g713:-_b[/g17]-_b[/g27]-_b[/g37]-_b[/g47]-_b[/g57]- _b[/g67]- _b[/g77]- _b[/g78] - _b[/g79] - _b[/g710] - _b[/g711] - _b[/g712]) ///
(g813:-_b[/g18]-_b[/g28]-_b[/g38]-_b[/g48]-_b[/g58]- _b[/g68]-_b[/g78]-_b[/g88] - _b[/g89] - _b[/g810] - _b[/g811] - _b[/g812]) ///
(g913:-_b[/g19]-_b[/g29]-_b[/g39]-_b[/g49]-_b[/g59] - _b[/g69]-_b[/g79]-_b[/g89] - _b[/g99] - _b[/g910] - _b[/g911] - _b[/g912]) ///
(g1013:-_b[/g110]-_b[/g210]-_b[/g310]-_b[/g410]-_b[/g510])- _b[/g610]- _b[/g710] - _b[/g810] - _b[/g910] - _b[/g1010]- _b[/g1011]- _b[/g1012] )///
(g1113:-_b[/g111]-_b[/g211]-_b[/g311]-_b[/g411]-_b[/g511] - _b[/g611]-_b[/g711]-_b[/g811] - _b[/g911] - _b[/g1011] - _b[/g1111] - _b[/g1112]) ///
(g1213:-_b[/g112]-_b[/g212]-_b[/g312]-_b[/g412]-_b[/g512])- _b[/g612]- _b[/g712] - _b[/g812] - _b[/g912] - _b[/g1012]- _b[/g1112]- _b[/g1212]) ///
(g131:-_b[/g11]-_b[/g12]-_b[/g13]-_b[/g14]-_b[/g15]- _b[/g16]-_b[/g17]- _b[/g18]-_b[/g19]- _b[/g110]- _b[/g111]- _b[/g112]) ///
(g132:-_b[/g12]-_b[/g22]-_b[/g23]-_b[/g24]-_b[/g25] - _b[/g26]- _b[/g27]-_b[/g28]- _b[/g29]-_b[/g210]- _b[/g211] - _b[/g212]) ///
(g133:-_b[/g13]-_b[/g23]-_b[/g33]-_b[/g34]-_b[/g35]- _b[/g36]- _b[/g37]-_b[/g38] - _b[/g39] -_b[/g310]- _b[/g311]-_b[/g312]) ///
(g134:-_b[/g14]-_b[/g24]-_b[/g34]-_b[/g44]-_b[/g45]-_b[/g46]- _b[/g47]-_b[/g48] - _b[/g49] -_b[/g410]- _b[/g411]-_b[/g412]) ///
(g135:-_b[/g15]-_b[/g25]-_b[/g35]-_b[/g45]-_b[/g55]- _b[/g56]- _b[/g57]-_b[/g58] - _b[/g59] -_b[/g510]- _b[/g511]-_b[/g512]) ///
(g136:-_b[/g16]-_b[/g26]-_b[/g36]-_b[/g46]-_b[/g56]- _b[/g66]- _b[/g67]-_b[/g68] - _b[/g69] -_b[/g610]- _b[/g611]-_b[/g612]) ///
(g137:-_b[/g17]-_b[/g27]-_b[/g37]-_b[/g47]-_b[/g57]- _b[/g67]- _b[/g77]- _b[/g78] - _b[/g79] - _b[/g710] - _b[/g711] - _b[/g712]) ///
(g138:-_b[/g18]-_b[/g28]-_b[/g38]-_b[/g48]-_b[/g58]- _b[/g68]-_b[/g78]-_b[/g88] - _b[/g89] - _b[/g810] - _b[/g811] - _b[/g812]) ///
(g139:-_b[/g19]-_b[/g29]-_b[/g39]-_b[/g49]-_b[/g59] - _b[/g69]-_b[/g79]-_b[/g89] - _b[/g99] - _b[/g910] - _b[/g911] - _b[/g912]) ///
(g1310:-_b[/g110]-_b[/g210]-_b[/g310]-_b[/g410]-_b[/g510] - _b[/g610]- _b[/g710] - _b[/g810] - _b[/g910] - _b[/g1010]- _b[/g1011]- _b[/g1012]) ///
(g1311:-_b[/g111]-_b[/g211]-_b[/g311]-_b[/g411]-_b[/g511] - _b[/g611]-_b[/g711]-_b[/g811] - _b[/g911] - _b[/g1011] - _b[/g1111] - _b[/g1112]) ///
(g1312:-(-_b[/g11]-_b[/g12]-_b[/g13]-_b[/g14]-_b[/g15]-_b[/g16]-_b[/g17]- _b[/g18]-_b[/g19] - _b[/g110] - _b[/g111]- _b[/g112])- ///
(-_b[/g12]-_b[/g22]-_b[/g23]-_b[/g24]-_b[/g25] - _b[/g26]- _b[/g27]-_b[/g28]- _b[/g29]-_b[/g210]- _b[/g211] - _b[/g212])-///
(-_b[/g13]-_b[/g23]-_b[/g33]-_b[/g34]-_b[/g35]- _b[/g36]- _b[/g37]-_b[/g38] - _b[/g39] -_b[/g310]- _b[/g311]-_b[/g312])-///
(-_b[/g14]-_b[/g24]-_b[/g34]-_b[/g44]-_b[/g45]-_b[/g46]- _b[/g47]-_b[/g48] - _b[/g49] -_b[/g410]- _b[/g411]-_b[/g412])-///
(-_b[/g15]-_b[/g25]-_b[/g35]-_b[/g45]-_b[/g55]- _b[/g56]- _b[/g57]-_b[/g58] - _b[/g59] -_b[/g510]- _b[/g511]-_b[/g512])-///
(-_b[/g16]-_b[/g26]-_b[/g36]-_b[/g46]-_b[/g56]- _b[/g66]- _b[/g67]-_b[/g68] - _b[/g69] -_b[/g610]- _b[/g611]-_b[/g612])-///
(-_b[/g17]-_b[/g27]-_b[/g37]-_b[/g47]-_b[/g57]- _b[/g67]- _b[/g77]- _b[/g78] - _b[/g79] - _b[/g710] - _b[/g711] - _b[/g712])-///
(-_b[/g18]-_b[/g28]-_b[/g38]-_b[/g48]-_b[/g58]- _b[/g68]-_b[/g78]-_b[/g88] - _b[/g89] - _b[/g810] - _b[/g811] - _b[/g812]) -///
(-_b[/g19]-_b[/g29]-_b[/g39]-_b[/g49]-_b[/g59] - _b[/g69]-_b[/g79]-_b[/g89] - _b[/g99] - _b[/g910] - _b[/g911] - _b[/g912])-///
(-_b[/g110]-_b[/g210]-_b[/g310]-_b[/g410]-_b[/g510] - _b[/g610]- _b[/g710] - _b[/g810] - _b[/g910] - _b[/g1010]- _b[/g1011]- _b[/g1012])-///
(-_b[/g111]-_b[/g211]-_b[/g311]-_b[/g411]-_b[/g511] - _b[/g611]-_b[/g711]-_b[/g811] - _b[/g911] - _b[/g1011] - _b[/g1111] - _b[/g1112])///

(r11:_b[/r11])(r12:_b[/r12])(r13:_b[/r13])(r14:_b[/r14])(r15:_b[/r15])(r16:_b[/r16])(r17:_b[/r17])(r18:_b[/r18])(r19:_b[/r19])(r110:_b[/r110])(r111:_b[/r111]) ///
(r112:-_b[/r11]-_b[/r12]-_b[/r13]-_b[/r14]-_b[/r15]-_b[/r16]- _b[/r17]-_b[/r18]-_b[/r19]-_b[/r110-_b[/r111]) ///
(r21:_b[/r21])(r22:_b[/r22])(r23:_b[/r23])(r24:_b[/r24])(r25:_b[/r25])(r26:_b[/r26])(r27:_b[/r27])(r28:_b[/r28])(r29:_b[/r29])(r210:_b[/r210])(r211:_b[/r211]) ///
(r212:-_b[/r21]-_b[/r22]-_b[/r23]-_b[/r24]-_b[/r25]-_b[/r26]-_b[/r27]-_b[/r28]-_b[/r29]-_b[/r210]-_b[/r211]) ///
(r31:_b[/r31])(r32:_b[/r32])(r33:_b[/r33])(r34:_b[/r34])(r35:_b[/r35])(r36:_b[/r36])(r37:_b[/r37])(r38:_b[/r38])(r39:_b[/r39])(r310:_b[/r310])(r311:_b[/r311]) ///
(r312:-_b[/r31]-_b[/r32]-_b[/r33]-_b[/r34]-_b[/r35]-_b[/r36]-_b[/r37]-_b[/r38]-_b[/r39]-_b[/r310]-_b[/r311]) ///
(r41:_b[/r41])(r42:_b[/r42])(r43:_b[/r43])(r44:_b[/r44])(r45:_b[/r45])(r46:_b[/r46])(r47:_b[/r47])(r48:_b[/r48])(r49:_b[/r49])(r410:_b[/r410])(r411:_b[/r411]) ///
(r412:-_b[/r41]-_b[/r42]-_b[/r43]-_b[/r44]-_b[/r45]-_b[/r46]-_b[/r47]-_b[/r48]-_b[/r49]-_b[/r410]-_b[/r411]) ///
(r51:_b[/r51])(r52:_b[/r52])(r53:_b[/r53])(r54:_b[/r54])(r55:_b[/r55])(r56:_b[/r56])(r57:_b[/r57])(r58:_b[/r58])(r59:_b[/r59])(r510:_b[/r510])(r511:_b[/r511]) ///
(r512:-_b[/r51]-_b[/r52]-_b[/r53]-_b[/r54]-_b[/r55]-_b[/r56]-_b[/r57]-_b[/r58]-_b[/r59]-_b[/r510]-_b[/r511]) ///
(r61:_b[/r61])(r62:_b[/r62])(r63:_b[/r63])(r64:_b[/r64])(r65:_b[/r65])(r66:_b[/r66])(r67:_b[/r67])(r68:_b[/r68])(r69:_b[/r69])(r610:_b[/r610])(r611:_b[/r611]) ///
(r612:-_b[/r61]-_b[/r62]-_b[/r63]-_b[/r64]-_b[/r65]-_b[/r66]-_b[/r67]-_b[/r68]-_b[/r69]-_b[/r610]-_b[/r611]) ///
(r71:_b[/r71])(r72:_b[/r72])(r73:_b[/r73])(r74:_b[/r74])(r75:_b[/r75])(r76:_b[/r76])(r77:_b[/r77])(r78:_b[/r78])(r79:_b[/r79])(r710:_b[/r710])(r711:_b[/r711]) ///
(r712:-_b[/r71]-_b[/r72]-_b[/r73]-_b[/r74]-_b[/r75]-_b[/r76]-_b[/r77]-_b[/r78]-_b[/r79]-_b[/r710]-_b[/r711]) ///
(r81:_b[/r81])(r82:_b[/r82])(r83:_b[/r83])(r84:_b[/r84])(r85:_b[/r85])(r86:_b[/r86])(r87:_b[/r87])(r88:_b[/r88])(r89:_b[/r89])(r810:_b[/r810])(r811:_b[/r811]) ///
(r812:-_b[/r81]-_b[/r82]-_b[/r83]-_b[/r84]-_b[/r85]-_b[/r86]-_b[/r87]-_b[/r88]-_b[/r89]-_b[/r810]-_b[/r811]) ///
(r91:_b[/r91])(r92:_b[/r92])(r93:_b[/r93])(r94:_b[/r94])(r95:_b[/r95])(r96:_b[/r96])(r97:_b[/r97])(r98:_b[/r98])(r99:_b[/r99])(r910:_b[/r910])(r911:_b[/r911]) ///
(r912:-_b[/r91]-_b[/r92]-_b[/r93]-_b[/r94]-_b[/r95]-_b[/r96]-_b[/r97]-_b[/r98]-_b[/r99]-_b[/r910]-_b[/r911]) ///
(r101:_b[/r101])(r102:_b[/r102])(r103:_b[/r103])(r104:_b[/r104])(r105:_b[/r105])(r106:_b[/r106])(r107:_b[/r107])(r108:_b[/r108])(r109:_b[/r109])(r1010:_b[/r1010])(r1011:_b[/r1011]) ///
(r1012:-_b[/r101]-_b[/r102]-_b[/r103]-_b[/r104]-_b[/r105]-_b[/r106]-_b[/r107]-_b[/r108]-_b[/r109]-_b[/r1010]-_b[/r1011]) ///
(l1:_b[/l1])(l2:_b[/l2])(l3:_b[/l3])(l4:_b[/l4])(l5:_b[/l5])(16:_b[/16])(17:_b[/17])(18:_b[/18])(19:_b[/19])(110:_b[/110])(111:_b[/111]) ///
(l12:-_b[/l1]-_b[/l2]-_b[/l3]-_b[/l4]-_b[/l5]-_b[/16]-_b[/17]-_b[/18]-_b[/19]-_b[/110]-_b[/111]) ///
(d1:_b[/d1]) (d2:_b[/d2]) (d3:_b[/d3]) (d4:_b[/d4])(d5:_b[/d5])(d6:_b[/d6])(d7:_b[/d7])(d8:_b[/d8])(d9:_b[/d9])(d10:_b[/d10])(d11:_b[/d11])(d12:_b[/d12]), post noheader


est store quaids

forv i=1(1)13 {
forv j=1(1)13 {
glo delta=cond(`i'==`j',1,0)
glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j'] ${gsum2`j'})-_b[l`i']*_b[b`j']/${bp}*(lnmmean - (${ap}))^2"
cap nlcom (elasm`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}) (d`i': _b[d`i']), post noheader
if _rc {
qui nlcom (elasm`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f1e+2)*(${mu`i'`j'})) (d`i': _b[d`i']) , post noheader
qui nlcom (elasm`i': _b[elasm`i']) (mu`i'`j':_b[mu`i'`j'f]/(1e+2)) (d`i': _b[d`i']) , post noheader
}

* Uncompensated price elasticity
if _b[d`i']!=0{
nlcom (elasm`i': _b[elasm`i']*_b[d`i']) (elu`i'`j'_b[mu`i'`j']/w`i'mean - ${delta})*_b[d`i']), post noheader
}
else {
display _b[d`i']
nlcom (elasm`i': _b[elasm`i']) (elu`i'`j'_b[mu`i'`j']/w`i'mean - ${delta})), post noheader
}
* Compensated price elasticity
nlcom (elc`i'`j': _b[elu`i'`j'] + _b[elasm`i']*w`j'mean), noheader
qui est restore quaids
}
}

[Ahalya Ramanathan]

David Achiles you had posted the commands for calculating elasticity with 6 products a few years back. It will be great help if you can provide me the corrected code for elasticity calculation. It repeatedly says g43 not found. Anyone who has the corrected codes, kindly help please.

[Juan Caro]

Hi Everyone,

I have been working on a generalized version of the code for censored QUAIDS, so I figured it is appropriate to post here. I have everything automated for J goods (including elasticities), however I cannot find a way to create a varlist that automatically completes the "args". I leave the evaluator function below, maybe someone can use it. Regards,

cap program drop nlsurquaidsNNc
program nlsurquaidsNNc
version 14.0
syntax varlist
tokenize `varlist'

*** needs to be updated manually
args w1 w2 w3 w4 w5 w6 w7 w8 w9 lnp1 lnp2 lnp3 lnp4 lnp5 lnp6 lnp7 lnp8 lnp9 lnexp ///
x1 x2 x3 elnexp pdf1 pdf2 pdf3 pdf4 pdf5 pdf6 pdf7 pdf8 pdf9 cdf1 cdf2 cdf3 cdf4 cdf5 cdf6 cdf7 cdf8 cdf9

*****************
*variables must enter the program in the following order and with the same names as quaids_.ado:
*args w1...wJ lnp1...lnpJ lnm x1 x2 x3 x4 pdf1...pdfJ cdf1...cdfJ
*****************

***Coefficients

qui describe `varlist'
local J= (r(k)-5)/4
local J1 =`J'-1
local aux=0

forvalues j=1(1)`J1' {
tempname a`j'
scalar `a`j'' = `at'[1,`j']
replace `aux' = `aux' + `a`j''
}
tempname a`J'
scalar `a`J'' = 1-`aux'

replace `aux'=0
forvalues j=1(1)`J1' {
tempname b`j'
scalar `b`j'' = `at'[1,`J1'+`j']
replace `aux' = `aux' + `b`j''
}
tempname b`J'
scalar `b`J'' = -`aux'

forvalues j=1(1)`J1' {
local aux`j'=0
local aux0=2*`J1'
forvalues k=1(1)`J1' {
tempname g`j'`k'
if `k'>=`j' {
scalar `g`j'`k'' = `at'[1,`aux0'+(`k'-`j'+1)]
}
else {
scalar `g`j'`k'' = `g`k'`j''
}
replace `aux`j'' = aux`j' + `g`j'`k''
replace `aux0'=`aux0'+`J1'-`j'
}
}

forvalues j=1(1)`J' {
tempname `g`j'`J''
scalar `g`j'`J'' = -`aux`j''
}

replace `aux'=0
replace `aux0'=`aux0'+1
forvalues j=1(1)`J1' {
tempname l`j'
scalar `l`j'' = `at'[1,`aux0'+`j']
replace `aux' = `aux' + `l`j''
}
tempname l`J'
scalar `l`J'' = -`aux'


//household demographics

replace `aux'=0
replace `aux0'=`aux0'+`J1'
forvalues j=1(1)`J1' {
tempname r`j'1
scalar `r`j'1' = `at'[1,`aux0'+`j']
replace `aux' = `aux' + `r`j'1'
}
tempname r`J'1
scalar `r`J'1' = -`aux'

replace `aux'=0
replace `aux0'=`aux0'+`J1'
forvalues j=1(1)`J1' {
tempname r`j'2
scalar `r`j'2' = `at'[1,`aux0'+`j']
replace `aux' = `aux' + `r`j'2'
}
tempname r`J'2
scalar `r`J'2' = -`aux'

replace `aux'=0
replace `aux0'=`aux0'+`J1'
forvalues j=1(1)`J1' {
tempname r`j'3
scalar `r`j'3' = `at'[1,`aux0'+`j']
replace `aux' = `aux' + `r`j'3'
}
tempname r`J'3
scalar `r`J'3' = -`aux'

replace `aux'=0
replace `aux0'=`aux0'+`J1'
forvalues j=1(1)`J1' {
tempname r`j'4
scalar `r`j'4' = `at'[1,`aux0'+`j']
replace `aux' = `aux' + `r`j'4'
}
tempname r`J'4
scalar `r`J'4' = -`aux'

replace `aux'=0
forvalues j=1(1)`J' {
tempname d`j'
scalar `d`j'' = `at'[1,`aux0'+`j']
}


***indexes

quietly {
forvalues i=1(1)`J' {
tempvar aa`i'
gen double `aa`i'' = `a`i''+`r`i'1'*`x1' +`r`i'2'*`x2' +`r`i'3'*`x3' +`r`i'4'*`x4'
}

tempvar lnpindex
gen double `lnpindex' = 5
forvalues i = 1(1)`J' {
replace `lnpindex' = `lnpindex' + `aa`i''*`lnp`i''
}
forvalues i = 1(1)`J' {
forvalues j = 1(1)`J' {
replace `lnpindex' = `lnpindex' + 0.5*`g`i'`j''*`lnp`i''*`lnp`j''
}
}

// The b(p) term in the QUAIDS model:
tempvar bofp
gen double `bofp' = 0
forvalues i = 1(1)`J' {
replace `bofp' = `bofp' + `lnp`i''*`b`i''
}
replace `bofp' = exp(`bofp')

// Finally, the expenditure shares

forvalues i = 1(1)`J' {
replace `w`i'' = `b`i''*(`lnm' - `lnpindex') + `l`i''/`bofp'*(`lnm' - `lnpindex')^2

forvalues j = 1(1)`J' {
replace `w`i'' = `w`i''+`aa`i'' + `g`i'`j''*`lnp`j''
}
replace `w`i''=`w`i''*`cdf`i'' + `d`i''*`pdf`i''
}

}

end

[Chutarat Noosuwan]

Dear Statalisters,

I am applying this QUAIDS code but still have an error message when I am trying to estimate the elasticities.

I got the error said:

expression too long
r(130);

after the command:

glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j']${gsum2`j'})- _b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2"

based on this code:

est store quaidsNNP2 
forv i=1(1)10 {
forv j=1(1)10 {
glo delta=cond(`i'==`j',1,0)
glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j']${gsum2`j'})- _b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2"
cap nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}) (d`i':_b[d`i']), post noheader
if _rc {
est restore quaidsNNP2
nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f:(100)*(${mu`i'`j'})) (d`i':_b[d`i']), post noheader
nlcom (elasexp`i': _b[elasexp`i']) (mu`i'`j':_b[mu`i'`j'f]/(100)) (d`i':_b[d`i']), post noheader
}
*
*Uncompensated price elasticity
if _b[d`i']!=0{
nlcom (elasexp`i': _b[elasexp`i']*_b[d`i']) (elu`i'`j':(_b[mu`i'`j']/w`i'mean - ${delta})*_b[d`i']), post noheader
}
else {
display _b[d`i']
nlcom (elasexp`i': _b[elasexp`i']) (elu`i'`j':(_b[mu`i'`j']/w`i'mean - ${delta})), post noheader
}
*
*Compensated price elasticity
nlcom (elc`i'`j': _b[elu`i'`j'] + _b[elasexp`i']*w`j'mean), noheader
qui est restore quaidsNNP2
}
}

I have tried to split it such as

glo muc`i'`j' "(mua`i'`j' : _b[g`i'`j'] - ${mu`i'}*(_b[a`j']${gsum2`j'})) (mub`i'`j' :  _b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2)"
glo mu`i'`j' "$mua`i'`j' - $mub`i'`j'"

However, it doesn't work. The Stata can estimated it but the elasticity values are not include the values of

(${mu`i'`j'}))

that mentioned in

glo mu`i'`j' "$mua`i'`j' - $mub`i'`j'"

Any suggestions will be much appreciated.

Thank you in advance.

Chutarat

[Jiao Chen]

Dear Statalisters,

I am applying this QUAIDS code but still have an error message when I am trying to estimate the elasticities.

I got the error said:

expression too long
r(130);

after the command:

glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j']${gsum2`j'})- _b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2"

based on this code:

est store quaidsNNP2
forv i=1(1)10 {
forv j=1(1)10 {
glo delta=cond(`i'==`j',1,0)
glo mu`i'`j' "_b[g`i'`j'] - ${mu`i'}*(_b[a`j']${gsum2`j'})- _b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2"
cap nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}) (d`i':_b[d`i']), post noheader
if _rc {
est restore quaidsNNP2
nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f:(100)*(${mu`i'`j'})) (d`i':_b[d`i']), post noheader
nlcom (elasexp`i': _b[elasexp`i']) (mu`i'`j':_b[mu`i'`j'f]/(100)) (d`i':_b[d`i']), post noheader
}
*
*Uncompensated price elasticity
if _b[d`i']!=0{
nlcom (elasexp`i': _b[elasexp`i']*_b[d`i']) (elu`i'`j':(_b[mu`i'`j']/w`i'mean - ${delta})*_b[d`i']), post noheader
}
else {
display _b[d`i']
nlcom (elasexp`i': _b[elasexp`i']) (elu`i'`j':(_b[mu`i'`j']/w`i'mean - ${delta})), post noheader
}
*
*Compensated price elasticity
nlcom (elc`i'`j': _b[elu`i'`j'] + _b[elasexp`i']*w`j'mean), noheader
qui est restore quaidsNNP2
}
}

I have tried to split it such as

glo muc`i'`j' "(mua`i'`j' : _b[g`i'`j'] - ${mu`i'}*(_b[a`j']${gsum2`j'})) (mub`i'`j' : _b[l`i']*_b[b`j']/${bp}*(lnexpmean - (${ap}))^2)"
glo mu`i'`j' "$mua`i'`j' - $mub`i'`j'"

However, it doesn't work. The Stata can estimated it but the elasticity values are not include the values of

(${mu`i'`j'}))

that mentioned in

glo mu`i'`j' "$mua`i'`j' - $mub`i'`j'"

Any suggestions will be much appreciated.

Thank you in advance.

Chutarat

Hello, Chutarat,
I meet the same problem as you(

expression too long r(130);

). I believe you must have solved this problem, could you share it with me?
or somebody else could help me ?
Thank you in advance

Best,
Jiao