GLS_chol function

GLS Estimate using Cholesky Factor

Computes the GLS estimate using the formula: [REMOVE_ME]$\mu_{GLS} = (X^\top \Sigma^{-1} X)^{-1}X^\top \Sigma^{-1} y. [REMOVE_ME_2]$

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.

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

Returns

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

Description

Computes the GLS estimate using the formula:

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.

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)

Author(s)

Jakob Dambon