set verbose off

# Illustrate the effects of heteroskedasticity, using four different
# estimators including feasible and infeasible GLS

nulldata 100
K = 2000

matrix slope = zeros(K,4)
matrix se = zeros(K,4)

# set seed 1411994449

series x = uniform()
set hc_version 1

loop i=1..K -q
   # generate artificial data
   series sdu = 2 + 10*x
   series u = normal() * sdu
   series y = 10 + 20*x + u

   # regular OLS
   ols y 0 x --quiet
   slope[i,1] = $coeff[2]
   se[i,1] = $stderr[2]

   # OLS with robust standard errors
   ols y 0 x --quiet --robust
   slope[i,2] = $coeff[2]
   se[i,2] = $stderr[2]

   # automatic h'sked correction (FGLS)
   hsk y 0 x --quiet
   slope[i,3] = $coeff[2]
   se[i,3] = $stderr[2]

   # WLS using known weights (IGLS)
   series w = 1.0/(sdu*sdu)
   wls w y 0 x --quiet
   slope[i,4] = $coeff[2]
   se[i,4] = $stderr[2]
endloop

strings S = defarray("OLS", "HC1", "HSK", "WLS")

printf "\n"
loop j=1..4
   sd1 = sdc(slope[,j])
   sd2 = meanc(se[,j])
   printf "%s: mean(b) %.3f, s.d.(b) %.3f, avg(se) %.3f, diff %+.3f\n", \
     S[j], meanc(slope[,j]), sd1, sd2, sd1 - sd2
endloop