*********************************************************

clear

*** Set sample size to be drawn, means and standard deviations for the X and error terms,
*** and the number of simulation replications
global SampSize 10
global MeanX 0
global MeanError 0
global StdDevX 1
global StdDevError .3
local RepNum 20

***  NO NEED TO CHANGE ANYTHING BELOW UNLESS YOU KNWO WHAT YOU ARE DOING  ***

drop _all	

*** Define program simttest
capture program drop simttest
program define simttest, rclass
	quietly {
		drawnorm x, n($SampSize) means($MeanX) sds($StdDevX)
		drawnorm error, n($SampSize) means($MeanError) sds($StdDevError)
		gen y = 1 + x + error

		*** uncommenting the following two lines to allow measurement error on x
*		drawnorm merror, n($SampSize) means(0) sds(.8)
*		replace x = x + merror

			*** linear regression
			regress y x
				return scalar beta0 = _b[_cons]
				return scalar beta1 = _b[x]
				return scalar SE0 = _se[_cons]
				return scalar SE1 = _se[x]
	}
	drop *
end

*** Start simulation ***
simulate "simttest" beta0=r(beta0) beta1=r(beta1), reps(`RepNum') dots

gen yatminus1 = beta0 - beta1
gen yatplus1 = beta0 + beta1
mkmat yatminus1 yatplus1, matrix(A)
matrix B=A'
svmat B
gen x = -1 in 1
replace x = 1 in 2
gen y = x + 1

local ypos = `RepNum'+1
twoway scatter B* x, c(l ..) ms(i ..) clcolor(cyan ..)  clwidth(thin) || ///
	scatter y x, c(l) ms(i) clcolor(black) clwidth(thick) ///
	title("y = 1 + x + error") xlabel(#10,labsize(small)) ///
	xtitle(x, size(large)) ylabel(#10, labsize(small)) ///
	ytitle(y, orientation(horizontal) size(large)) ///
	subtitle((Sample size for each blue line = $SampSize), size(small)) ///
	legend(order(1 `ypos') label(1 "fitted lines") label(`ypos' "real line"))