# 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
 # 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
 # cov matrix for gls is sigmasq * inv(X'inv(V)X)
  covgls = sigmasq*solve( t(X) %*% solve(V) %*% X )
  covgls                      # write it out
 # 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
 # difference
  covols - covgls
 # is the difference positive definite?
  eigen(covols - covgls)$values
 # done
 rm(list=ls())                # clean up
 q()                          # quit R