R/GRS.test.R

In PengWu12245/EcoPleio: Testing Factor Models for Asset Returns: Based on Pleiotropy Model

#' A Test of the Efficiency of a Given Portfolio (Gibbons, Ross, Shanken, 1989)
#'
#' The test is for the following multivariate linear regression model:
#' \deqn{y_j = \alpha_j + \sum_{k=1}^q \beta_{jk}x_k + \varepsilon_j,   j=1,...,p.}
#' where \eqn{y_j} denotes the excess return on asset \eqn{j};
#' \eqn{(x_1,...,x_q)} is the excess return on the porfolio whose
#' efficiency is being tested; and \eqn{\varepsilon_j} is the disturbance
#' term for asset \eqn{j}. The disturbances are assumed to be jointly
#' normally distributed with mean zero and nonsingular covariance matrix
#' \eqn{\Sigma}, conditional on the excess returns for portfolios \eqn{(x_1,...,x_q)}.
#'
#' The test of the efficiency of a given portfolio is equivalent to
#' the following hypothesis test.
#' \itemize{
#'     \item \code{H0}: all the Intercepts \eqn{\alpha} are zero, i.e.
#'                      \eqn{\alpha_j=0,   \forall j=1,...,p.}
#'     \item \code{H1}: otherwise.
#' }
#'@param x samples of predictor which is a \eqn{n*q} matrix.
#'@param y samples of response which is a \eqn{n*p} vector.
#'@return A list with the following elements:
#' \itemize{
#'    \item \code{p.value}: the \eqn{p}-value of the GRS test.
#'    \item \code{grs.stat}: the test statistic of the GRS test.
#' }
#'@examples
#' ## Quick example for the GRS test
#'
#'
#' set.seed(1)
#' n = 200; q = 2; p = 3
#' x <- matrix( rnorm(n*q), nrow = n)
#'
#' # Generate data under H0
#' y <- matrix(NA, nrow = n, ncol = p)
#' y[,1] <- x[,1] + x[,2] + rnorm(n, sd = 0.5)
#' y[,2] <- x[,1] + 2 * x[,2] + rnorm(n, sd = 0.5)
#' y[,3] <- x[,1] + 3 * x[,2] + rnorm(n, sd = 0.5)
#' GRS.test(x, y)$p.value
#'
#' # Generate data under H1
#' y <- matrix(NA, nrow = n, ncol = p)
#' y[,1] <- 1 + x[,1] + x[,2] + rnorm(n, sd = 0.5)
#' y[,2] <- 1 + x[,1] + 2 * x[,2] + rnorm(n, sd = 0.5)
#' y[,3] <- x[,1] + 3 * x[,2] + rnorm(n, sd = 0.5)
#' GRS.test(x, y)$p.value
#'
#' @references
#' Gibbons, M. R., Ross, S. A. and Shanken, J. (1989).
#' A test of the efficiency of a given portfolio. Econometrica, 57(5), 1121-1152.

GRS.test <- function(x, y){

  # Initial Values
  p <- ncol(y)
  y <- unname(as.matrix(y))

  # alpha = 0.05   # the significant level of given.
  Y <- y
  mean.y <- matrix(colMeans(y), ncol = 1)
  dev.y <- mean.y - mean(mean.y)

  tmp.x <- x          # factors involved
  mean.x <- matrix(colMeans(tmp.x), ncol = 1)   # sample means
  omega.x <- var(tmp.x)   # sample variance-covariance matrix
  n <- nrow(tmp.x)    # the number of time points
  q <- ncol(tmp.x)    # q assets
  # X <- .Internal(cbind(1, 1, tmp.x)) # X <- cbind(rep(1, n), x
  X <- cbind(1, tmp.x)
  coeff <- solve(t(X) %*% X) %*% (t(X) %*% Y)
  intercept <- matrix(coeff[1,], ncol = 1)    # intercept
  model <- lm(Y~X[,2:(q+1)])   # linear model
  sigma <- t(model$residuals) %*% model$residuals / (n-q-1)

  # GRS test statistics ~ F(p, n-p-q)
  grs.stat <- n*(n-p-q) / (p*(n-q-1)) * (t(intercept) %*% solve(sigma) %*% intercept) /
    (1 + (t(mean.x) %*% solve(omega.x) %*% mean.x))
  # criterion <- qf(1-alpha, p, n-p-q)
  p.value <- 1 - pf(grs.stat, p, n-p-q)
  return(list(p.value = p.value, grs.stat = grs.stat))
  # return(p.value)
  # return(grs.stat)
  # return(mean(abs(intercept)))
  # return(mean(abs(intercept)) / mean(abs(dev.y)))
  # return(mean((intercept)^2) / mean((dev.y)^2))
}