# Set the sample size and save it in n
set hac_lag nw1
nulldata 400
scalar n = $nobs
setobs 1 1 --time-series

# generate n observations on x   
series x = uniform()

set seed 3213799 # set a seed if you want to get same results

# Set the values of the parameters
scalar slope = 1
scalar sigma = 1
scalar rho = .880
scalar alpha = .05

# initialize variables u
series u = normal()

# start the loop, indicating the desired number of samples.
loop 400 --progressive --quiet
   # generate normal errors
   series e = normal(0,sigma)
   # generate autocorrelated errors and x
   series u=rho*u(-1)+e
   series x=rho*x(-1)+normal()
   # generate sample of y
   series y = x + u
   # Estimate the model using OLS, save the slope estimates
   ols y const x
   genr b2 = $coeff(x)
   genr s = $stderr(x)
   # generate the lower and upper bounds for the confidence interval
   genr c2L = b2 - critical(t,$df,alpha)*s
   genr c2R = b2 + critical(t,$df,alpha)*s
   # count the number of instances when coefficient is inside inverval
   genr p2 = (slope>c2L && slope<c2R)
   # Repeat for the HAC covariance estimator --robust 
   ols y const x --robust 
   genr sr = $stderr(x)
   # generate the lower and upper bounds for the confidence interval
   genr c2Lr = b2 - critical(t,$df,alpha)*sr
   genr c2Rr = b2 + critical(t,$df,alpha)*sr
   # count the number of instances when coefficient is inside inverval
   genr p2r = (slope>c2Lr && slope<c2Rr)
   print p2 p2r
endloop