*
* Monte Carlo Analysis of OLS vs GLS, page 165
*
all 100
*
compute ndraws=20000
*
* The stderrs matrix will keep the intercept standard errors in the first row and
* the slope standard errors in the second. The columns are OLS, OLS corrected and
* GLS respectively. These are accumulated across the draws, then converted to
* mean values at the end by dividing by the number of draws.
*
dec rect stderrs(2,3)
declare real alpha
dofor alpha = .5 1.0 2.0 3.0
*
*     The OLS standard errors are computed using the unconditional value for the
*     variance, rather than the sample estimate. Since x is U(0,1), the expected
*     value of x**alpha is 1/(1+alpha).
*
*     The uncorrected OLS standard errors are computed by taking the square roots of
*     the diagonal elements of sigma**2*inv(X'X). The corrected OLS standard errors
*     are the diagonal elements of inv(X'X) X'VX inv(X'X). The center matrix can be
*     computed by using the CMOMENT instruction with a spread equal to the
*     reciprocal of the x**alpha series.
*
   compute ucvar=1.0/(1+alpha)
   compute stderrs=%const(0.0)
   do draw=1,ndraws
      set x   = %uniform(0,1)
      set xs  = x**alpha
      set xss = 1./xs
      set y   = 1+x+%ran(sqrt(x))
      linreg(noprint) y
      # constant x
      compute stderrs(1,1)=stderrs(1,1)+sqrt(ucvar*%xx(1,1)),stderrs(2,1)=stderrs(2,1)+sqrt(ucvar*%xx(2,2))
      cmom(spread=xss)
      # constant x
      compute %xx=%xx*%cmom*%xx
      compute stderrs(1,2)=stderrs(1,2)+sqrt(%xx(1,1)),stderrs(2,2)=stderrs(2,2)+sqrt(%xx(2,2))
      linreg(noprint,spread=xs) y
      # constant x
      compute stderrs(1,3)=stderrs(1,3)+sqrt(%xx(1,1)),stderrs(2,3)=stderrs(2,3)+sqrt(%xx(2,2))
   end do draw
   compute stderrs=stderrs*(1.0/ndraws)
   disp "Std Errors for alpha=" alpha
   disp stderrs
   disp
end dofor alpha