cscript 

set seed 112 

** Simulation parameters 

local N     = 1000    // Number of observations 
local R     = 2000    // Number of simulation repetitions

** True regression parameters 

local tb1 = 1.0

program define mkdata
	syntax [anything], [n(integer 100) sigma(real .5)]
	
	clear
	set obs `n'	
	generate d        = runiform()>.9
	generate e        = rnormal(0, `sigma') 
	quietly replace e = rnormal() if d==1
	generate y = 1 + d + e 
end

postfile sims rep vce b1 se_b1 r_b1 bias1 hatmax hatmin N rc	///
	      using harmless_1000_sim, replace 
	      
forvalues i=1/`R' {	
	quietly mkdata, n(`N')	
	// Default regresssion and hat matrix 
	capture regress y i.d
	quietly predict double hat if e(sample), hat
	local rank =  e(rank)
	local rc = _rc 
	if (!`rc' & `rank'==2) {
		local NS         = e(N)
		local b1         =  _b[1.d]
		local se1        =  _se[1.d]
		qui test _b[1.d] == `tb1'
		local b1r        = (r(p)<.05) 
		local bs1        = `b1' - `tb1'
		summarize hat, meanonly 
		local hmx        = r(max)
		local hmn        = r(min)
	}
	else {
		local NS   = .
		local b1   = .
		local se1  = .
		local b1r  = .
		local bs1  = .
		local hmx  = .
		local hmn  = .		
	}
	post sims (`i') (1) (`b1') (`se1') (`b1r') (`bs1') (`hmx')	///
	          (`hmn') (`NS') (`rc')	
		  
	capture regress y i.d, vce(robust)
	local rc = _rc 
	if (!`rc' & `rank'==2) {
		local NS         = e(N)
		local b1         =  _b[1.d]
		local se1        =  _se[1.d]
		qui test _b[1.d] == `tb1'
		local b1r        = (r(p)<.05) 
		local bs1        = `b1' - `tb1'
		summarize hat, meanonly 
		local hmx        = r(max)
		local hmn        = r(min)
	}
	else {
		local NS   = .
		local b1   = .
		local se1  = .
		local b1r  = .
		local bs1  = .
		local hmx  = .
		local hmn  = .		
	}
	post sims (`i') (2) (`b1') (`se1') (`b1r') (`bs1') (`hmx')	///
	          (`hmn') (`NS') (`rc')	
		  
	capture regress y i.d, vce(hc2)
	local rc = _rc 
	
	if (!`rc' & `rank'==2) {
		local NS         = e(N)
		local b1         =  _b[1.d]
		local se1        =  _se[1.d]
		qui test _b[1.d] == `tb1'
		local b1r        = (r(p)<.05) 
		local bs1        = `b1' - `tb1'
		summarize hat, meanonly 
		local hmx        = r(max)
		local hmn        = r(min)
	}
	else {
		local NS   = .
		local b1   = .
		local se1  = .
		local b1r  = .
		local bs1  = .
		local hmx  = .
		local hmn  = .		
	}
	post sims (`i') (3) (`b1') (`se1') (`b1r') (`bs1') (`hmx')	///
	          (`hmn') (`NS') (`rc')	
		  
	capture regress y i.d, vce(hc3)
	local rc = _rc 
	if (!`rc' & `rank'==2) {
		local NS         = e(N)
		local b1         =  _b[1.d]
		local se1        =  _se[1.d]
		qui test _b[1.d] == `tb1'
		local b1r        = (r(p)<.05) 
		local bs1        = `b1' - `tb1'
		summarize hat, meanonly 
		local hmx        = r(max)
		local hmn        = r(min)
	}
	else {
		local NS   = .
		local b1   = .
		local se1  = .
		local b1r  = .
		local bs1  = .
		local hmx  = .
		local hmn  = .		
	}
	post sims (`i') (4) (`b1') (`se1') (`b1r') (`bs1') (`hmx')	///
	          (`hmn') (`NS') (`rc')	
		  
	local rc = _rc 
	if (!`rc' & `rank'==2) {
		local NS         = e(N)
		local b1         = _b[1.d]
		local se1        = _se[1.d]
		quietly boottest 1.d, nograph 
		matrix CI        = r(CI)
		local lo         = CI[1,1]
		local up         = CI[1,2]
		local b1r        = 1 - (`tb1'<=`up')*(`tb1'>=`lo')  
		local bs1        = `b1' - `tb1'
	}
	else {
		local NS   = .
		local b1   = .
		local se1  = .
		local b1r  = .
		local bs1  = .
		local hmx  = .
		local hmn  = .		
	}
	post sims (`i') (5) (`b1') (`se1') (`b1r') (`bs1') (`hmx')	///
	          (`hmn') (`NS') (`rc')	
		  
	if (`i'/50) == int(`i'/50) {
		di ".                 `i'"
	}
	else {
		di _c "."
	}
}

postclose sims 

use harmless_1000_sim, clear 

label define estvce 1 "OLS" 2 "HC1" 3 "HC2" 4 "HC3" 5 "WILDBOOT"
label values vce estvce