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 3rd 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
}
}
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 3rd 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
}
*
*
*
*
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
}
}
