* Tricks with Hicks: The EASI demand system * Arthur Lewbel and Krishna Pendakur * 2008, American Economic Review * Herein, find Stata code to estimate a demand system with neq equations, nprice prices, * ndem demographic characteristics and npowers powers of implicit utility * This Stata code estimates Lewbel and Pendakur's EASI demand system using approximate * OLS estimation and iterated linear 3SLS estimation. Note that iterated linear 3SLS is * not formally equivalent to fully nonlinear 3SLS (which does not exist in Stata). * However, in some contexts they are asymptotically equivalent (see, e.g., Blundell and * Robin 1999 and Dominitz and Sherman 2005), and we have verified in our data that * coefficients estimated using iterated linear 3SLS are within 0.001 of those * estimated using fully nonlinear 3SLS. * Code to estimate the fully nonlinear 3SLS/GMM version in TSP is available on request * from the authors. * This model includes pz,py,zy interactions. See 'iterated 3sls without pz,py,zy.do' for * shorter code to estimate the model without interactions. set more off macro drop _all use "C:\projects\hixtrix\revision\hixdata.dta", clear * set number of equations and prices and demographic characteristics and convergence criterion global neqminus1 "7" global neq "8" global nprice "9" global ndem 5 global npowers "5" * set a convergence criterion and choose whether or not to base it on parameters global conv_crit "0.00000000000001" scalar conv_param=1 scalar conv_y=0 *note set the matrix size big enough to do constant,y,z,p,zp,yp,yz global matsize_value=100+\$neq*(1+\$npowers+\$ndem+\$neq*(1+\$ndem+1)+\$ndem) set matsize \$matsize_value *data labeling conventions: * data weights: wgt (replace with 1 if unweighted estimation is desired) * budget shares: s1 to sneq * prices: p1 to nprice * log total expenditures: x * implicit utility: y, or related names * demographic characteristics: z1 to zndem g obs_weight=wgt g s1=sfoodh g s2=sfoodr g s3=srent g s4=soper g s5=sfurn g s6=scloth g s7=stranop g s8=srecr g s9=spers g p1=pfoodh g p2=pfoodr g p3=prent g p4=poper g p5=pfurn g p6=pcloth g p7=ptranop g p8=precr g p9=ppers * polynomial systems are easier to estimate if you normalise the variable in the polynomial g x=log_y *egen mean_log_y=mean(log_y) *replace x=log_y-mean_log_y * normalised prices are what enter the demand system * generate normalised prices, backup prices (they get deleted), and pAp, pBp global nplist "" forvalues j=1(1)\$neq { g np`j'=p`j'-p\$nprice global nplist "\$nplist np`j'" } forvalues j=1(1)\$neq { g np`j'_backup=np`j' g Ap`j'=0 g Bp`j'=0 } g pAp=0 g pBp=0 *list demographic characteristics: fill them in, and add them to zlist below g z1=age g z2=hsex g z3=carown g z4=tran g z5=time global zlist "z1 z2 z3 z4 z5" *make pz interactions global npzlist "" forvalues j=1(1)\$neq { forvalues k=1(1)\$ndem { g np`j'z`k'=np`j'*z`k' global npzlist "\$npzlist np`j'z`k'" } } *make y_stone=x-p'w, and gross instrument, y_tilda=x-p'w^bar g y_stone=x g y_tilda=x forvalues num=1(1)\$nprice { egen mean_s`num'=mean(s`num') replace y_tilda=y_tilda-mean_s`num'*p`num' replace y_stone=y_stone-s`num'*p`num' } * make list of functions of (implicit) utility, y: fill them in, and add them to ylist below * alternatively, fill ylist and yinstlist with the appropriate variables and instruments g y=y_stone g y_inst=y_tilda global ylist "" global yinstlist "" global yzlist "" global yzinstlist "" global ynplist "" global ynpinstlist "" forvalues j=1(1)\$npowers { g y`j'=y^`j' g y`j'_inst=y_inst^`j' global ylist "\$ylist y`j'" global yinstlist "\$yinstlist y`j'_inst" } forvalues k=1(1)\$ndem { g yz`k'=y*z`k' g yz`k'_inst=y_inst*z`k' global yzlist "\$yzlist yz`k'" global yzinstlist "\$yzinstlist yz`k'_inst" } forvalues k=1(1)\$neq { g ynp`k'=y*np`k' g ynp`k'_inst=y_inst*np`k' global ynplist "\$ynplist ynp`k'" global ynpinstlist "\$ynpinstlist ynp`k'_inst" } *set up the equations and put them in a list global eqlist "" forvalues num=1(1)\$neq { global eq`num' "(s`num' \$ylist \$zlist \$yzlist \$nplist \$ynplist \$npzlist)" macro list eq`num' global eqlist "\$eqlist \\$eq`num'" } *create linear constraints and put them in a list, called conlist global conlist "" forvalues j=1(1)\$neq { local jplus1=`j'+1 forvalues k=`jplus1'(1)\$neq { constraint `j'`k' [s`j']np`k'=[s`k']np`j' global conlist "\$conlist `j'`k'" } } *add constraints for yp interactions forvalues j=1(1)\$neq { local jplus1=`j'+1 forvalues k=`jplus1'(1)\$neq { constraint `j'`k'0 [s`j']ynp`k'=[s`k']ynp`j' global conlist "\$conlist `j'`k'0" } } * add constraints for pz interactions forvalues h=1(1)\$ndem { forvalues j=1(1)\$neq { local jplus1=`j'+1 forvalues k=`jplus1'(1)\$neq { constraint `j'`k'`h' [s`j']np`k'z`h'=[s`k']np`j'z`h' global conlist "\$conlist `j'`k'`h'" } } } *an approximate model would use one of: *reg3 \$eqlist [aweight=obs_weight], constr(\$conlist) endog(\$ylist \$ynplist \$yzlist) exog(\$yinstlist \$ynpinstlist \$yzinstlist) *sureg \$eqlist, constr(\$conlist) *sureg \$eqlist *the exact model requires two steps: step 1) get a pre-estimate to construct the intrument, step 2) use the instrument to estimate the model *first get a pre-estimate to create the instrument: *run three stage least squares on the model with py, pz or yz interactions, and then iterate to convergence, * constructing y=(y_stone+0.5*p'A(z)p)/(1-0.5*p'Bp) at each iteration * note that the difference in predicted values for y=1 between p and p=0 is A(z)p, and * that the difference in difference in predicted values for y=1 vs y=0 between p and p=0 is Bp replace y=y_stone g y_backup=y_stone g y_old=y_stone g y_change=0 scalar crit_test=1 scalar iter=0 while crit_test>\$conv_crit { scalar iter=iter+1 quietly reg3 \$eqlist [aweight=obs_weight], constr(\$conlist) endog(\$ylist \$ynplist \$yzlist) exog(\$yinstlist \$ynpinstlist \$yzinstlist) if (iter>1) { matrix params_old=params } matrix params=e(b) quietly replace pAp=0 quietly replace pBp=0 quietly replace y_old=y quietly replace y_backup=y *predict with y=1 *generate rhs vars,interactions with y=1 forvalues j=1(1)\$npowers { quietly replace y`j'=1 } forvalues j=1(1)\$neq { quietly replace ynp`j'=np`j' } forvalues j=1(1)\$ndem { quietly replace yz`j'=z`j' } *generate predicted values forvalues j=1(1)\$neq { quietly predict s`j'hat_y1, equation(s`j') } *set all p, pz, py to zero foreach yvar in \$nplist \$ynplist \$npzlist { quietly replace `yvar'=0 } forvalues j=1(1)\$neq { quietly predict s`j'hat_y1_p0, equation(s`j') } *refresh p,pz forvalues j=1(1)\$neq { quietly replace np`j'=np`j'_backup forvalues k=1(1)\$ndem { quietly replace np`j'z`k'=np`j'_backup*z`k' } } *generate rhs vars,interactions with y=0 foreach yvar in \$ylist \$ynplist \$yzlist { quietly replace `yvar'=0 } *generate predicted values forvalues j=1(1)\$neq { quietly predict s`j'hat_y0, equation(s`j') } *set all p, pz, py to zero foreach yvar in \$nplist \$ynplist \$npzlist { quietly replace `yvar'=0 } forvalues j=1(1)\$neq { quietly predict s`j'hat_y0_p0, equation(s`j') } *refresh p only forvalues j=1(1)\$neq { quietly replace np`j'=np`j'_backup } *fill in pAp and pBp forvalues j=1(1)\$neq { quietly replace Ap`j'=s`j'hat_y0-s`j'hat_y0_p0 quietly replace pAp=pAp+np`j'*Ap`j' quietly replace Bp`j'=(s`j'hat_y1-s`j'hat_y1_p0)-(s`j'hat_y0-s`j'hat_y0_p0) quietly replace pBp=pBp+np`j'*Bp`j' quietly drop s`j'hat_y0 s`j'hat_y0_p0 s`j'hat_y1 s`j'hat_y1_p0 } *round pAp and pBp to the nearest millionth, for easier checking quietly replace pAp=int(1000000*pAp+0.5)/1000000 quietly replace pBp=int(1000000*pBp+0.5)/1000000 *recalculate y,yz,py,pz quietly replace y=(y_stone+0.5*pAp)/(1-0.5*pBp) forvalues j=1(1)\$npowers { quietly replace y`j'=y^`j' } forvalues j=1(1)\$ndem { quietly replace yz`j'=y*z`j' } *refresh py,pz forvalues j=1(1)\$neq { quietly replace ynp`j'=y*np`j'_backup forvalues k=1(1)\$ndem { quietly replace np`j'z`k'=np`j'_backup*z`k' } } if (iter>1 & conv_param==1) { matrix params_change=(params-params_old) matrix crit_test_mat=(params_change*(params_change')) svmat crit_test_mat, names(temp) scalar crit_test=temp drop temp } quietly replace y_change=abs(y-y_old) quietly summ y_change if(conv_y==1) { scalar crit_test=r(max) } display "iteration " iter scalar list crit_test summ y_change y_stone y y_old pAp pBp } *now, create the instrument, and its interactions yp and yz quietly replace y_inst=(y_tilda+0.5*pAp)/(1-0.5*pBp) forvalues j=1(1)\$npowers { quietly replace y`j'_inst=y_inst^`j' } forvalues j=1(1)\$neq { replace ynp`j'_inst=y_inst*np`j' } forvalues j=1(1)\$ndem { replace yz`j'_inst=y_inst*z`j' } *with nice instrument in hand, run three stage least squares on the model, and then iterate to convergence replace y_old=y replace y_change=0 scalar iter=0 scalar crit_test=1 while crit_test>\$conv_crit { scalar iter=iter+1 quietly reg3 \$eqlist [aweight=obs_weight], constr(\$conlist) endog(\$ylist \$ynplist \$yzlist) exog(\$yinstlist \$ynpinstlist \$yzinstlist) if (iter>1) { matrix params_old=params } matrix params=e(b) quietly replace pAp=0 quietly replace pBp=0 quietly replace y_old=y quietly replace y_backup=y *predict with y=1 *generate rhs vars,interactions with y=1 forvalues j=1(1)\$npowers { quietly replace y`j'=1 } forvalues j=1(1)\$neq { quietly replace ynp`j'=np`j' } forvalues j=1(1)\$ndem { quietly replace yz`j'=z`j' } *generate predicted values forvalues j=1(1)\$neq { quietly predict s`j'hat_y1, equation(s`j') } *set all p, pz, py to zero foreach yvar in \$nplist \$ynplist \$npzlist { quietly replace `yvar'=0 } forvalues j=1(1)\$neq { quietly predict s`j'hat_y1_p0, equation(s`j') } *refresh p,pz forvalues j=1(1)\$neq { quietly replace np`j'=np`j'_backup forvalues k=1(1)\$ndem { quietly replace np`j'z`k'=np`j'_backup*z`k' } } *generate rhs vars,interactions with y=0 foreach yvar in \$ylist \$ynplist \$yzlist { quietly replace `yvar'=0 } *generate predicted values forvalues j=1(1)\$neq { quietly predict s`j'hat_y0, equation(s`j') } *set all p, pz, py to zero foreach yvar in \$nplist \$ynplist \$npzlist { quietly replace `yvar'=0 } forvalues j=1(1)\$neq { quietly predict s`j'hat_y0_p0, equation(s`j') } *refresh p only forvalues j=1(1)\$neq { quietly replace np`j'=np`j'_backup } *fill in pAp and pBp forvalues j=1(1)\$neq { quietly replace Ap`j'=s`j'hat_y0-s`j'hat_y0_p0 quietly replace pAp=pAp+np`j'*Ap`j' quietly replace Bp`j'=(s`j'hat_y1-s`j'hat_y1_p0)-(s`j'hat_y0-s`j'hat_y0_p0) quietly replace pBp=pBp+np`j'*Bp`j' quietly drop s`j'hat_y0 s`j'hat_y0_p0 s`j'hat_y1 s`j'hat_y1_p0 } *round pAp and pBp to the nearest millionth, for easier checking quietly replace pAp=int(1000000*pAp+0.5)/1000000 quietly replace pBp=int(1000000*pBp+0.5)/1000000 *recalculate y,yz,py,pz quietly replace y=(y_stone+0.5*pAp)/(1-0.5*pBp) forvalues j=1(1)\$npowers { quietly replace y`j'=y^`j' } forvalues j=1(1)\$ndem { quietly replace yz`j'=y*z`j' } *refresh py,pz forvalues j=1(1)\$neq { quietly replace ynp`j'=y*np`j'_backup forvalues k=1(1)\$ndem { quietly replace np`j'z`k'=np`j'_backup*z`k' } } if (iter>1 & conv_param==1) { matrix params_change=(params-params_old) matrix crit_test_mat=(params_change*(params_change')) svmat crit_test_mat, names(temp) scalar crit_test=temp drop temp } quietly replace y_change=abs(y-y_old) quietly summ y_change if(conv_y==1) { scalar crit_test=r(max) } display "iteration " iter scalar list crit_test summ y_change y_stone y y_old pAp pBp } *note that reported standard errors are wrong for iterated estimates reg3 \$eqlist [aweight=obs_weight], constr(\$conlist) endog(\$ylist \$ynplist \$yzlist) exog(\$yinstlist \$ynpinstlist \$yzinstlist)