R/PleioSeq.test.R

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

#' Testing Factor Models for Asset Returns: Based on Pleiotropy Model
#'
#' The factor models for asset return through a linear function are
#' commonly used in the description of asset returns in finance.
#' Also consider 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)}.
#'
#' Sequential testing provides a rigorous procedure to evaluate the
#' number of economic quantities associated with a fundamental economic factor.
#' Thus we consider a sequential hypothesis test:
#' \itemize{
#'     \item \code{H0k}: there is no more than \code{k} nonzero \eqn{\alpha_j (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.
#'    \item \code{k.test}: the number of nonzero \eqn{\alpha_j (j=1,...,p)}.
#'    \item \code{index}: the corresponding indexs of nonzero \eqn{\alpha_j (j=1,...,p)}.
#'}
#'@examples
#' ## Quick example for the sequential test of the efficiency of a given portfolio
#'
#' set.seed(1)
#' n = 500; 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)
#' PleioSeq.test(x,y)
#'
#' # Generate data under H01
#' y <- matrix(NA, nrow = n, ncol = p)
#' y[,1] <- 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)
#' PleioSeq.test(x,y)
#'
#' # Generate data under H02
#' y <- matrix(NA, nrow = n, ncol = 3)
#' 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)
#' PleioSeq.test(x,y)

PleioSeq.test <- 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. for a given combinations
  # k.test = 0
  V.tmp <- V0
  T.tmp <-  T0C(V.tmp, An_inv, beta.hat)
  pval <- 1 - pchisq(T.tmp, df = p-0 )
  if(pval > alpha){
    return( list( index = 0, pvalue = pval, k.test = 0) )
  }

  for(k.test in 1:p){
    result <- combn(p, k.test)
    stat <-  apply(result, 2, function(com){
      V.tmp <- V0[-com, , drop = F]
      T.tmp <-  T0C(V.tmp, An_inv, beta.hat)
      return(T.tmp)
    })

    pval <- 1 - pchisq(min(stat), df = p-k.test )  # intercept test
    index = result[, which.min(stat)]

    if(pval > alpha){  break }
  }


  return( list( index = index, pvalue = pval,
                k.test = k.test) )

}