GLS_chol {varycoef}

GLS Estimate using Cholesky Factor

Description

Computes the GLS estimate using the formula:

\mu_{GLS} = (X^\top \Sigma^{-1} X)^{-1}X^\top \Sigma^{-1} y.

The computation is done depending on the input class of the Cholesky factor R. It relies on the classical solve or on using forwardsolve and backsolve functions of package spam, see solve. This is much faster than computing the inverse of \Sigma, especially since we have to compute the Cholesky decomposition of \Sigma either way.

Usage

GLS_chol(R, X, y)

## S3 method for class 'spam.chol.NgPeyton'
GLS_chol(R, X, y)

## S3 method for class 'matrix'
GLS_chol(R, X, y)

Arguments

R

(spam.chol.NgPeyton or matrix(n, n))
Cholesky factor of the covariance matrix \Sigma. If covariance tapering and sparse matrices are used, then the input is of class spam.chol.NgPeyton. Otherwise, R is the output of a standard chol, i.e., a simple matrix

X

(matrix(n, p))
Data / design matrix.

y

(numeric(n))
Response vector

Value

A numeric(p) vector, i.e., the mean effects.

Author(s)

Jakob Dambon

Examples

# generate data
n <- 10
X <- cbind(1, 20+1:n)
y <- rnorm(n)
A <- matrix(runif(n^2)*2-1, ncol=n)
Sigma <- t(A) %*% A
# two possibilities
## using standard Cholesky decomposition
R_mat <- chol(Sigma); str(R_mat)
mu_mat <- GLS_chol(R_mat, X, y)
## using spam
R_spam <- chol(spam::as.spam(Sigma)); str(R_spam)
mu_spam <- GLS_chol(R_spam, X, y)
# should be identical to the following
mu <- solve(crossprod(X, solve(Sigma, X))) %*%
      crossprod(X, solve(Sigma, y))
## check
abs(mu - mu_mat)
abs(mu - mu_spam)