capture log close
log using VtFsim, text replace

clear all

local n=1000                      //sample size
local R=2000                      //number of replications

local beta=1
local f0=3
local pi=(`f0')/sqrt(`n')         //Argument of MC()

//Taken from Mostly Harmless Econometrics
capture program drop MC
program define MC, rclass
{
  syntax [, n(integer 1000) rho(real 0.8) pi(real 0.04472136) beta(real 1)]
  drop _all                                     //NOTE: This drops the data, but not the locals/globals (pi,rho,n)
  
  set obs `n'
  gen v=rnormal()                               //v, Z, and epsilon are all N(0,1)
  gen Z=rnormal()
  gen epsilon=rnormal()                         //epsilon is used to construct u
  gen u = `rho'*v + sqrt(1-(`rho')^2)*epsilon   //Var(u)=1
  gen X = `pi'*Z + v 
  gen Y = `beta'*X + u
  keep X Y Z
  vtf Y (X=Z), robust beta(`beta')              //run IV, test true null
  local bhat=e(beta_hat)                        //IV coefficient itself
  local se_bhat=e(unadj_se)                     //"regular" IV standard error (Wald)
  local F=e(F)                                  //capture F statistic
  local rhat = e(r_hat)                         //extra information needed: rhat
  local VtFLB=e(VtF_LB_05)                      //if you want to do 1% test, use the -level- option of -vtf-
  local VtFUB=e(VtF_UB_05)
  local VtFcritval=e(VtF_critical_value_05)
  
  //What follows is Stata gobbledy-gook for getting the locals passed to the -simulate- command
  return scalar bhat=`bhat'
  return scalar se_bhat=`se_bhat'
  return scalar F=`F'
  return scalar rhat=`rhat'
  return scalar VtFLB=`VtFLB'
  return scalar VtFUB=`VtFUB'
  return scalar VtFcritval=`VtFcritval'
  ereturn clear
}
end

//**************************************
//*** Run the simulation ***************
//**************************************
set seed 123456789
//set rmsg on
simulate, reps(`R'): MC, n(`n') pi(`pi') beta(`beta')

gen byte cover=0
replace cover=1 if VtFLB==.|VtFUB==.
replace cover=1 if VtFLB!=. & VtFLB<1 & VtFUB>1 & VtFUB!=.
gen byte reject=0
replace reject=1 if abs((bhat-(`beta'))/se_bhat)>VtFcritval & VtFcritval!=.

summ cover reject

gen tmp=0
replace tmp=1 if VtFLB==. & VtFUB!=.
replace tmp=1 if VtFUB==. & VtFLB!=.
assert tmp==0

log close