*AIB Presentation 11/16/22
*This file is used to demonstrate how to use simulations. The first
*step specifies a data generation process (DGP) where y is a function of x's and e.
*The correlation matrix specifies the correlations between the x's and e. 
*Then, this runs a regression and saves the outputs. The second step uses the 
*-simulate- command to run the program across many iterations. 

*STEP 1*
clear
*setting the seed allows for the replication of simulation results
set seed 12345
*below we name the program "regsims"
program regsims, rclass
drop _all
*'set obs' allows us to set the sample size for each regression
set obs 500 
* below we set the correlation matrix for x1-e. notice there are 
* 4 rows/columns and 4 variables in the following line.
matrix regintro =  (1., .0, .5, .0 \ 			///		
					.0, 1., .2, .0 \				///
					.5, .2, 1., .0 \ 				///
					.0, .0, .0, 1.)				
*'drawnorm' refers to drawing randomly generated variables based on the 
* correlation matrix with the specified means and SDs of each variable.
drawnorm x1 x2 x3 e, corr(regintro) means(0 0 0 0) sds(1 1 1 1)
* the ID command gives each row an ID (e.g., like a firm identifier)
gen ID=_n
* here is our ultimate data generating process (DGP) that defines our y.
gen y = .2*x1 + .2*x2 + e

* with our dataset, we now run the regression and output a series of variables
* FOR EACH REGRESSION. We compile them later.
ivreg y x1 (x2 = x3)

*first we store stats for the coefficients
return scalar bx12sls= _b[x1] //this stores the beta for x1
return scalar sex12sls= _se[x1] //this stores the SE for x1
return scalar bx22sls= _b[x2] //this stores the beta for x2
return scalar sex22sls= _se[x2] //this stores the SE for x2

*run OLS regression and store estimates
reg y x1 x2
return scalar bx1ols= _b[x1] //this stores the beta for x1
return scalar sex1ols= _se[x1] //this stores the SE for x1
return scalar bx2ols= _b[x2] //this stores the beta for x2
return scalar sex2ols= _se[x2] //this stores the SE for x2


end

*STEP 2*
*now that we have the program set to create data and run regs, we can run sims. below this runs "reps" iterations and*
*collapses to give average bx1, sex1, and mr2 (model r2) for the conditions specified. 
simulate  bx1olsf=r(bx1ols) bx12slsf=r(bx12sls)  ///
bx2olsf=r(bx2ols) bx22slsf=r(bx22sls)  ///
sex1olsf=r(sex1ols) sex12slsf=r(sex12sls) ///
sex2olsf=r(sex2ols) sex22slsf=r(sex22sls)  ///
, reps(500) saving(regsimsdata, replace) nolegend nodots: regsims


*the next few lines summarize our findings
sum bx1olsf bx12slsf  ///
sex1olsf sex12slsf  ///
bx2olsf bx22slsf  ///
sex2olsf sex22slsf    

*this compares the betas
graph bar (mean) bx1olsf (mean) bx12slsf  (mean) bx2olsf (mean) bx22slsf

*this compares the SEs
graph bar (mean) sex1olsf (mean) sex12slsf (mean) sex2olsf (mean) sex22slsf 


program drop regsims