*! mtefe_gendata v1.82 28may2018 * Author: Martin Eckhoff Andresen * This program is part of the mtefe package. cap program drop mtefe_gendata { program define mtefe_gendata, rclass version 13.0 syntax , [parameters(namelist min=4 max=5) fixedeffects(namelist min=3 max=3) /* */ obs(integer 5000) POLynomial(integer 0) link(string) districts(integer 0)] qui { clear //Check input if "`link'"!="" { if !inlist("`link'","lpm","probit") { noi di in red "Link() must be probit or lpm." exit 301 } } else loc link probit tempname U0 U1 V index gamma beta0 beta1 pi pi0 pi1 if "`parameters'"!="" { gettoken gammaname rest: parameters confirm matrix `gammaname' mat `gamma'=`gammaname' gettoken beta0name rest: rest confirm matrix `beta0name' mat `beta0'=`beta0name' gettoken beta1name rest: rest confirm matrix `beta1name' mat `beta1'=`beta1name' if colsof(`gamma')!=colsof(`beta1')+1|colsof(`beta1')!=colsof(`beta0')|colsof(`beta1')<`=`numX'+1' { noi di in red "Matrix beta0, beta1 or gamma of wrong dimensions - should be 1 x K for betas and 1 x K+1 for gamma." exit 301 } local numpi: word count `rest' if `numpi'==1&`polynomial'!=0|`numpi'==2&`polynomial'<=0 { noi di in red "Specify either variance-covariance matrix for joint normal errors or two matrices for the pi-parameters if using the polynomial model." exit 301 } if `polynomial'==0 { confirm matrix `rest' mat `pi'=`rest' if colsof(`pi')!=3|rowsof(`pi')!=3 { noi di in red "When using joint normal errors, matrix pi must be 3 x 3." exit 301 } } if `polynomial'!=0 { tempname pi0name pi1name gettoken pi0name pi1name: rest confirm matrix `pi0name' `pi1name' mat `pi0'=`pi0name' mat `pi1'=`pi1name' if colsof(`pi1')!=`polynomial'|colsof(`pi0')!=`polynomial' { noi di in red "Columns of pi1 and pi0 must match degree of polynomial." exit 301 } } } else { if !inlist(`polynomial',0,2) { noi di in red "When using default parameters, use polynomial of degree 0 (normal) or 2." exit 301 } if "`link'"=="lpm" mat `gamma'=[ 0.015 , -0.01 , 0.00025 , `=0.5-0.015*40+0.01*15-0.00025*305'] else mat `gamma'=[ -0.125 , -0.08 , 0.002 , `=0.125*40+0.08*15-0.002*305'] mat `beta0'=[ 0.025 , -0.0004 , 3.2] mat `beta1'=[ 0.01 , 0 , 3.6 ] mat `pi'= [ 0.5 , 0.3 , -0.1 \ 0.3 , 0.5, -0.5 \ -0.1 , -0.5 , 1 ] mat `pi1'= [ 0.5 , -0.1 ] mat `pi0'= [ 2 , -1 ] } if "`fixedeffects'"!="" { tempname fix0 fix1 avgDist gettoken fix0name rest: fixedeffects confirm matrix `fix0name' mat `fix0'=`fix0name' loc districts=colsof(`fix0') gettoken fix1name rest: rest confirm matrix `fix1name' mat `fix1'=`fix1name' if colsof(`fix1')!=`districts' { noi di in red "Fixed effects matrices must be same dimensions - 1 x G." exit 301 } confirm matrix `rest' mat `avgDist'=`rest' if colsof(`avgDist')!=`districts' { noi di in red "Fixed effects matrices must be same dimensions - 1 x G." exit 301 } } else if `districts'>0 { set obs `districts' tempname fix0 fix1 avgDist cov means mat `cov'= [ 0.1 , 0.05 , -0.05 \ 0.05 , 0.1 , -0.1 \ -0.05 , -0.1 , 10 ] drawnorm `fix0' `fix1' `avgDist', cov(`cov') foreach state in fix0 fix1 avgDist { replace ``state''=0 in 1 mkmat ``state'', matrix(``state'') mat ``state''=``state''' } } //Generate X, Z clear set obs `obs' if `districts'>0 { tempname fix gen district=ceil(`districts'*runiform()) fvexpand i.district loc districtlist `r(varlist)' foreach state in fix0 fix1 avgDist { mat colnames ``state''=`districtlist' } mat score `avgDist'=`avgDist' gen distCol=`avgDist'+rnormal(40,10) } else gen distCol=rnormal(0,10)+40 gen exp=floor(31*runiform()) gen exp2=exp^2 //Combine coefs, include fixed effects if `districts'>0 { foreach state in 0 1 { mat `beta`state''=`beta`state''[1,1..colsof(`beta`state'')-1],`fix`state'',`beta`state''[1,colsof(`beta`state'')] } } mat colnames `gamma'=distCol exp exp2 _cons mat colnames `beta0'=exp exp2 `districtlist' _cons mat colnames `beta1'=exp exp2 `districtlist' _cons //Create errors tempname U0 U1 V Ud if `polynomial'==0 { drawnorm `U0' `U1' `V', cov(`pi') gen `Ud'=normal(`V') } if `polynomial'>0 { gen `V'=rnormal() gen `Ud'=normal(`V') loc cons1=0 loc cons0=0 forvalues p=1/`polynomial' { tempname p`p' loc names `names' `p`p'' loc cons1=`cons1'-`pi1'[1,`p']*(1/(`p'+1)) loc cons0=`cons0'-`pi0'[1,`p']*(1/(`p'+1)) gen `p`p''=`Ud'^`p' } mat `pi1'=`pi1',`cons1' mat `pi0'=`pi0',`cons0' mat colnames `pi1'=`names' _cons mat colnames `pi0'=`names' _cons mat score `U1'=`pi1' mat score `U0'=`pi0' replace `U1'=`U1'+rnormal(0,0.2) replace `U0'=`U0'+rnormal(0,0.2) } //determine first stage tempname index mat score `index'=`gamma' if "`link'"=="probit" gen col=`index'>`V' if "`link'"=="lpm" gen col=`index'>`Ud' //Determine outomes tempname lwage0 lwage1 mat score `lwage0'=`beta0' mat score `lwage1'=`beta1' replace `lwage0'=`lwage0'+`U0' replace `lwage1'=`lwage1'+`U1' gen lwage=col*`lwage1'+(1-col)*`lwage0' //generate MTE tempname support temp mtexs beta10pi mte forvalues i=1/99 { mat `support'=[nullmat(`support'),`=round(`i'/100,0.01)' ] loc mtenames `mtenames' u`i' } mat `mtexs'=[ 15 \ 305 ] if `districts'>0 { forvalues i=1/`districts' { mat `mtexs'=[nullmat(`mtexs') \ `=1/`districts'' ] } } mat `mtexs'=[nullmat(`mtexs') \ 1] if `polynomial'==0 { mat `beta10pi'=`beta1'-`beta0',`pi'[3,2]-`pi'[3,1] mata: mtefecalc_small(st_matrix("`beta10pi'"),invnormal(st_matrix("`support'")),st_matrix("`mtexs'"),"`mte'") } else { tempname S tempsup forvalues k=1/`polynomial' { if `k'==1 mat `tempsup'=`support' else mat `tempsup'=hadamard(`tempsup',`support') mat `S'=[nullmat(`S') \ `tempsup'] } mat `beta10pi'=`beta1'[1,1..`=colsof(`beta1')-1']-`beta0'[1,1..`=colsof(`beta0')-1'],`beta1'[1,`=colsof(`beta1')']-`beta0'[1,`=colsof(`beta0')']+`pi1'[1,`=`polynomial'+1']-`pi0'[1,`=`polynomial'+1'],`pi1'[1,1..`polynomial']-`pi0'[1,1..`polynomial'] mata: mtefecalc_small(st_matrix("`beta10pi'"),st_matrix("`S'"),st_matrix("`mtexs'"),"`mte'") } mat colnames `mte'=`mtenames' //post parameters return matrix gamma=`gamma' return matrix beta1=`beta1' return matrix beta0=`beta0' return matrix mte=`mte' if `polynomial'>0 { forvalues state=0/1 { tempname retpi`state' mat `retpi`state''=`pi`state''[1,1..`polynomial'] return matrix pi`state'=`retpi`state'' } } else return matrix pi=`pi' /*gen lwage0=`lwage0' gen lwage1=`lwage1' gen U1=`U1' gen U0=`U0' gen Ud=`Ud' */ } end } mata: mata clear void mtefecalc_small(real matrix beta10pi, real matrix S, real matrix mtexs, /// string scalar mtename) { real matrix fullS real scalar i real vector mte fullS=J(rows(mtexs),cols(S),.) for (i=1;i<=cols(S);i++) fullS[.,i]=mtexs fullS=fullS \ S mte=beta10pi*fullS st_matrix(mtename,mte) } end