Nothing
#' 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)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.