R/rejection_rates.R

In asael697/nortstest: Assessing Normality of Stationary Process

#' Rejection rates estimates.
#'
#' Performs a rejection rate estimate for a stationary AR(1) process.
#'
#' @param reps an integer with the total bootstrap repetitions.
#' @param n an integer with the length of the simulated AR(1) process.
#' @param dist optional: a function to generate the innovations. By default rnorm.
#' @param seed An optional \code{\link[=set.seed]{seed}} to use.
#' @param htest A character string naming the desired test for checking normality. Valid values are
#' \code{"epps"} for the Epps, \code{"lobato"} for Lobato and Velasco's,\code{"vavra"} for the Psaradakis
#' and Vavra, and \code{"rp"} for the random projections. The default value is \code{"lobato"} test.
#' @param alpha Level of the test, possible values range from 0.01 to 0.1. By default \code{alpha = 0.05}.
#' @param k an integer with the number of random projections to be used, by default \code{k = 2}.
#' @param phi a real value with the auto regressive coefficient to simulate an stationary AR(1) process.
#' @param sd optional value for the standard deviation in a normal distribution.
#' @param ... additional arguments for rand.gen. Most usefully, the standard deviation of the innovations
#' generated by \code{rnorm} can be specified by \code{sd}.
#'
#' @return a real value with the proportional number of rejection for a defined test.
#'
#' @noRd
#' @importFrom stats arima.sim
#'
#' @author Asael Alonzo Matamoros.
#'
#' @references
#' Nieto-Reyes, A., Cuesta-Albertos, J. & Gamboa, F. (2014). A random-projection
#' based test of Gaussianity for stationary processes. \emph{Computational
#' Statistics & Data Analysis, Elsevier}, vol. 75(C), pages 124-141.
#'
#' @examples
#' # rejection rate for a Gaussian AR(1) process with phi = 0.5 and Lobato's test
#' nortsTest:::rejection_rate(reps = 1,n = 100,phi = 0.5)
#'
#'
rejection_rate = function(reps = 1000, n = 100, dist = rnorm, seed = NULL, htest = "lobato",
                          alpha = 0.05, k = 1, phi = 0.5, sd = 1, ...){

  if (!is.null(seed))
    set.seed(seed)

  start_time = Sys.time()
  pvalues = parallel::mclapply(1:reps, FUN = function(i){

    x = stats::arima.sim(n = n, n.start = 500, list(ar = phi), rand.gen = dist, ...)

    if(htest == "vavra"){
      suppressWarnings(vavra.test(x, reps = 100, h = 50)$p.value)
    }
    else if(htest == "lobato"){
      suppressWarnings(lobato.test(x)$p.value)
    }
    else if(htest == "epps"){
      suppressWarnings(epps.test(x)$p.value)
    }
    else if(htest == "rp"){
      suppressWarnings(rp.test(x, k = k)$p.value)
    }
    else if(htest == "epps_bootstrap"){
      suppressWarnings(epps_bootstrap.test(x, reps = 250, h = 50)$p.value)
    }
    else if(htest == "lobato_bootstrap"){
      suppressWarnings(lobato_bootstrap.test(x, reps = 250, h = 50)$p.value)
    }
    else if(htest == "elbouch"){
      suppressWarnings(elbouch.test(y = x)$p.value)
    }
  })
  end_time = Sys.time()

  rr = c(mean(unlist(pvalues) < alpha), difftime(end_time, start_time, units="secs"))
  return(rr)
}
#' Rejection rates table.
#'
#' Performs a rejection rate estimates tables for different stationary AR(1) processes, simulated
#' with different overaggressive coefficients and distributions, the used distributions are normal
#' \code{X ~ N(0,1)} log normal\code{X ~ logN(0,1)},t student \code{X ~ t(df = 7)}, Chi square
#' \code{X ~ Chisq(10)}, and a beta \code{X ~ beta(7,1)} distributions.
#'
#' @param reps an integer with the total bootstrap repetitions.
#' @param n an integer with the length of the simulated AR(1) process.
#' @param phi a real vectors with the auto regressive coefficients to simulate an stationary AR(1) process.
#' @param seed An optional \code{\link[=set.seed]{seed}} to use.
#' @param htest A character string naming the desired test for checking normality. Valid values are
#' \code{"epps"} for the Epps, \code{"lobato"} for Lobato and Velasco's,\code{"vavra"} for the Psaradakis
#' and  Vavra,  and \code{"rp"} for the random projections. The default value is \code{"lobato"} test.
#' @param alpha  Level of the test, possible values range from 0.01 to 0.1. By default \code{alpha = 0.05}.
#' @param k an integer with the number of random projections to be used, by default \code{k = 2}.
#'
#' @return a matrix with rejection rates estimates, where the columns represent different \code{AR(1)}
#' coefficient values, and the row different distributions.
#'
#' @noRd
#'
#' @author Asael Alonzo Matamoros.
#'
#' @references
#' Nieto-Reyes, A., Cuesta-Albertos, J. & Gamboa, F. (2014). A random-projection
#' based test of Gaussianity for stationary processes. \emph{Computational
#' Statistics & Data Analysis, Elsevier}, vol. 75(C), pages 124-141.
#'
#' @examples
#' # rejection rate for a Gaussian AR(1) process with phi = 0.5 and Lobato's test
#' nortsTest:::rejection_table(reps = 1,n = 100)
#'
rejection_table = function(reps = 1000, n = 100, phi = c(-0.4,-0.25,0.0,0.25,0.4),
                           seed = NULL, htest = "lobato", alpha = 0.05, k = 1){

  if (!is.null(seed))
    set.seed(seed)

  df = matrix(0, nrow =5, ncol = length(phi) + 1)
  row.names(df) = c("N","logN","t_3","Chisq_10","Gamma(1,1)")
  colnames(df) = c(phi, "max-time")

  #normal
  x1 = parallel::mclapply(phi, FUN = function(i){
    rejection_rate(reps = reps,n = n,phi = i,dist = rnorm,
                   htest = htest,alpha = alpha,k = k,sd = 1,seed = seed)
  })

  x1 = matrix(unlist(x1),nrow = 2,byrow = FALSE)
  df[1, ] = c(x1[1,],max(x1[2,]))

  #logN
  x1 = parallel::mclapply(phi, FUN = function(i){
    rejection_rate(reps = reps,n = n,phi = i, dist = rlnorm,
                   htest = htest,alpha = alpha,k = k,sdlog = 1,seed = seed)
  })

  x1 = matrix(unlist(x1), nrow = 2, byrow = FALSE)
  df[2, ] = c(x1[1, ], max(x1[2, ]))

  #t_3
  x1 = parallel::mclapply(phi, FUN = function(i){
    rejection_rate(reps = reps, n = n, phi = i, dist = rt, df = 3,
                   htest = htest, alpha = alpha, k = k, seed = seed)
  })

  x1 = matrix(unlist(x1), nrow = 2, byrow = FALSE)
  df[3, ] = c(x1[1, ], max(x1[2, ]))

  #chi_10
  x1 = parallel::mclapply(phi, FUN = function(i){
    rejection_rate(reps = reps, n = n, phi = i, dist = rchisq, df = 10,
                   htest = htest, alpha = alpha, k = k, seed = seed)
  })

  x1 = matrix(unlist(x1), nrow = 2, byrow = FALSE)
  df[4, ] = c(x1[1, ], max(x1[2, ]))

  #gamma(1,1)

  x1 = parallel::mclapply(phi, FUN = function(i){
    rejection_rate(reps = reps, n = n,phi = i, dist = rgamma,shape = 1,
                   htest = htest, alpha = alpha, k = k, seed = seed)
  })

  x1 = matrix(unlist(x1), nrow = 2, byrow = FALSE)
  df[5, ] = c(x1[1, ], max(x1[2, ]))

  return(df)
}