R/testKron_A.R

In KOFM: Test the Kronecker Product Structure in Tensor Factor Models

#' Testing Kronecker product structure along specified modes
#'
#' @description Testing the Kronecker product structure of a tensor time series with a specified set of mode indices.
#' @param ten An array representing an order-(K+1) tensor which is an order-K tensor time series with mode-1 being the time mode.
#' @param AA A vector with elements at least 1 and at most K, representing the tensor modes to test.
#' @param r A vector representing the rank for 'ten'.
#' @param alpha A vector representing the desired significance levels. Default is c(0.01, 0.05).
#' @param q A number between 0 and 1, representing the quantile for the decision statistic. Default is 0.05.
#' @param r.exact Logical. Perform the test only using 'r' if TRUE, otherwise using all divisor combinations of the last mode rank of the reshaped data. Default is FALSE.
#'
#'
#' @return A list containing the following:
#' level: a matrix with each entry reporting an example test size, corresponding to its rank used (row) and level of alpha (column); the alpha is reported in order of parameter 'alpha'.
#' decision: a matrix with each entry reporting the decision statistic aggregated by quantile of parameter 'q', corresponding to its rank used (row) and level of alpha (column); the alpha is reported in order of parameter 'alpha'.
#' rank: a matrix with K columns and each row represents a different rank used, corresponding to the rows in the 'level' and 'decision' matrices.
#'
#'
#' @export
#' @import tensorMiss
#' @import MEFM
#' @import RSpectra
#' @import stats
#'
#'
#'
#'
testKron_A <- function(ten, AA, r, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE){
  if (requireNamespace(c('tensorMiss', 'RSpectra'), quietly = TRUE)){

    if ( length(dim(ten)) <3 ){
      message("The length of dimension in `ten` should be at least 3 (i.e., matrix time series).")
      return()
    }

    AA <- sort(unique(AA))
    if ( (min(AA)<1)|(max(AA)> ( length(dim(ten)) -1)) ){
      message("The elements in `AA` can only be in between 1 and the order of the data.")
      return()
    }

    ten.reshape <- ten_reshape(ten, AA)

    if ( length(dim(ten)[-1]) != length(r) ){
      message("The length of `r` does not match the order of the data.")
      return()
    } else if ( sum(r > (dim(ten)[-1])) != 0 ){
      message("Some elements in `r` are larger than the dimension of the data.")
      return()
    }

    r.reshape <- c(r[-AA], prod(r[AA]))

    ### estimation: reshaped data
    loading.reshape <- list()
    for (j in 2:length(dim(ten.reshape))){
      loading.reshape[[j-1]] <- RSpectra::eigs_sym(tensorMiss::unfold(ten.reshape, j) %*% t(tensorMiss::unfold(ten.reshape, j)), r.reshape[j-1])$vectors
    }

    C.reshape <- tensorMiss::ttm(ten.reshape, (loading.reshape[[1]] %*% t(loading.reshape[[1]])), 2)
    if (length(loading.reshape)>1){
      for (j in 2:length(loading.reshape)){
        C.reshape <- tensorMiss::ttm(C.reshape, (loading.reshape[[j]] %*% t(loading.reshape[[j]])), j+1)
      }
    }

    e.reshape <- (ten.reshape - C.reshape)^2
    e.reshape <- ten_reshape_back(e.reshape , AA, dim(ten))

    d.min <- which.min(dim(ten)[-1])
    e.reshape.test <- tensorMiss::ttm(e.reshape, matrix(1, nrow=1, ncol=dim(ten)[-1][d.min]), d.min+1)/(dim(ten)[-1][d.min])
    e.reshape.test <- tensorMiss::unfold(e.reshape.test, 1)


    ### estimation: original data
    loading.original <- list()
    for (j in 2:length(dim(ten))){ loading.original[[j-1]] <- tensorMiss::unfold(ten, j) %*% t(tensorMiss::unfold(ten, j)) }

    ## matrix of divisor combination
    if (r.exact){
      search_mat <- matrix(r, nrow=1, ncol=length(r))
    } else {
      search_mat <- divisors(r.reshape[length(r.reshape)], dim(ten)[-1][AA])
      search_mat <- cbind(matrix(r[-AA], byrow = TRUE, nrow=nrow(search_mat), ncol=length(r[-AA])), search_mat)
    }

    test_vec.alpha <- matrix(nrow=nrow(search_mat), ncol=length(alpha))
    test_vec.quantile <- matrix(nrow=nrow(search_mat), ncol=length(alpha))
    for (j in seq_len(nrow(search_mat))){
      loading.j <- list()
      for (k in 1:length(loading.original)){
        loading.j[[k]] <- RSpectra::eigs_sym(loading.original[[k]], search_mat[j,k])$vectors
      }

      C.j <- tensorMiss::ttm(ten, (loading.j[[1]] %*% t(loading.j[[1]])), 2)
      for (k in 2:length(loading.j)){
        C.j <- tensorMiss::ttm(C.j, (loading.j[[k]] %*% t(loading.j[[k]])), k+1)
      }
      e.j <- (ten - C.j)^2
      e.j.test <- tensorMiss::ttm(e.j, matrix(1, nrow=1, ncol=dim(ten)[-1][d.min]), d.min+1)/(dim(ten)[-1][d.min])
      e.j.test <- tensorMiss::unfold(e.j.test, 1)

      test_vec.j <- matrix(nrow=ncol(e.j.test), ncol=length(alpha))
      for (i in 1:nrow(test_vec.j)){
        for (h in 1:ncol(test_vec.j)){
          test_vec.j[i,h] <- sum(e.j.test[,i] >= MEFM::qHat(e.reshape.test[,i], 1-alpha[h]))/nrow(e.j.test)
        }
      }

      test_vec.alpha[j,] <- test_vec.j[1,]
      test_vec.quantile[j,] <- apply(test_vec.j, 2, stats::quantile, probs=q)
    }

    return(list(level=test_vec.alpha, decision=test_vec.quantile, rank=search_mat))

  }
}