R/rejection_rates.R

Defines functions rejection_table rejection_rate

#' 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}
#' is used
#' @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 autoregressive coefficent to simulate an stationary AR(1) process.
#' @param sd optional value for the standard desviation 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 lobatos 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 = 64,phi = 0.5,sd = 1,...){

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

  pvalues = rep(0,reps)

  for(i in 1:reps){
    x = stats::arima.sim(n=n,n.start = 500,list(ar = phi),rand.gen = dist,...)
    if(htest == "vavra"){
      pvalues[i] = vavra.test(x,reps = 100,h = 50)$p.value
    }
    else if(htest == "lobato"){
      pvalues[i] = suppressWarnings(lobato.test(x)$p.value)
    }
    else if(htest == "epps"){
      pvalues[i] = suppressWarnings(epps.test(x)$p.value)
    }
    else if(htest == "rp"){
      pvalues[i] = suppressWarnings(rp.test(x,k = k)$p.value)
    }
  }

  rr = mean(pvalues < alpha)
  return(rr)
}
#' Rejection rates table
#'
#' Performs a rejection rate estimates tables for different stationary AR(1) processes, simulated
#' with different autorregressive coefficents 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 autoregressive coefficents 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}
#' is used.
#' @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 lobatos 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 = 64){

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

  nphi = length(phi)
  df = matrix(0,nrow =5,ncol = nphi)
  row.names(df) = c("N","logN","t_3","Chisq_10","beta(7,1)")
  colnames(df) = phi
  #normal
  for(i in 1:nphi) df[1,i] = rejection_rate(reps = reps,n=n,phi = phi[i],dist = rnorm,
                                         htest = htest,alpha = alpha,k = k,sd = 1,seed = seed)
  #logN
  for(i in 1:nphi) df[2,i] = rejection_rate(reps = reps,n=n,phi = phi[i],dist = rlnorm,
                                         htest = htest,alpha = alpha,k = k,sdlog = 1,seed = seed)
  #t_3
  for(i in 1:nphi) df[3,i] = rejection_rate(reps = reps,n=n,phi = phi[i],dist = rt,df = 3,
                                         htest = htest,alpha = alpha,k = k,seed = seed)
  #chi_10
  for(i in 1:nphi) df[4,i] = rejection_rate(reps = reps,n=n,phi = phi[i],dist = rchisq,df = 10,
                                         htest = htest,alpha = alpha,k = k,seed = seed)
  #beta(7,1)
  for(i in 1:nphi) df[5,i] = rejection_rate(reps = reps,n=n,phi = phi[i],dist = rbeta,shape1 = 7,shape2 = 1,
                                         htest = htest,alpha = alpha,k = k,seed = seed)

  return(df)
}

Try the nortsTest package in your browser

Any scripts or data that you put into this service are public.

nortsTest documentation built on July 27, 2020, 5:07 p.m.