R/testKron_auto.R

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

#' Testing Kronecker product structure without specified modes
#'
#' @description Testing the Kronecker product structure of a tensor time series without 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 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.
#' @param all Logical. when TRUE, all sets of mode indices are tested; otherwise each reshaped data using leave-one-out sets is sequentially tested. Default is FALSE.
#'
#'
#' @return A list containing the following:
#' decision: a data frame with each row reporting the decision rule aggregated by quantile of parameter 'q' for different levels of alpha; the first column represents all non-trivial sets of mode indices to test if 'all' is TRUE, otherwise represents each mode to identify.
#' level: a data frame with each row reporting an example test size aggregated by quantile of parameter 'q' for different levels of alpha; the first column represents all non-trivial sets of mode indices to test if 'all' is TRUE, otherwise represents each mode to identify.
#' rank: a vector of integers representing either the input rank or the estimated rank used in testing.
#'
#'
#' @export
#' @import tensorMiss
#' @import utils
#'
#'
#'
#'
testKron_auto <- function(ten, r=0, alpha=c(0.01, 0.05), q=0.05, r.exact=FALSE, all=FALSE){
  if (requireNamespace(c('tensorMiss'), 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()
    }

    if ( r[1]!=0 ){
      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()
      }
    }

    ten_order <- length(dim(ten))-1

    ### estimate rank for r[1] being 0;
    if ( r[1]==0 ){
      r <- numeric(ten_order)
      for (j in 2:length(dim(ten))){
        ej_val <- svd( tensorMiss::unfold(ten, j) %*% t(tensorMiss::unfold(ten, j)) )$d[1:ceiling(dim(ten)[j]/2)]
        r[j-1] <- which.min( ((ej_val[-1])/(ej_val[-length(ej_val)])) )
      }
    }

    ### create a list containing sets of mode indices
    if (ten_order==2){
      AA_list <- list(c(1,2))
    } else {
      if (all){
        AA_list <- do.call("c", lapply(2:ten_order, function(i) utils::combn(1:ten_order, i, FUN = list)))
      } else {
        ## list containing the mode to exclude each time
        AA_list <- lapply(1:ten_order, function(i) (1:ten_order)[i])
        # AA_list <- lapply(1:ten_order, function(i) (1:ten_order)[-i])
      }
    }

    ### testing for each set of mode indices
    test_A_mat <- matrix(nrow=length(AA_list), ncol=length(alpha)+1)
    test_A_mat_decision <- matrix(nrow=length(AA_list), ncol=length(alpha)+1)
    for (j in 1:length(AA_list)){
      mode_name <- ""
      for (i in 1:length(AA_list[[j]])){ mode_name <- paste0(mode_name, ",", AA_list[[j]][i]) }
      if ( (!all) & (ten_order!=2) ){
        test_A_mat[j,1] <- substr(mode_name, 2, nchar(mode_name))
        test_A_mat_decision[j,1] <- substr(mode_name, 2, nchar(mode_name))
      } else {
        test_A_mat[j,1] <- paste0("{", substr(mode_name, 2, nchar(mode_name)), "}")
        test_A_mat_decision[j,1] <- paste0("{", substr(mode_name, 2, nchar(mode_name)), "}")
      }

      if ( (!all) & (ten_order!=2) ){
        testKron_j <- testKron_A(ten_reshape(ten, seq_len(length(AA_list))[-AA_list[[j]]] ), c(1,2), c(r[AA_list[[j]]], prod(r[-AA_list[[j]]])), alpha, q=q, r.exact)
      } else {
        testKron_j <- testKron_A(ten, AA_list[[j]], r, alpha, q, r.exact)
      }

      test_A_mat[j,2:ncol(test_A_mat)] <- apply(testKron_j$decision, 2, min)
      test_A_mat_decision[j,2:ncol(test_A_mat_decision)] <- as.numeric(test_A_mat[j,2:ncol(test_A_mat)] > alpha)
    }

    test_A_mat <- as.data.frame(test_A_mat)
    test_A_mat_decision <- as.data.frame(test_A_mat_decision)
    if ( (!all) & (ten_order!=2) ){
      colnames(test_A_mat) <- c("Mode to be identified", as.character(alpha))
      colnames(test_A_mat_decision) <- c("Mode to be identified", as.character(alpha))
    } else {
      colnames(test_A_mat) <- c("Set of Modes to be tested", as.character(alpha))
      colnames(test_A_mat_decision) <- c("Set of Modes to be tested", as.character(alpha))
    }

    return( list(decision=test_A_mat_decision, level=test_A_mat, rank=r) )

  }
}