// Monte Carlo exercise to show sampling distribution of s2
clear all
global sample_size=50
global num_regressors=3
global sigma2=1

// make sample_size observations, with regressors xk, each centered at k, having a beta of k
global regressors ""
range count 1 100000 100000
g e=rnormal(0,$sigma2)
g y=e
forval k=1/$num_regressors {
	g x`k'=rnormal(`k',1)
	replace y=y+x`k'*`k'
	global regressors "$regressors x`k'"
}
summ y $regressors e

regress y $regressors

preserve
sample $sample_size, count
regress y $regressors
matrix coefs_keep=e(rss),$sample_size,e(b)
restore

forval i=1/19 {
	preserve
	sample $sample_size, count
	qui regress y $regressors
	matrix coefs_keep=coefs_keep\(e(rss),$sample_size,e(b))
	restore
}

svmat coefs_keep
g s2=coefs_keep1/(coefs_keep2-1-$num_regressors)
summ
kdensity s2