如何将搜索AR(1)模拟匹配种子的R代码重写为单个可调用函数

背景问题 #

使用arima.sim()函数模拟AR(1)时间序列后，用forecast::auto.arima()识别模型时，经常出现识别出的阶数、系数和模拟时的设定值不匹配的情况。
比如设置种子为1时，识别结果和设定的AR(1)、系数0.8完全不符：

set.seed(1)
ar1 <- arima.sim(n=10, model=list(ar=0.8, order=c(1, 0, 0)), sd=1)
auto_ar <- forecast::auto.arima(ar1, ic="aicc")
auto_ar

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

#sigma^2 estimated as 1.008:  log likelihood=-14.23
#AIC=30.46   AICc=30.96   BIC=30.76

而设置种子为289805时，识别结果和设定完全匹配：

set.seed(289805)
ar1 <- arima.sim(n=10, model=list(ar=0.8, order=c(1, 0, 0)), sd=1)
auto_ar <- forecast::auto.arima(ar1, ic="aicc")
auto_ar

#Series: ar1 
#ARIMA(1,0,0) with zero mean 

#Coefficients:
#         ar1
#      0.8000
#s.e.  0.1458

#sigma^2 estimated as 1.849:  log likelihood=-17.25
#AIC=38.49   AICc=40.21   BIC=39.1

原有种子搜索代码 #

原有的非封装版本实现如下：

FUN <- function(i) {
  set.seed(i)
  ar1 <- arima.sim(n=100, model=list(ar=0.8, order=c(1, 0, 0)), sd=1)
  ar2 <- auto.arima(ar1, ic="aicc")
  cf <- ar2$coef
  ## case handling
  if (length(cf) == 0) rep(NA, 2)  ## sometimes result is `character(0)` -> NA
  else if (substr(cf[1], 1, 5) %in% "0.800") c(cf, i)  ## hit, that's what we want
  else rep(NA, 2)  ## all other cases -> NA
}

R <- 1e3  ## this would be your 1e5
seedv <- 1:R  ## or use custom seed vector

library(parallel)
cl <- makeCluster(detectCores() - 1)  ## for all cores remove `- 1`
clusterExport(cl, c("FUN"), envir=environment())
clusterEvalQ(cl, suppressPackageStartupMessages(library(forecast)))

res <- `colnames<-`(t(parSapply(cl, seedv, "FUN")), c("cf", "seed"))

stopCluster(cl)

修正后的单函数实现 #

你编写的版本存在几个问题：不需要将内部校验函数作为入参、系数匹配时错误将参数名作为字符串匹配、并行计算时没有把外部参数导出到子进程。修正后的实现如下：

search_ar1_seed <- function(n, ar, sd, target_ar, R) {
  # 内部种子校验函数
  check_seed <- function(i) {
    set.seed(i)
    # 模拟AR(1)序列
    ar1 <- arima.sim(n = n, model = list(ar = ar, order = c(1, 0, 0)), sd = sd)
    # 自动识别模型
    ar_fit <- forecast::auto.arima(ar1, ic = "aicc")
    cf <- ar_fit$coef
    # 异常情况处理：没有识别到ar1系数直接返回空
    if (length(cf) == 0 || is.null(cf["ar1"])) {
      return(rep(NA, 2))
    }
    # 匹配目标系数精度（保留前4位小数）
    target_str <- substr(as.character(target_ar), 1, 5)
    actual_str <- substr(as.character(cf["ar1"]), 1, 5)
    if (actual_str == target_str) {
      return(c(ar1_coef = unname(cf["ar1"]), seed = i))
    } else {
      return(rep(NA, 2))
    }
  }
  
  # 并行计算初始化
  library(parallel)
  cl <- makeCluster(detectCores() - 1)
  # 导出所有需要的变量到集群节点
  clusterExport(cl, c("n", "ar", "sd", "target_ar", "check_seed"), envir = environment())
  # 节点加载forecast包
  clusterEvalQ(cl, suppressPackageStartupMessages(library(forecast)))
  
  # 批量校验种子
  res <- parSapply(cl, 1:R, check_seed)
  # 整理结果，过滤无效值
  res_df <- na.omit(as.data.frame(t(res), stringsAsFactors = FALSE))
  rownames(res_df) <- NULL
  
  # 关闭集群
  stopCluster(cl)
  return(res_df)
}

参数说明 #

• n：模拟序列的样本量
• ar：模拟AR(1)序列时设定的φ系数值
• sd：模拟序列的误差项标准差
• target_ar：需要auto.arima识别结果匹配的φ系数值
• R：从1开始搜索的最大种子数

调用示例 #

# 搜索样本量为10、系数匹配0.8的种子，搜索范围1到30万
result <- search_ar1_seed(n = 10, ar = 0.8, sd = 1, target_ar = 0.8, R = 300000)
# 查看匹配到的种子
print(result)

内容的提问来源于stack exchange，提问作者Daniel James