I am trying to write a function with arima.sim() function in R such that I only need to input the variables of any of the ARMA families to get the seeds that will produce my type of ARMA simulation.
Two out of my three (2/3) cases are true when I use the following
arma_sim_search <- function(a, z, n, ar11, ma11, ar22, ma22, ar33, ma33, p, d, q, sd = 1, j1, k1, j2, k2, j3, k3, arr1, maa1, arr2, maa2, arr3, maa3){
output <- if (p == 1) {
ar1_sim_search <- function(a, z, n, ar11, p, d, q, sd = 1, j1, arr1){
future::plan(future::multisession)
n_cores <- parallel::detectCores()
cl <- parallel::makeCluster(n_cores)
doParallel::registerDoParallel(cores = n_cores)
message('processing...')
`%dopar%` <- foreach::`%dopar%`
i <- a:z
res <- foreach::foreach(i = a:z, .packages = c('foreach', 'forecast')) %dopar% {
set.seed(i)
mod <- stats::arima.sim(n = n, model = list(ar = c(ar11), order = c(p, d, q)), sd = sd)
best.mod <- forecast::auto.arima(mod, ic = "aicc")
(cf <- best.mod$coef)
if (length(cf) == 0) {
rep(NA, 2)
} else if (all(grepl(c("ar1|intercept"), names(cf))) &
substr(cf["ar1"], 1, j1) %in% arr1) {
c(cf, seed = I)
} else {
rep(NA, 2)
}
}
message(' done!\n')
res1 = res[!sapply(res, anyNA)]
parallel::stopCluster(cl)
options(max.print = .Machine$integer.max)
res2 <- tibble::tibble(Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x)))))
res2[order(res2$seed), ]
res2 <- Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x))))
res2[order(res2$seed), ]
}
ar1_sim_search(a = a, z = z, n = n, p = 1, d = d, q = q, ar11 = ar11, sd = sd, j1 = j1, arr1 = arr1)
#######################################################################
} else if (p == 1 && q == 1) {
arma1_sim_search <- function(a, z, n, ar11, ma11, p, d, q, sd = 1, j1, k1, arr1, maa1){
future::plan(future::multisession)
n_cores <- parallel::detectCores()
cl <- parallel::makeCluster(n_cores)
doParallel::registerDoParallel(cores = n_cores)
message('processing...')
`%dopar%` <- foreach::`%dopar%`
i <- a:z
res <- foreach::foreach(i = a:z, .packages = c('foreach', 'forecast')) %dopar% {
set.seed(i)
mod <- stats::arima.sim(n = n, model = list(ar = ar11, ma = ma11, order = c(p, d, q)), sd = sd)
best.mod <- forecast::auto.arima(mod, ic = "aicc")
(cf <- best.mod$coef)
if (length(cf) == 0) {
rep(NA, 2)
} else if (all(grepl(c("ar1|ma1|intercept"), names(cf))) &
substr(cf["ar1"], 1, j1) %in% arr1 & substr(cf["ma1"], 1, k1) %in% maa1) {
c(cf, seed = I)
} else {
rep(NA, 2)
}
}
message(' done!\n')
res1 = res[!sapply(res, anyNA)]
parallel::stopCluster(cl)
options(max.print = .Machine$integer.max)
res2 <- tibble::tibble(Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x)))))
res2[order(res2$seed), ]
res2 <- Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x))))
res2[order(res2$seed), ]
}
arma1_sim_search(a = a, z = z, n = n, p = 1, d = d, q = 1, ar11 = ar11, ma11 = ma11, sd = sd, j1 = j1, k1 = k1, arr1 = arr1, maa1 = maa1)
##############################################################
} else {
ma1_sim_search <- function(a, z, n, ma11, p, d, q, sd = 1, k1, maa1){
future::plan(future::multisession)
n_cores <- parallel::detectCores()
cl <- parallel::makeCluster(n_cores)
doParallel::registerDoParallel(cores = n_cores)
message('processing...')
`%dopar%` <- foreach::`%dopar%`
i <- a:z
res <- foreach::foreach(i = a:z, .packages = c('foreach', 'forecast')) %dopar% {
set.seed(i)
mod <- stats::arima.sim(n = n, model = list(ma = c(ma11), order = c(p, d, q)), sd = sd)
best.mod <- forecast::auto.arima(mod, ic = "aicc")
(cf <- best.mod$coef)
if (length(cf) == 0) {
rep(NA, 2)
} else if (all(grepl(c("ma1|intercept"), names(cf))) &
substr(cf["ma1"], 1, k1) %in% maa1) {
c(cf, seed = I)
} else {
rep(NA, 2)
}
}
message(' done!\n')
res1 = res[!sapply(res, anyNA)]
parallel::stopCluster(cl)
options(max.print = .Machine$integer.max)
res2 <- tibble::tibble(Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x)))))
res2[order(res2$seed), ]
res2 <- Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x))))
res2[order(res2$seed), ]
}
ma1_sim_search(a = a, z = z, n = n, p = p, d = d, q = 1, ma11 = ma11, sd = sd, k1 = k1, maa1 = maa1)
}
##############################################################################
return(output)
}
The first case works well with this output.
arma_sim_search(a = 1, z = 1000, n = 25, p = 1, d = 0, q = 0, ar11 = 0.8, sd = 1, j1 = 3, arr1 = "0.8")
#processing...
#done!
# ar1 seed
#64 0.8835177 16
#2 0.8026257 30
#9 0.8115264 42
#19 0.8246183 43
#67 0.8971222 49
The second case is giving me error.
arma_sim_search(a = 1, z = 1000, n = 25, p = 1, d = 0, q = 1, ar11 = 0.4, ma11 = 0.5, sd = 1, j1 = 3, k1 = 3, arr1 = "0.4", maa1 = "0.5")
I got this error message: Error in { : task 1 failed - "inconsistent specification of 'ma' order"
The third case works well like the first.
arma_sim_search(a = 1, z = 1000, n = 25, p = 0, d = 0, q = 1, ma11 = 0.8, sd = 1, k1 = 3, maa1 = "0.8")
#processing...
#done!
# ma1 seed intercept
#81 0.8551335 24 NA
#14 0.8090642 28 NA
#12 0.8080051 29 NA
#43 0.8314141 33 NA
This is the skeletal image of my function: the first case is when the expression is an AR1 with p == 1 the second is when the expression is an ARMA11 with p ==1 && q ==1 and the last is a case when the expression is an MA1 case when q == 1.
arma_sim_search <- function(a, z, n, ar11, ma11, ar22, ma22, ar33, ma33, p, d, q, sd = 1, j1, k1, j2, k2, j3, k3, arr1, maa1, arr2, maa2, arr3, maa3){
output <- if (p == 1) {
ar1_sim_search <- function(a, z, n, ar11, p, d, q, sd = 1, j1, arr1){
AR1 EXPRESION
}
ar1_sim_search(a = a, z = z, n = n, p = 1, d = d, q = q, ar11 = ar11, sd = sd, j1 = j1, arr1 = arr1)
########################################################
} else if (p == 1 && q == 1) {
arma1_sim_search <- function(a, z, n, ar11, ma11, p, d, q, sd = 1, j1, k1, arr1, maa1){
ARMA11 EXPRESION
}
arma1_sim_search(a = a, z = z, n = n, p = 1, d = d, q = 1, ar11 = ar11, ma11 = ma11, sd = sd, j1 = j1, k1 = k1, arr1 = arr1, maa1 = maa1)
##################################################
} else {
ma1_sim_search <- function(a, z, n, ma11, p, d, q, sd = 1, k1, maa1){
MA1 EXPRESION
}
ma1_sim_search(a = a, z = z, n = n, p = p, d = d, q = 1, ma11 = ma11, sd = sd, k1 = k1, maa1 = maa1)
}
##############################################################################
return(output)
}
arma_sim_search(a = 1, z = 1000, n = 25, p = 1, d = 0, q = 1, ar11 = 0.4, ma11 = 0.5, sd = 1, j1 = 3, k1 = 3, arr1 = "0.4", maa1 = "0.5")
Meanwhile, when I ran the function of case two alone it ran well.
arma1_sim_search <- function(a, z, n, ar11, ma11, p, d, q, sd = 1, j1, k1, arr1, maa1){
future::plan(future::multisession)
n_cores <- parallel::detectCores()
cl <- parallel::makeCluster(n_cores)
doParallel::registerDoParallel(cores = n_cores)
message('processing...')
`%dopar%` <- foreach::`%dopar%`
i <- a:z
res <- foreach::foreach(i = a:z, .packages = c('foreach', 'forecast')) %dopar% {
set.seed(i)
mod <- stats::arima.sim(n = n, model = list(ar = ar11, ma = ma11, order = c(p, d, q)), sd = sd)
best.mod <- forecast::auto.arima(mod, ic = "aicc")
(cf <- best.mod$coef)
if (length(cf) == 0) {
rep(NA, 2)
} else if (all(grepl(c("ar1|ma1|intercept"), names(cf))) &
substr(cf["ar1"], 1, j1) %in% arr1 & substr(cf["ma1"], 1, k1) %in% maa1) {
c(cf, seed = I)
} else {
rep(NA, 2)
}
}
message(' done!\n')
res1 = res[!sapply(res, anyNA)]
parallel::stopCluster(cl)
options(max.print = .Machine$integer.max)
res2 <- tibble::tibble(Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x)))))
res2[order(res2$seed), ]
res2 <- Reduce(function(...) merge(..., all = T), lapply(res1, function(x) as.data.frame(t(x))))
res2[order(res2$seed), ]
}
arma1_sim_search(a = 1, z = 1000, n = 50, p = 1, d = 0, q = 1, ar11 = 0.2, ma11 = 0.7, sd = 1, j1 = 3, k1 = 3, arr1 = 0.2, maa1 = 0.7)
#processing...
# done!
# ar1 ma1 seed
#5 0.2887975 0.7388537 212
#4 0.2871108 0.7819748 229
#3 0.2858437 0.7425638 249
#2 0.2728310 0.7659739 310
#1 0.2574142 0.7674407 935
