version 11 clear clear matrix clear mata set more off set mem 500m capture log close /*-----------------------------------------------------* * Title - Elasticity Estimator Simulation * * Authors - Matt Trombley/Joe Terza * * Date Created - 7/12/11 * * Purpose - Simulation to determine the difference * * between true and estimated elasticity * *------------------------------------------------------*/ /************************************************* ** Set up the output file. ** *************************************************/ *Opening MATA mata /*-----------------------------------------------------* * Generating x-variables and default coefficient values* * per Terza presentation on blackboard. * *------------------------------------------------------*/ /***************************************************** ** Choose number of repititions and number of ** ** observations. ** *****************************************************/ reps = 1000 obs = 25000 trueobs=1000000 /***************************************************** ** Set the true parameter values. ** *****************************************************/ betap1 = -0.7518 betap2 = -0.0573 beta1oconst= 3.0358 beta2oconst= 0.9756 /* beta1oconst=1 beta2oconst=2 */ beta1o = -0.482 \ beta1oconst beta2o = -1.997 \ beta2oconst /***************************************************** ** Pick the means and the variances of the x's. ** *****************************************************/ xomean = 3 xovar = 2 xpmean = 4 xpvar = 3 /***************************************************** ** Generate the x's for the calculation of the true ** ** elasticity value. ** *****************************************************/ xoc1=sqrt(12:*xovar) xoc2=xomean:-(xoc1:/2) xo = xoc1:*runiform(trueobs,cols(xomean)):+xoc2 xpc1=sqrt(12:*xpvar) xpc2=xpmean:-(xpc1:/2) xp = ln(xpc1:*runiform(trueobs,cols(xpmean)):+xpc2) /***************************************************** ** Create the xo matrix, with the constant term for ** ** the calculation of the true elasticity value. ** *****************************************************/ xoc=xo,J(trueobs,1,0) /***************************************************** ** Create the x matrix, with xp and NO constant ** ** term for the calculation of the true ** ** elasticity value. ** *****************************************************/ xpxo=xp,xo /***************************************************** ** Create the x matrix, with xp and the constant ** ** term for the calculation of the true ** ** elasticity value. ** *****************************************************/ xpxoc=xp,xoc /***************************************************** ** Concatenate the true vector of parameters for the** ** calculation of the true elasticity value. ** *****************************************************/ truebeta1=betap1 \ beta1o truebeta2=betap2 \ beta2o /************************************************* ** Compute x1b1 and x1b2 multiplying the matrix ** ** of exogenous variables (xpxoc) by the ** ** coefficient vectors for the ** ** calculation of the true elasticity value. ** *************************************************/ x1b1=xpxoc*truebeta1 x2b2=xpxoc*truebeta2 /************************************************* ** Compute the true elasticity measure with ** ** log price. ** *************************************************/ numlog=mean(invlogit(x1b1):*(1:-invlogit(x1b1)):*exp(x2b2):*betap1) denomlog=mean(invlogit(x1b1):*exp(x2b2)) etatrue=numlog/denomlog:+betap2 /***************************************************** ** Generate the x's for the simulation loop. ** *****************************************************/ xoc1=sqrt(12:*xovar) xoc2=xomean:-(xoc1:/2) xo = xoc1:*runiform(obs,cols(xomean)):+xoc2 xpc1=sqrt(12:*xpvar) xpc2=xpmean:-(xpc1:/2) xp = ln(xpc1:*runiform(obs,cols(xpmean)):+xpc2) /***************************************************** ** Create the xo matrix, with the constant term. ** *****************************************************/ xoc=xo,J(obs,1,0) /***************************************************** ** Create the x matrix, with xp and NO constant ** ** term. ** *****************************************************/ xpxo=xp,xo /***************************************************** ** Create the x matrix, with xp and the constant ** ** term. ** *****************************************************/ xpxoc=xp,xoc /***************************************************** ** Appended xp to the list of x variables (x1 with ** ** last variable now changed to x21). ** *****************************************************/ xvars = st_addvar("float",("x1", "x2")) /***************************************************** ** Set the number of observations to be included in ** ** the temporary Stata dataset. ** *****************************************************/ st_addobs(obs) /***************************************************** ** xpxo is now the matrix of x's. ** *****************************************************/ st_store(.,xvars,xpxo) /*----------------------------------------------------* * Initializing sample variable, and creating empty * * matrix for results * *-----------------------------------------------------*/ samp = 0 results = J(reps,4,.) /*----------------------------------------------------* * Initializing do-loop * *-----------------------------------------------------*/ do{ samp = samp + 1 /*-----------------------------------------------------* * Generating data for first stage (logit) and second * * stage (log-OLS) * *------------------------------------------------------*/ /***************************************************** ** Generate h. ** *****************************************************/ /***************************************************** ** Used xo matrix with the constant term, i.e. xoc. ** *****************************************************/ /* expLinIndex = exp(-(xp*betap1 + xoc*beta1o)) logisticPr = expLinIndex :/ (1 :+ expLinIndex) h = runiform(obs, 1) :> logisticPr */ h=(xp*betap1 :+ xoc*beta1o:+logit(runiform(obs, 1))):>0 /***************************************************** ** Generate generate y. ** *****************************************************/ /***************************************************** ** Used xo matrix with the constant term, i.e. xoc. ** *****************************************************/ sigma = 5 signorm = sigma*rnormal(obs,1,0,1) y = exp(xp*betap2 + xoc*beta2o + signorm) /***************************************************** ** Send h and y to Stata. ** *****************************************************/ hvar = st_addvar("float",("h")) st_store(.,hvar,h) yvar = st_addvar("float",("y")) st_store(.,yvar,y) /*-----------------------------------------------------* * Compute first stage logit regression and save results* * in Mata. * *------------------------------------------------------*/ stata("qui logit h x*") beta1=st_matrix("e(b)")' bp1=beta1[1] /*-----------------------------------------------------* * Compute second stage ols regression and save results * * in Mata. * *------------------------------------------------------*/ stata("gen y2 = log(y)") /***************************************************** ** Set y2 missing if h=0 ** *****************************************************/ stata("replace y2=. if h==0") stata("qui reg y2 x*") beta2=st_matrix("e(b)")' bp2=beta2[1] /************************************************* ** Compute x1b1 and x1b2 multiplying the matrix ** ** of exogenous variables (xpxoc) by the ** ** coefficient vectors. ** *************************************************/ x1b1=xpxoc*beta1 x2b2=xpxoc*beta2 /************************************************* ** Compute the correct elasticity measure with ** ** log price. ** *************************************************/ numlog=mean(invlogit(x1b1):*(1:-invlogit(x1b1)):*exp(x2b2):*bp1) denomlog=mean(invlogit(x1b1):*exp(x2b2)) etahat=numlog/denomlog:+bp2 /************************************************* ** Compute the conventional (Manning et al. 1995)* ** elasticity measure. ** *************************************************/ etatilde=mean((1:-mean(h)):*bp1:+bp2) /************************************************* ** Compute the proportional absolute biases. ** *************************************************/ abspbiashat=abs(etahat-etatrue)/abs(etatrue) abspbiastilde=abs(etatilde-etatrue)/abs(etatrue) stata("drop h y y2") /************************************************* ** Save the respective results in the "results" ** ** matrix. ** *************************************************/ etahat,etatilde,abspbiashat,abspbiastilde results[samp,1] = etahat results[samp,2] = etatilde results[samp,3] = abspbiashat results[samp,4] = abspbiastilde /************************************************* ** Report the sample number. ** *************************************************/ "This is sample #" samp } while (samp