nulldata 100

scalar iters = 100

series b1ols = 0
series b1hc1 = 0
series b1hsk = 0
series b1wls = 0

series se1ols = 0
series se1hc1 = 0
series se1hsk = 0
series se1wls = 0

genr x = uniform()
set hc_version 1

loop i=1..iters -q
  genr sdu = 2 + 10*x
  genr u = normal() * sdu
  genr y = 10 + 20*x + u

  ols y 0 x -q
  b1ols[i] = $coeff(x)
  se1ols[i] = $stderr(x)

  ols y 0 x -q -r
  b1hc1[i] = $coeff(x)
  se1hc1[i] = $stderr(x)

  hsk y 0 x -q
  b1hsk[i] = $coeff(x)
  se1hsk[i] = $stderr(x)

  genr w = 1.0/(sdu*sdu)
  wls w y 0 x -q
  b1wls[i] = $coeff(x)
  se1wls[i] = $stderr(x)
endloop 

printf "\n"
printf "OLS: mean(b) = %.4f, s.d.(b) = %.4f, avg(se) = %.4f\n", \
  mean(b1ols), sd(b1ols), mean(se1ols)
printf "HC1: mean(b) = %.4f, s.d.(b) = %.4f, avg(se) = %.4f\n", \
  mean(b1hc1), sd(b1hc1), mean(se1hc1)
printf "HSK: mean(b) = %.4f, s.d.(b) = %.4f, avg(se) = %.4f\n", \
  mean(b1hsk), sd(b1hsk), mean(se1hsk)
printf "WLS: mean(b) = %.4f, s.d.(b) = %.4f, avg(se) = %.4f\n", \
  mean(b1wls), sd(b1wls), mean(se1wls)