Value-Added Model: data simulation

vam_datasim.ado

*! version 1.0  27jul2013  Michael Stepner, stepner@mit.edu

program define vam_datasim
	
	version 11

	syntax , var_mu(real) var_w(real) delta(real) var_lambda(real) var_ind(real) num_teachers(integer) max_T(integer) max_C(integer) max_I(integer) ///
		[teacher(name) year(name) class(name) score(name) ///
		randomize_num_classes randomize_class_sizes randomize_class_years ///
		debug]

	*Note that in steady state with large T, var_w = var_psi/(1-delta^2)
	local var_psi=`var_w' * (1-`delta'^2)
	
	***** Default variable names *****
	
	if ("`teacher'"=="") local teacher="teacher"
	if ("`year'"=="") local year="year"
	if ("`class'"=="") local class="class"
	if ("`score'"=="") local score="score"
	
	**************************************
	***** Generate teacher component *****
	**************************************

	clear
	qui set obs `=`num_teachers'*`max_T''

	* Teacher id's
	g `teacher' = 1+floor((_n-1)/`max_T')

	* Years
	g `year' = _n - (`teacher'-1)*`max_T'

	sort `teacher' `year'

	*Generate persistent component of teacher VA
	qui g temp = sqrt(`var_mu')*rnormal() if `year'==1
	by `teacher': egen mu_j = min(temp)
	drop temp

	*Generate idiosyncatic shock to VA in each period
	g psi_jt = sqrt(`var_psi')*rnormal()

	*Set initial value of transitory shock
	qui g temp = sqrt(`var_psi')*rnormal() if `year'==1
	qui g w_jt = `delta'*psi_jt + temp
	drop temp

	*Propagate transitory shock forward
	forvalues s = 2/`max_T' {
		qui replace w_jt = `delta'*w_jt[_n-1] + psi_jt if `year'==`s' & `teacher'==`teacher'[_n-1] & `year'==`year'[_n-1]+1
	}

	*Calculate teacher effect in each period: sum of permanent and transitory component
	g tfx_jt = mu_j + w_jt

	tempfile teachers_sim
	qui save `teachers_sim'

	**********************************************
	***** Generate class component and merge *****
	**********************************************

	clear
	qui set obs `=`num_teachers'*`max_C'*`max_I''

	*Teacher id's
	g `teacher' = 1+floor((_n-1)/(`max_C'*`max_I'))

	*** Classes

	*Class id's
	g `class' = 1+floor((_n-1)/`max_I')

	if "`randomize_num_classes'"!="" {
		* Set random number of classes per teacher, distributed U[1,max_C]
		quietly {
			bys `teacher': g temp=1+int(`max_C'*runiform()) if _n==1
		}
		by `teacher': egen num_classes=min(temp)
		drop temp
		
		* Drop extra classes
		g temp=`class'-`max_C'*(`teacher'-1)
		qui drop if temp>num_classes
		drop temp
	}
	
	if ("`randomize_class_sizes'"!="") quietly {
		* Set random number of students per class, distributed N[max_I/2,max_I/3] and truncated at 1 and max_I.
		bys `teacher' `class': g temp=min(max(1,round(rnormal(`max_I'/2,`max_I'/3))),`max_I') if _n==1
		by `teacher' `class': egen classsize=min(temp)
		drop temp
		
		* Drop extra students
		by `teacher' `class': drop if _n>classsize
	}


	*** Years
	
	if "`randomize_class_years'"=="" {
		* Spread classes evenly between years.
		*   note: if max_C is not a multiple of max_T, then the number of classes taught per year will vary by 1 for each teacher.  the years with an extra class will vary across teachers.
		g `year'=mod(`class'+`max_T'-1,`max_T')+1
	}
	else quietly {
		* Random years
		bys `teacher' `class': g temp=1+int(`max_T'*runiform()) if _n==1
		by `teacher' `class': egen `year'=min(temp)
		drop temp
	}

	*** Merge in teacher-level shocks
	quietly merge m:1 `teacher' `year' using `teachers_sim', keep(3) nogen

	*** Create class/individual-level shocks

	* Simulate idiosyncatic class-level shocks
	quietly {
		bys `teacher' `class': g temp = sqrt(`var_lambda')*rnormal() if _n==1
	}
	by `teacher' `class': egen lambda_jt = min(temp)
	drop temp

	* Simulate individual shocks
	g ind_ijt= sqrt(`var_ind')*rnormal()

	* Calculate individuals' scores
	g `score' = tfx_jt + lambda_jt + ind_ijt


	*** Prepare final output

	if ("`debug'"=="debug") {	
		* Calculate class means
		by `teacher' `class': egen cfx_jt=mean(`score')
	}
	else keep `teacher' `year' `class' `score'
	
	* Order & sort
	order `teacher' `year' `class' `score'
	sort `teacher' `year' `class'

end