Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Welcome To Ask or Share your Answers For Others

Categories

0 votes
816 views
in Technique[技术] by (71.8m points)

search - Searching for the Right AR(2) Seed Like I Did for AR(1) in R

Using an arima.sim() function to simulate time series data that follows a particular ARIMA model requires a lot of trials (especially for a ridiculously small sample size) of changing set.seed value like the example bellow:

library(forecast)
set.seed(1)
ar1 <- arima.sim(n = 15, model=list(ar=0.2, order = c(1, 0, 0)), sd = 1)
(ar2 <- auto.arima(ar1, ic ="aicc"))

set.seed(3)
ar1 <- arima.sim(n = 15, model=list(ar=0.2, order = c(1, 0, 0)), sd = 1)
(ar2 <- auto.arima(ar1, ic ="aicc"))

set.seed(3)
ar1 <- arima.sim(n = 15, model=list(ar=0.2, order = c(1, 0, 0)), sd = 1)
(ar2 <- auto.arima(ar1, ic ="aicc"))

until I get my desired result I have found an R code that prints out the right seed to set for AR(1):

library(future.apply)
FUN <- function(i) { 
  set.seed(i) 
  ar1 <- arima.sim(n=15, model=list(ar=0.6, order=c(1, 0, 0)), sd=1) 
  ar2 <- auto.arima(ar1, ic="aicc") 
  (cf <- ar2$coef) 
  if (length(cf) == 0) { 
    rep(NA, 2) 
  } 
  else if (all(grepl(c("ar1|intercept"), names(cf))) & 
           substr(cf["ar1"], 1, 6) %in% "0.6000") { 
    c(cf, seed=i) 
  } 
  else { 
    rep(NA, 2) 
  } 
}

R <- 1e4 
seedv <- 1:R 

library(parallel) 
cl <- makeCluster(detectCores() - 1 + 1) 
clusterExport(cl, c("FUN"), envir=environment()) 
clusterEvalQ(cl, suppressPackageStartupMessages(library(forecast))) 

res <- parLapply(cl, seedv, "FUN") 

(res1 <- res[!sapply(res, anyNA)]) 

stopCluster(cl) 

res2 <- Reduce(function(...) merge(..., all=T), lapply(res1, 
function(x) as.data.frame(t(x)))) 
res2[order(res2$seed), ]

which produces this result:

        #ar1 seed
#1 0.8000417 4195

which when I check as bellow:

ar1 <- arima.sim(n=15, model=list(ar=0.6, order=c(1, 0, 0)), sd=1) 
(ar2 <- forecast::auto.arima(ar1, ic="aicc"))

it will be true

Series: ar1 ARIMA(1,0,0) with zero mean Coefficients: ar1 0.6000 s.e. 0.1923

sigma^2 estimated as 2.051: log likelihood=-26.38 AIC=56.75 AICc=57.75 BIC=58.17

What I want

I want a case of AR(2) Bellow is my trial:

FUN <- function(i) { 
  set.seed(i) 
  ar1 <- arima.sim(n=15, model=list(ar=c(-0.9, -0.9), order=c(2, 0, 0)), sd=1) 
  ar2 <- auto.arima(ar1, ic="aicc") 
  (cf <- ar2$coef) 
  if (length(cf) == 0) { 
    rep(NA, 2) 
  } 
  else if (all(grepl(c("ar1|ar2|intercept"), names(cf))) & 
           substr(cf["ar1"], 1, 5) %in% "-0.90" & substr(cf["ar2"], 1, 5) %in% "-0.90") { 
    c(cf, seed=i) 
  } 
  else { 
    rep(NA, 2) 
  } 
}

R <- 2e5 
seedv <- 1e5:R 

library(parallel) 
cl <- makeCluster(detectCores() - 1 + 1) 
clusterExport(cl, c("FUN"), envir=environment()) 
clusterEvalQ(cl, suppressPackageStartupMessages(library(forecast))) 

res <- parLapply(cl, seedv, "FUN") 

(res1 <- res[!sapply(res, anyNA)]) 

stopCluster(cl) 

res2 <- Reduce(function(...) merge(..., all=T), lapply(res1, 
function(x) as.data.frame(t(x)))) 
res2[order(res2$seed), ]

The search stops here ###################################################################### I confirm my result as bellow:

set.seed(1509)
ar1 <- arima.sim(n=20, model=list(ar=c(-0.9, -0.9), order=c(2, 0, 0)), sd=1)
library(forecast)
(ar2 <-auto.arima(ar1, ic="aicc"))

My output on coefficients of AR(2) are not exactly -0.90, -0.90:

Series: ar1 ARIMA(2,0,0) with zero mean

Coefficients: ar1 ar2 -0.9369 -0.8997 s.e. 0.0940 0.0752

sigma^2 estimated as 0.9048: log likelihood=-28.12 AIC=62.24 AICc=63.74 BIC=65.23

question from:https://stackoverflow.com/questions/65894592/searching-for-the-right-ar2-seed-like-i-did-for-ar1-in-r

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
Welcome To Ask or Share your Answers For Others

1 Reply

0 votes
by (71.8m points)
Waitting for answers

与恶龙缠斗过久,自身亦成为恶龙;凝视深渊过久,深渊将回以凝视…
OGeek|极客中国-欢迎来到极客的世界,一个免费开放的程序员编程交流平台!开放,进步,分享!让技术改变生活,让极客改变未来! Welcome to OGeek Q&A Community for programmer and developer-Open, Learning and Share
Click Here to Ask a Question

1.4m articles

1.4m replys

5 comments

57.0k users

...