mtefe_gendata.ado

*! 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