R/PleioSeq.test.H00.R

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

#' A sequential test of the Efficiency of a Given Portfolio under \code{H00}.
#'
#' 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{H00}: 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{pvalue}: the \eqn{p}-value of the sequential test for \code{H00}.
#'    \item \code{stat}: the test statistic for \code{H00}.
#' }
#'@examples
#' ## Quick example for the sequential test for \code{H00}
#'
#'
#' set.seed(1)
#' n = 500; q = 2; p = 3
#' x <- matrix( rnorm(n*q), nrow = n)
#'
#' # 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
#' PleioSeq.test.H00(x,y)$pvalue
#'
#' # Under H1
#' y <- matrix(NA, nrow = n, ncol = p)
#' y[,1] <- 0.5 + x[,1] + x[,2] + rnorm(n, sd = 0.5)
#' y[,2] <- 0.5 + 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
#' PleioSeq.test.H00(x,y)$pvalue


PleioSeq.test.H00 <- function(x,y){

  p <- ncol(y)
  y <- unname(as.matrix(y))

  alpha = 0.05   # the significant level of given.
  Y <- y

  tmp.x <- x
  n <- nrow(tmp.x)
  q <- ncol(tmp.x)
  # X <- .Internal(cbind(1, 1, tmp.x)) # X <- cbind(rep(1, n), x
  X <- cbind(1, tmp.x)
  
  # step 1. Compute the parameter matrix and An
  Pi <- .lm.fit(X,Y)$coefficients   # lmC(X, Y)
  Omega <- OmegaC(X, Y, Pi, n)
  An.hat <- An.hat.est(Omega, tmp.x, y, n, p, q)
  An_inv = solve(An.hat)

  beta.hat <- matrix(c(Pi[1, ],  as.vector( Pi[-1, ]) ), ncol=1)
  # V0 <- .Internal(diag(1, p*(q+1), p*(q+1)))[(1:p), ]   # intercept test
  V0 <- diag(1, p*(q+1), p*(q+1))[(1:p), ] 
  
  # step 2. Compute the stat. and pvalue. for H00
  # k.test = 0
  V.tmp <- V0
  T.tmp <-  T0C(V.tmp, An_inv, beta.hat)
  pval <- 1 - pchisq(T.tmp, df = p-0 )

  return( list( pvalue = pval,
                stat = T.tmp) )

}