R/FastMultiHorizonMCS.R

In lucabarbaglia/MultiHorizonSPA: Multi Horizon Superior Predictive Ability

#' @title Fast Multihorizon model confidence set and p-values
#'
#' @description Produces Multi-Horizon Model Confidence Set and p-values as described in section 2.2. of Quadvlieg (2021).
#' @param Losses A list containing M matrices of dimension T by H. Each matrix contains the losses for a forecasting method. Rows correspond to time periods (T rows), columns correspond to forecast horizons (H columns).
#' @param alpha_t Alpha level of critical values for pairwise tests. See Quadvlieg (2021).
#' @param alpha_mcs Alpha level of the critical value for the MCS test. Denoted by alpha tilde in Quadvlieg (2021) section 2.2
#' @param weights the 1 x H vector of weights for the losses at different horizons. For instance \code{weights <- matlab::ones(1,20)/20}
#' @param L integer, the parameter for the moving block bootstrap
#' @param B integer, the number of bootstrap iterations. Default 999
#' @param unif_or_average If equals "a", then the confidence set is calculated based on average SPA. If equals "u", then the confidence set is calculated based on uniform SPA.
#' @param num_cores integer, the number of processor cores to use in parallel.
#' @param seed integer, the seed for random number generation
#' @export
#' @return Returns a list of length 2. The elements of the list are:
#' \item{MCS_set}{A M by 2 matrix. The first column contains the indices of the methods for which the p-value is greater than or equal to 1-alpha_mcs. The second column contains the p-values that are greater than or equal to 1-alpha_mcs.}
#' \item{p_values}{A vector containing the p-values for all methods.}
#' @examples
#' 
#' 
#' Trow <- 20 
#' H <- 12
#' Mmethods <- 9
#' weights <- rep(1/H,H)
#' 
#' loss_list <- vector(mode = "list", length = Mmethods)
#' 
#' loss_list[[1]] <- matrix(rnorm(Trow*H, mean = 1), nrow = Trow, ncol = H)
#' loss_list[[2]] <- matrix(rnorm(Trow*H, mean = 2), nrow = Trow, ncol = H)
#' loss_list[[3]] <- matrix(rnorm(Trow*H, mean = 3), nrow = Trow, ncol = H)
#' loss_list[[4]] <- matrix(rnorm(Trow*H, mean = 2), nrow = Trow, ncol = H)
#' loss_list[[5]] <- matrix(rnorm(Trow*H, mean = 1), nrow = Trow, ncol = H)
#' loss_list[[6]] <- matrix(rnorm(Trow*H, mean = 1), nrow = Trow, ncol = H)
#' loss_list[[7]] <- matrix(rnorm(Trow*H, mean = 2), nrow = Trow, ncol = H)
#' loss_list[[8]] <- matrix(rnorm(Trow*H, mean = 3), nrow = Trow, ncol = H)
#' loss_list[[9]] <- matrix(rnorm(Trow*H, mean = 2), nrow = Trow, ncol = H)
#' loss_list[[10]] <- matrix(rnorm(Trow*H, mean = 1), nrow = Trow, ncol = H)
#' 
#' 
#' num_cores <- 1
#' 
#' seed <- 42
#' FastMultiHorizonMCS(loss_list, #
#'                                       0.05, # alpha_t
#'                                       0.05, # alpha_mcs
#'                                       weights, #
#'                                       3,#l
#'                                       5,#b
#'                                       "u", # 1 means uniform
#'                                       num_cores,
#'                                       seed)
#' 
#' 
#' @export
FastMultiHorizonMCS <- function(Losses,
                                alpha_t = 0.05,
                                alpha_mcs = 0.05,
                                weights = NULL, 
                                L, 
                                B=999,
                                unif_or_average = "u",
                                num_cores = 1,
                                seed = stats::runif(1, 0, .Machine$integer.max)){
  
  
  if(!(unif_or_average %in% c("a","u"))){
    stop("unif_or_average must equal 'u' for uniform or 'a' for average.")
  }
  
  if ((!is.matrix(Losses[[1]])) & (!is.data.frame(Losses[[1]]))  ){
    stop("Elements of Losses list must be matrices or data frames.")
  }
  
  #note: this 
  if ((!is.list(Losses)  )   ){
    stop("Losses must be a list")
  }
  if(is.data.frame(Losses)){
    stop("Losses cannot be a data frame")
  }
  
  unif_avg  <- 1
  if(unif_or_average =="a"){
	  unif_avg <- 0
  }
  
  
  # Set <- Losses
  #p_values <- rep(0, length(Losses))
  #IDs <- 1:length(Losses)
  
  #Trows <- nrow(Losses[[1]])
  Hcols <- ncol(Losses[[1]])
  
  Mmethods <- length(Losses)
  
  if(is.null(weights)){
    weights <- rep(1/Hcols, Hcols)
  }
  if(unif_or_average =="a" &  length(weights)!=Hcols){
	print("length of weights not equal to number of columns (horizons), setting equal weights")
    weights <- rep(1/Hcols, Hcols)
  }
  
  #in case list contains data frames, convert to matrices
	loss_list <- vector(mode = "list", length = Mmethods)

	
	for(i in 1:Mmethods){
	  loss_list[[i]] <- as.matrix(Losses[[i]])
	}

  
  
  
  ret <- MultiHorizonMCS_cpp(loss_list, 
                         alpha_t,
                         alpha_mcs,
                         weights,
                         L, 
                         B,
                         unif_avg,
                         num_cores,
                         seed
                           )
  
  
  names(ret) <- c("MCS_set", "p_values")
  
  return(ret)
  
  }