Simulate bias of random effects in TS model

random_effects_bias.do

* Simulation: Broken or Fixed Effects?

set matsize 5000
clear all

global graph_opts ///
	title(, justification(left) color(black) span pos(11)) ///
	graphregion(color(white) lc(white) lw(med) la(center)) /// <- Delete la(center) for version < 15
	ylab(,angle(0) nogrid) xtit(,placement(left) justification(left)) ///
	yscale(noline) xscale(noline) legend(region(lc(none) fc(none)))

* Run a simulation

	foreach sampleSize of numlist 1 2 3 4 {
		qui forvalues iteration = 1/50 {

			* Create data

				clear 
					// Need an empty dataset
				local nuids = 10^`sampleSize' 
					// Number of units
				set obs `nuids' 
					// Create units
				
				gen fe = rnormal() 
					// Gen fixed effects
				gen uid = _n 
					// Gen unit IDs
				expand 10 
					// Create panel structure
				bys uid: gen period = _n
					// Time structure
				
				gen x_vary = rnormal() + fe 
					// Time-variant characteristics correlated with FEs
				gen y = x_vary + fe + rnormal() 
					// True DGP includes correlated X's and FE's – violates RE assumption

			* Run a regression

				xtset uid
					// Prep for [xtreg]
			
				xtreg y x_vary
					// The biased regression
				local err_re = abs(1-_b[x_vary])
					// Store the error
			
				xtreg y x_vary , fe
					// The unbiased regression
				local err_fe = abs(1-_b[x_vary])
					// Store the error
			
				mat results = /// All data to matrix
					nullmat(results) \ [`err_re',`err_fe',`nuids']

		}
	}

* Visualize the simulation

	clear
	svmat results

	tw ///
		(lpolyci results1 results3 , lw(thick)) ///
		(lpolyci results2 results3 , lw(thick)) ///
	, ${graph_opts} ///
		xscale(log) xlab(10 100 1000 10000) ///
		xtitle("Observations{&rarr}") ///
		ytitle(" ") yline(0 , lw(thin)) ///
		legend(c(1) ring(0) pos(3) ///
			order(0 "Average Error:" 2 "Random Effects" 4 "Fixed Effects"))

* Have a lovely day!