how to count how many times an arima order is not true before the first true order is obtained in r for different combo of arima simulation

Question

Most times one runs arima.sim() function to simulate a particular order of arima mosel but when one check such simulated time series data through auto.arima() function, it will not often time be the same order of ARIMA one desire and specified in the arima.sim().

In my bid to know how many times one may need to run arima.sim() function for a different combination of its parameter (sample size, standard deviation and coefficient of the model) before obtaining the true order of the model sought for, I want this R script to count how many time it will run an arima.sim() before it get the exert ARIMA-order specified in the arima.sim() function.

**Here is my trial**

library(forecast)
N <- c(10, 20, 30)
SD <- c(1, 2, 3, 4, 5) ^ 2
phi <- c(0.2, 0.4, 0.6)

## generate all combos
all_combos <- expand.grid(N = N, SD = SD, phi = phi)

## create function
set.seed(123)
res2 <- by(all_combos, all_combos["N"], function(DF){
  res <- mapply(function(N, SD, phi){
    cnt <- 0
    repeat {
      x <- arima.sim(n=N, model = list(ar=phi, order = c(1, 0, 0)), sd = SD)
      if(all(arimaorder(auto.arima(x), ic = "aicc"))) != c(1, 0, 0) cnt <- cnt + 1){
      }
        {else(all(arimaorder(auto.arima(x), ic = "aicc"))) == c(1, 0, 0) cnt <- cnt + 1)}
        break
    }
    cnt
  }, DF[["N"]], DF[["SD"]], DF[["phi"]])
  names(res) <- paste("SD", DF[["SD"]], "phi", DF[["phi"]], sep = "-")
  res
})
res2

I am interested in knowing how many trials of arima.sim() will one make before obtaining the first ARIMA(1, 0, 0).

Answer 1

It seems odd to me that you are running by + mapply. I think only mapply is enough. Moreover, arimaorder does not have ic argument, maybe you meant to use it for auto.arima function.

Since you want to know how many trials are needed to get c(1, 0, 0), I add an additional column (index) which is the row number in all_combos. As soon as you get output as c(1, 0, 0) the loop is broken and it prints the index. The code doesn't run for rest of the combinations.

library(forecast)
N <- c(10, 20, 30)
SD <- c(1, 2, 3, 4, 5) ^ 2
phi <- c(0.2, 0.4, 0.6)

## generate all combos
all_combos <- expand.grid(N = N, SD = SD, phi = phi)
all_combos$index <- seq_len(nrow(all_combos))

mapply(function(N, phi, SD, index) {
  x <- with(all_combos, arima.sim(n=N[1], 
             model = list(ar=phi[1], order = c(1, 0, 0)), sd = SD[1]))
  if(all(arimaorder(auto.arima(x, ic = "aicc")) == c(1, 0, 0))) {
    print(index)
    break
  }
}, all_combos$N, all_combos$SD, all_combos$phi, all_combos$index)