** Program ass4w07q4.do for Stata version 9
** STATA Program by A. Colin Cameron 

********** SETUP
clear
capture log close 
log using ass4w07q4.txt, text replace
set mem 10m
set more off
version 9.0
set scheme s1mono  /* Graphics scheme */

********** OVERVIEW OF ASS4W07Q4.DO 

* This program does Monte Carlo study for OLS estimator with bivariate regression
* where y is redrawn each time while x is fixed.
* - No data need to be read in
* - Sample size N is defined in the global numobs
* - Number of simulation replications in defined in the global numsims

********** MONTE CARLO OVERVIEW

* The data generating process is 
* - y = b1 + b2*x2 + u 
* - where b1 = 10 and b2 = 5
* - where regressor x = 1, 2, 3, ..., N where N is sanmple size
* - and where u is either 
*     - N[0, 2^2] 
*     - chisquare(2) - 2 

* The simulation is done using stata command simulate

********** INITIAL SIMULATION SET UP

set seed 10101         /* set seed so get same results in future */
global numobs "5"      /* sample size N */
global numsims "1000"  /* number of simulations */

****** PROGRAM TO BE SIMULATED 

* The program is rclass, so results returned by program are put into r( )
* Here we return b2 and se2 and t2

* The program has no arguments (makes program simpler).
* It returns one value: 
*    - the slope coeff
* labelled b2

program simols, rclass
    version 9.0  
    /* No arguments to define in this example */ 
    /* Generate the data used in this example   */ 
    drop _all
    set obs $numobs
    gen x = _n                                 // x is a time trend
    gen y = 10 + 5*x + 2*invnorm(uniform())    // error is N[0,2^2]
    gen z = 10 + 5*x + (invnorm(uniform()))^2 + (invnorm(uniform()))^2 -2
                                               // error is chisquare(2)-2
    /* Do the analysis for this example */
    regress y x
    return scalar b2normal =_b[x]
    regress z x
    return scalar b2chisquare =_b[x]
end

****** RUN PROGRAM ONCE TO CHECK

simols
sum y z x
sum

****** RUN THE SIMULATIONS

simulate "simols" b2normal=r(b2normal) b2chisquare=r(b2chisquare), reps($numsims) 

****** SUMMARIZE SIMULATION RESULTS

summarize
summarize, detail

kdensity b2normal, normal
graph export b2normal.wmf, replace

kdensity b2chisquare, normal
graph export b2chisquare.wmf, replace

********** CLOSE OUTPUT
log close