Problem with QUAIDS elasticities using nlsur command

[Mariko Wijekoon]

Hi Stata listers

I am trying to estimate a 13 goods QUAIDS model that corrects for the problem of censoring in the dependent variable (zero expenditure problem) using the Shonkwiler and Yen method (1999): Shonkwiler, J. Scott, and Steven T. Yen. "Two-step estimation of a censored system of equations." American Journal of Agricultural Economics 81.4 (1999): 972-982.

Brian Poi’s (2012) “quaids” command does not accommodate this correction for censoring, and therefore, I am using his older “nlsur” command to estimate the QUAIDS model. I am very grateful to Jorge, who helped me by debugging the problems I had in my commands up to the estimation of expenditure share equations.However, now I am having problems with estimating the elasticities and therefore seeking help from the Stata community.

The problem seem to occur right after this command in the elasticities part of the file: qui nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f: (1e+2)*(${mu`i'`j'})), post noheader

I get the following error message: [b1] not found

These are some suggestions were made by Jorge in the past for the same problem: http://www.stata.com/statalist/archi.../msg00572.html

1. "It might be that you have an extra : somewhere in a -nlcom- call
2. or a missing _b (all the calls to b coefficients should read _b[b1], _b[b2] ... or _b[b`i'] when inside the loop.
3. or that the calls to nlcom are not finding the estimates previously stored, and that would be why _b[b1] is not found.

If your problem is among the first two, you have to check the code again. As for the third potential problem , type -ereturn list- before running the elasticity calculations and make sure e(b) is listed under matrices"

I tried typing ereturn list and e(b) is listed and therefore it cannot be the 3^rd problem. I am not sure how to proceed with the first two suggestions.

Any advice/suggestions will be appreciated.

Thank you very much in advance.

Mariko Wijekoon.


Here are the Stata codes:


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)

est store quaidsmariko

set trace on
set tracedepth 4

* Share means and price 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)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 "6.11 + ${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}*(lnexpmean-(${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"
}
}
}
*
*ereturn list
*

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}*(lnexpmean - (${ap}))^2"
* If expression is too long, split it
cap nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}), post noheader
if _rc {
*****the problem begins here**********
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']) (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 quaidsmariko
}
}

[Jorge Eduardo Perez Perez]

Hi Mariko.

My guess is that one of the "mu`i'" variables is not well defined. I can't find an error of the code, but maybe another set of eyes will.

Could you paste the output of ereturn list when executred right before the elasticity calculations as well as the output of

forv i=1(1)13 {
di "${mu`i'}"
}

[Nick Cox]

1e+2 is a strange way to write 100! Perhaps the + sign is being misparsed. Unlikely, but may be worth a try.

[Mariko Wijekoon]

Thanks Jorge

This is what ereturn list gives me just before the elasticities calculation.

The output for for code you sent me is on a word document attached here.

Nick - I will try running the codes using 100 instead of (1e+2) and see what happens.


ereturn list

scalars:
e(mss_5) = 1.867691324673476
e(rss_5) = 1.236914032588635
e(rmse_4) = .0633976062532114
e(r2_4) = .7030670676658199
e(mss_4) = 4.339592516919401
e(rss_4) = 1.832780954258579
e(rmse_3) = .0200351415223624
e(r2_3) = .7246737420932414
e(mss_3) = .481775338159782
e(rss_3) = .1830415444944174
e(rmse_2) = .0226295978859409
e(r2_2) = .7982219443802627
e(mss_2) = .9237793432741365
e(rss_2) = .2335170074185974
e(rmse_1) = .0802441606605177
e(r2_1) = .6520306927189259
e(mss_1) = 5.501977640949134
e(rss_1) = 2.936241145970609
e(ic) = 129
e(N) = 456
e(k_eq) = 13
e(rank) = 244
e(converge) = 1
e(k_eq_model) = 0
e(k) = 273
e(n_eq) = 13
e(ll) = 10718.65304795093
e(rmse_13) = .2256336569177784
e(r2_13) = .245036525696731
e(mss_13) = 7.534899993369606
e(rss_13) = 23.21520949314491
e(rmse_12) = .0142522421861945
e(r2_12) = .2146557173818122
e(mss_12) = .0253170794218859
e(rss_12) = .0926256417442776
e(rmse_11) = .1092733286633729
e(r2_11) = .6833665787108133
e(mss_11) = 11.75141515784558
e(rss_11) = 5.444941122871126
e(rmse_10) = .0121006864408387
e(r2_10) = .5372624892480251
e(mss_10) = .0775240889939594
e(rss_10) = .0667705352268108
e(rmse_9) = .0826751215129467
e(r2_9) = .51191240015064
e(mss_9) = 3.2689810411254
e(rss_9) = 3.116840127034308
e(rmse_8) = .039066706074255
e(r2_8) = .3575902643370135
e(mss_8) = .3873932105732222
e(rss_8) = .6959506307124587
e(rmse_7) = .0346154224780695
e(r2_7) = .5534019259445773
e(mss_7) = .6770612147104713
e(rss_7) = .5463917278408679
e(rmse_6) = .0551499957945507
e(r2_6) = .4644055629046585
e(mss_6) = 1.202588830064792
e(rss_6) = 1.386934048479368
e(rmse_5) = .0520819620715339
e(r2_5) = .6015873548322919

macros:
e(_estimates_name) : "quaidsmariko"
e(cmdline) : "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.."
e(cmd) : "nlsur"
e(predict) : "nlsur_p"
e(estat_cmd) : "nlsur_estat"
e(covariates) : "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 .."
e(rhs) : "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 .."
e(funcprog) : "nlsurquaids"
e(vce) : "gnr"
e(method) : "ifgnls"
e(constants) : ". . . . . . . . . . . . ."
e(type) : "3"
e(params) : "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 g.."
e(depvar) : "w1 w2 w3 w4 w5 w6 w7 w8 w9 w10 w11 w12 w13"
e(depvar_13) : "w13"
e(depvar_12) : "w12"
e(depvar_11) : "w11"
e(depvar_10) : "w10"
e(depvar_9) : "w9"
e(depvar_8) : "w8"
e(depvar_7) : "w7"
e(depvar_6) : "w6"
e(depvar_5) : "w5"
e(depvar_4) : "w4"
e(depvar_3) : "w3"
e(depvar_2) : "w2"
e(depvar_1) : "w1"
e(properties) : "b V"

matrices:
e(b) : 1 x 273
e(V) : 273 x 273
e(init) : 1 x 273
e(Sigma) : 13 x 13

functions:
e(sample)

[Jorge Eduardo Perez Perez]

Try adding

est restore quaidsmariko

Right before the last piece of code that calculates the elasticities

...

*ereturn list
*
est restore quaidsmariko

forv i=1(1)13 {
forv j=1(1)13 {

...

[Mariko Wijekoon]

Jorge - I added the code: est restore quaidsmariko

However I get the same error message: [b1] not found.

Nick - I tried substituting 100 instead of (1e+2) however it did not change anything.

[Jorge Eduardo Perez Perez]

Try adding

est restore quaidsmariko

right after

*****the problem begins here**********

[Mariko Wijekoon]

Jorge - I did that - but I still get the same message.

Just to clarify that I did exactly what you suggested- this is how I placed the commands that I ran after the budget share estimations:


est store quaidsmariko

set trace on
set tracedepth 4

* Share means and price 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)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 "6.11 + ${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}*(lnexpmean-(${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"
}
}
}
*
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}*(lnexpmean - (${ap}))^2"
* If expression is too long, split it
cap nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j': ${mu`i'`j'}), post noheader
if _rc {
*****the problem begins here**********
*
est restore quaidsmariko
*
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']) (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 quaidsmariko
*
}
}

[Mariko Wijekoon]

I experimented with a simple 5 good case and before estimating the elasticities I posted the coefficients using the nlcom command as given below and the missing b1 problem was solved and I was able to get the elasticities estimates.

However for the 13 goods case - I got an error message that says: expression too long

The error appears right after the command:

qui nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f: (1e+2)*(${mu`i'`j'})), post noheader

I read that the solution for this is to split the command -but I am not sure how to do it. Any advice/suggestions will be appreciated.

Thank you in advance

Mariko Wijekoon

nlcom (a1:_b[/a1]) (a2:_b[/a2]) (a3:_b[/a3]) ///
(a4:_b[/a4]) (a5:_b[/a5]) ///
(b1:_b[/b1]) (b2:_b[/b2]) (b3:_b[/b3]) ///
(b4:_b[/b4]) (b5:_b[/b5]) ///
(g11:_b[/g11]) (g12:_b[/g12]) (g13:_b[/g13]) ///
(g14:_b[/g14]) (g15:_b[/g15]) ///
(g21:_b[/g12]) (g22:_b[/g22]) (g23:_b[/g23]) ///
(g24:_b[/g24]) (g25:_b[/g25]) ///
(g31:_b[/g13]) (g32:_b[/g23]) (g33:_b[/g33]) ///
(g34:_b[/g34]) (g35:_b[/g35]) ///
(g41:_b[/g14]) (g42:_b[/g24]) (g43:_b[/g34]) ///
(g44:-_b[/g44]) (g45:-_b[/g45]) ///
(g51:_b[/g15]) (g52:_b[/g25]) (g53:_b[/g35]) ///
(g54:-_b[/g45]) (g55:-_b[/g55]) ///
(l1:_b[/l1]) (l2:_b[/l2]) (l3:_b[/l3]) ///
(l4:_b[/l4]) (l5:_b[/l5]), noheader post

[Vanclei Zanin]

Dear Mariko, I am trying to use the same approach - [Shonkwiler and Yen method (1999): Shonkwiler, J. Scott, and Steven T. Yen. "Two-step estimation of a censored system of equations." American Journal of Agricultural Economics 81.4 (1999): 972-982.] - to estimate a demand system with 16 products. I am in the beginning, and a I have some silly questions. In your example here you don´t exclud one good to avoid singularity, Can I do the same?
The second one is if you tried to estimate elasticities for demographics components, if yes, could you share the code, to help me?
thanks a lot

[Nigussie Tefera]

Dear Mariko,
Does the code finally work for you? I was tried to replicate it but it is still pop up with a message of "b[1] not found". Please let me know.

[Jorge Eduardo Perez Perez]

Dear Nigussie.

I can't tell from here if you're using the same code as Mariko. Please post the code you were using, exactly as you ran it. Also, from this post, keep in mind that you need all -nlcom- calls should be in the same line (i.e. if you copy pasted this code and the lines were split, you need to rejoin them) and that you need to restore the estimates before calculating elasticities, as suggested in this thread.

[Nigussie Tefera]

Thanks Jorge,

Please find, here is the do file. I also included the results as it is reported (just below do file).

FYI: The is the stata code I modified to adjust for bias due zero consumption expenditure and socioeconomic variables. I already posted my original ado file on stata website in 2013. Thus, this a test file whether the results are consistent with Poi's QUAIDS estimation results by considering all CDF's=1 and PDF's=0 in the system. The results are consistent with Poi's quaids estimation but I want to estimates all elasticity using yours simplified versions of the estimation techniques. Besides, it is worth to note that all elasticity parameters should have to be adjusted for CDF and PDF parameters once the system is working properly. I will work on the modification of the last section but really want to know how to solve the problem that reports "b[1] is not found".

Best

Nigussie

use "G:\ERHS\do\Nigussie\Demand\Examples\food.dta" , clear

* TEST modified model gives same results

set trace off


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

*** create additional variables
generate x1 = int(runiform()*4) // generate number of kinds
generate x2 = (runiform() > 0.7) // generate rural vs urban binary variables

* Nigussie's approaches
gen cdf1=1
gen cdf2=1
gen cdf3=1

gen pdf1=0
gen pdf2=0
gen pdf3=0

nlsur quaidsNNP @ w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnexp ///
x1 x2 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 quaidsNNP2

*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)3 {
glo asum "${asum} + _b[a`i']*lnp`i'mean"
}
glo gsum ""
forv i=1(1)3 {
forv j=1(1)3 {
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)3 {
glo bp "${bp} + _b[b`i']*lnp`i'mean"
}
glo bp "(exp(${bp}))"
* Mus
forv i=1(1)3 {
glo mu`i' "_b[b`i'] + 2*_b[l`i']/${bp}*(lnexp-(${ap}))"
}
forv j=1(1)3 {
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: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]), post noheader

forv i=1(1)3 {
forv j=1(1)3 {
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 {
est restore quaidsNNP2
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']) (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
}
}

The result is given below
. use "G:\ERHS\do\Nigussie\Demand\Examples\food.dta" , clear

.
. set trace off

.
. * TEST modified model gives same results
. ****************************************
.
. *** create additional variables
. generate x1 = int(runiform()*4) // generate number of kinds

. generate x2 = (runiform() > 0.7) // generate rural vs urban binary variables

.
. * Nigussie's approaches
. gen cdf1=1

. gen cdf2=1

. gen cdf3=1

.
. gen pdf1=0

. gen pdf2=0

. gen pdf3=0

.
. nlsur quaidsNNP @ w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnexp ///
> x1 x2 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
(obs = 4,048)
Calculating NLS estimates...
Calculating FGNLS estimates...
FGNLS iteration 2...
FGNLS iteration 3...

FGNLS regression
-----------------------------------------------------------------------
Equation | Obs Parms RMSE R-sq Constant
----------------+------------------------------------------------------
1 w1 | 4,048 . .1333049 0.9017* (none)
2 w2 | 4,048 . .1024036 0.8481* (none)
3 w3 | 4,048 . .0537615 0.7907* (none)
-----------------------------------------------------------------------
* Uncentered R-sq

------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
/a1 | .3145464 .0092701 33.93 0.000 .2963772 .3327155
/a2 | .2708841 .0071607 37.83 0.000 .2568494 .2849187
/a3 | .1063753 .003666 29.02 0.000 .0991901 .1135604
/g11 | .0098714 .0107752 0.92 0.360 -.0112477 .0309905
/g21 | -.0262766 .0083007 -3.17 0.002 -.0425457 -.0100074
/g31 | .0042753 .0043604 0.98 0.327 -.004271 .0128217
/g22 | .1215136 .0057199 21.24 0.000 .1103028 .1327244
/g32 | -.0522457 .0039323 -13.29 0.000 -.0599529 -.0445385
/g33 | -.0351974 .0021802 -16.14 0.000 -.0394705 -.0309243
/b1 | .0643824 .0044601 14.44 0.000 .0556407 .0731241
/b2 | -.001214 .0019786 -0.61 0.539 -.005092 .0026639
/b3 | .0425595 .0017629 24.14 0.000 .0391043 .0460146
/l1 | -.0025805 .0041987 -0.61 0.539 -.0108098 .0056489
/l2 | -.0000699 .0032462 -0.02 0.983 -.0064323 .0062925
/l3 | .0011486 .0017053 0.67 0.501 -.0021937 .0044908
/r1 | -.0005774 .0018703 -0.31 0.758 -.0042432 .0030883
/r2 | .0000853 .0045997 0.02 0.985 -.0089299 .0091006
/m11 | -.0005068 .0014367 -0.35 0.724 -.0033228 .0023092
/m12 | .003366 .0035336 0.95 0.341 -.0035598 .0102918
/m13 | -.0002668 .0007543 -0.35 0.724 -.0017452 .0012117
/m21 | -.0023503 .0018552 -1.27 0.205 -.0059865 .0012859
/m22 | 0 (constrained)
/m23 | 0 (constrained)
/d1 | 0 (constrained)
/d2 | 0 (constrained)
/d3 | 0 (constrained)
------------------------------------------------------------------------------

.
. est store quaidsNNP2

.
. *set trace on
. *set tracedepth 4
.
. *Share, price and total expenditure means
. quietly {

. 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]), post noheader

------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
a1 | .3145464 .0092701 33.93 0.000 .2963772 .3327155
a2 | .2708841 .0071607 37.83 0.000 .2568494 .2849187
a3 | .1063753 .003666 29.02 0.000 .0991901 .1135604
a4 | .3081943 .0068902 44.73 0.000 .2946897 .3216989
b1 | .0643824 .0044601 14.44 0.000 .0556407 .0731241
b2 | -.001214 .0019786 -0.61 0.539 -.005092 .0026639
b3 | .0425595 .0017629 24.14 0.000 .0391043 .0460146
b4 | -.1057278 .0045213 -23.38 0.000 -.1145894 -.0968662
g11 | .0098714 .0107752 0.92 0.360 -.0112477 .0309905
g12 | -.0262766 .0083007 -3.17 0.002 -.0425457 -.0100074
g13 | .0042753 .0043604 0.98 0.327 -.004271 .0128217
g21 | -.0262766 .0083007 -3.17 0.002 -.0425457 -.0100074
g22 | .1215136 .0057199 21.24 0.000 .1103028 .1327244
g23 | -.0522457 .0039323 -13.29 0.000 -.0599529 -.0445385
g31 | .0042753 .0043604 0.98 0.327 -.004271 .0128217
g32 | -.0522457 .0039323 -13.29 0.000 -.0599529 -.0445385
g33 | -.0351974 .0021802 -16.14 0.000 -.0394705 -.0309243
g14 | .0121298 .0086606 1.40 0.161 -.0048446 .0291042
g24 | -.0429913 .0094628 -4.54 0.000 -.0615381 -.0244446
g34 | .0831678 .0062559 13.29 0.000 .0709064 .0954292
g41 | .0121298 .0086606 1.40 0.161 -.0048446 .0291042
g42 | -.0429913 .0094628 -4.54 0.000 -.0615381 -.0244446
g43 | .0831678 .0062559 13.29 0.000 .0709064 .0954292
g44 | -.0523063 .0149227 -3.51 0.000 -.0815542 -.0230583
l1 | -.0025805 .0041987 -0.61 0.539 -.0108098 .0056489
l2 | -.0000699 .0032462 -0.02 0.983 -.0064323 .0062925
l3 | .0011486 .0017053 0.67 0.501 -.0021937 .0044908
l4 | .0015018 .0033825 0.44 0.657 -.0051278 .0081314
------------------------------------------------------------------------------

.
. forv i=1(1)3 {
2. forv j=1(1)3 {
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. est restore quaidsNNP2
8. qui nlcom (elasexp`i': ${mu`i'}/w`i'mean + 1) (mu`i'`j'f: (1e+2)*(${mu`i'`j'})), post noheader
9. qui nlcom (elasexp`i': _b[elasexp`i']) (mu`i'`j':_b[mu`i'`j'f]/(1e+2)), post noheader
10. }
11. * Uncompensated price elasticity
. nlcom (elasexp`i': _b[elasexp`i']) (elu`i'`j':_b[mu`i'`j']/w`i'mean - ${delta}) , post noheader
12. * Compensated price elasticity
. nlcom (elc`i'`j': _b[elu`i'`j'] + _b[elasexp`i']*w`j'mean), noheader
13. qui est restore quaidsNNP2
14. }
15. }
(results quaidsNNP2 are active now)
[b1] not found
r(111);

end of do-file

r(111);

[Jorge Eduardo Perez Perez]

Can you send me the code for the quaidsNNP routine? I can't test this code without it.

[Nigussie Tefera]

Dear Jorge,

Please find the code for quaidsNNP

*
program nlsurquaidsNNP

version 12.1

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

tokenize `varlist'
args w1 w2 w3 lnp1 lnp2 lnp3 lnp4 lnm x1 x2 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'*(`lnm' - `lnpindex') + ///
`l1'/`bofp'*(`lnm' - `lnpindex')^2 + ///
`r11'*`x1' +`r12'*`x2')*`cdf1'+ ///
`d1'*`pdf1'

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

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

end