nulldata 100

setobs 1 1900 --time-series

scalar iters = 1000

matrix B = zeros(iters, 4)
matrix SE = zeros(iters, 4)

genr x = uniform()

scalar r = 0.8
scalar b0 = 10
scalar b1 = 20
scalar b2 = -r*b1*2

set messages off

loop i=1..iters -q
  genr e = normal()
  genr y = b0 + b1*x
  if 1
    genr y = (1-r)*b0 + r*y(-1) + b1*x - r*b1*x(-1) + e
  else
    genr y = (1-r)*b0 + r*y(-1) + b1*x + b2*x(-1) + e
  endif

  smpl 1970 ;	

  ols y 0 x -q
  B[i,1] = $coeff(x)
  SE[i,1] = $stderr(x)

  ols y 0 x -qr
  B[i,2] = $coeff(x)
  SE[i,2] = $stderr(x)

  corc y 0 x -q
  B[i,3] = $coeff(x)
  SE[i,3] = $stderr(x)

  ols y 0 y(-1) x x(-1) -q
  B[i,4] = $coeff(x)
  SE[i,4] = $stderr(x)

  smpl --full
endloop 

matrix Bmean = meanc(B)
matrix Bsd = sd(B)
matrix SEmean = meanc(SE)

printf "\n"
printf "OLS1: mean(b1) = %.4f, sd(b1) = %.4f, avg(se)/sd(b1) = %.3f\n", \
  Bmean[1], Bsd[1], SEmean[1]/Bsd[1]
printf "HAC:  mean(b1) = %.4f, sd(b1) = %.4f, avg(se)/sd(b1) = %.3f\n", \
  Bmean[2], Bsd[2], SEmean[2]/Bsd[2]
printf "CORC: mean(b1) = %.4f, sd(b1) = %.4f, avg(se)/sd(b1) = %.3f\n", \
  Bmean[3], Bsd[3], SEmean[3]/Bsd[3]
printf "OLS2: mean(b1) = %.4f, sd(b1) = %.4f, avg(se)/sd(b1) = %.3f\n", \
  Bmean[4], Bsd[4], SEmean[4]/Bsd[4]