>  # glsexample.r
>  #
>  # simple example of generalized least squares
>  #
>   sigmasq = 1                 # KISS
>   N <- 6                      # number of observations
>  # design matrix for simple linear regression
>   X <- matrix( c(rep(1,N),(1:N)),N,2)
>   X                           # write it out
     [,1] [,2]
[1,]    1    1
[2,]    1    2
[3,]    1    3
[4,]    1    4
[5,]    1    5
[6,]    1    6
>  # cov matrix for AR(1) process with rho = 0.4
>   rho <- .4
>   V <- rho^abs( matrix(1:N,N,N) - matrix(1:N,N,N,byrow=T) )
>   V                           # write it out
        [,1]   [,2]  [,3]  [,4]   [,5]    [,6]
[1,] 1.00000 0.4000 0.160 0.064 0.0256 0.01024
[2,] 0.40000 1.0000 0.400 0.160 0.0640 0.02560
[3,] 0.16000 0.4000 1.000 0.400 0.1600 0.06400
[4,] 0.06400 0.1600 0.400 1.000 0.4000 0.16000
[5,] 0.02560 0.0640 0.160 0.400 1.0000 0.40000
[6,] 0.01024 0.0256 0.064 0.160 0.4000 1.00000
>  # cov matrix for gls is sigmasq * inv(X'inv(V)X)
>   covgls = sigmasq*solve( t(X) %*% solve(V) %*% X )
>   covgls                      # write it out
          [,1]        [,2]
[1,]  1.228801 -0.26017699
[2,] -0.260177  0.07433628
>  # cov matrix for ols is sigmasq * inv(X'X) X'V X inv(X'X)
>   covols = sigmasq*solve(t(X)%*%X) %*% t(X) %*% V %*% X %*% solve(t(X)%*%X)
>   covols                      # write it out
           [,1]        [,2]
[1,]  1.2754133 -0.27085714
[2,] -0.2708571  0.07738776
>  # difference
>   covols - covgls
            [,1]         [,2]
[1,]  0.04661205 -0.010680152
[2,] -0.01068015  0.003051472
>  # is the difference positive definite?
>   eigen(covols - covgls)$values
[1] 0.0490896762 0.0005738418
>  # done
>  rm(list=ls())                # clean up
>  q()                          # quit R