我如何在 R 中为多个产品做自动 arima？

Posted 2023-03-25

【中文标题】我如何在 R 中为多个产品做自动 arima？

【英文标题】：how can i do auto arima for multiple products in R?

【发布时间】：2021-11-15 01:59:46

【问题描述】：

我正在 R 中创建自动 arima 模型来预测我的需求。我为 1 个产品及其工作做这件事。我以 xlsx 格式导出，分列：

Sku（产品），

日期预测（未来 3 个月）

点预测，低 95% 和高 95%。

我的代码是：

ps：葡萄牙语中的变量名，因为我来自巴西。

bdvendas <- read.csv("Pedidos+PedidosItem.csv", header = T, sep = ";")

vendas <- bdvendas %>% 
  dplyr::select(dataPedido,SkuRaiz,quantidadeItemReal)

vendas$dataPedido <- dmy(vendas$dataPedido)

vendas <- subset(vendas, vendas$dataPedido > "2018-12-31")
vendas <- subset(vendas, vendas$SkuRaiz!="")
vendas <- na.omit(vendas)

teste <-  data.frame(as.yearmon(vendas$dataPedido))
teste <- cbind(vendas,teste)
names(teste)[1:length(teste)] <- c("dataPedido","SkuRaiz","Pedidos","MesPedido")

vendas <- teste %>% 
  group_by(MesPedido,SkuRaiz) %>% 
  summarise(Pedidos = sum(Pedidos))

analisesku <- vendas %>% 
  filter(SkuRaiz == 1081) ## <- HERE I SELECT MY PRODUCT

analisesku <- analisesku[-length(analisesku$Pedidos),]

ano_inicial <- as.numeric(format(analisesku$MesPedido,'%Y'))[1]
mes_inicial <- as.numeric(format(analisesku$MesPedido,'%m'))[1]

ano_final <- as.numeric(format(analisesku$MesPedido,'%Y'))[length(analisesku$MesPedido)]
mes_final <- as.numeric(format(analisesku$MesPedido,'%m'))[length(analisesku$MesPedido)]

tsbanco <- ts(analisesku$Pedidos, start = c(ano_inicial,mes_inicial), end = c(ano_final,mes_final), frequency = 12)
autoplot(tsbanco)

modelo <- auto.arima(tsbanco, stepwise = FALSE, approximation = FALSE, trace = TRUE)
previsao <- forecast(modelo, h=2, level = c(95))
print(previsao)
autoplot(previsao)
accuracy(previsao)
output <- print(summary(previsao))
output <- cbind(analisesku$SkuRaiz[1],output)
names(output) <- c("SkuRaiz","pointForecast","low95","high95")
mesprevisao <- data.frame(seq(as.Date(Sys.Date()), by = "month", length = 3))
names(mesprevisao) <- "mesPrevisao"
output <- cbind(mesprevisao,output)

write.table(output, file = "previsao.csv", sep = ";", dec = ',', row.names = F, col.names = T)

效果很好。

但是，我的问题是：我需要自动为多个产品（大约 3000 个产品）执行此操作。

ps：每个产品都有独特的系列。他们是独立的。

我该怎么做？我需要使用循环或类似的东西吗？

【问题讨论】：

将 fable 包与 ARIMA 函数一起使用（它实现了与 forecast 包中的 auto.arima 相同的算法）。它旨在处理多个时间序列，并与您已经在使用的 tidyverse 集成。

【参考方案1】：

您没有提供任何数据，因此我将模拟一些数据并逐步演示如何预测多个时间序列。

负载预测库

library(forecast)

让我们从 ARIMA 模型中模拟 5 个时间序列

bts <- ts(dplyr::tibble(AA = arima.sim(list(order=c(1,0,0), ar=.5),
                                       n=100, mean = 12),
                        AB = arima.sim(list(order=c(1,0,0), ar=.5),
                                       n=100, mean = 12),
                        AC = arima.sim(list(order=c(1,0,0), ar=.5),
                                       n=100, mean = 11),
                        BA = arima.sim(list(order=c(1,0,0), ar=.5),
                                       n=100, mean = 10),
                        BB = arima.sim(list(order=c(1,0,0), ar=.5),
                                       n=100, mean = 14)), start = c(2000, 1),
          frequency = 12)

绘制所有 ts

autoplot(bts)

将模型拟合到所有 ts

fit <- sapply(bts, FUN = auto.arima, simplify = FALSE, USE.NAMES = TRUE,
              # auto.arima arguments
              max.p = 5,
              max.q = 5,
              max.P = 2,
              max.Q = 2 # other arguments passed to auto arima
              )

预测所有模型

fc <- sapply(fit, FUN = forecast, simplify = FALSE, USE.NAMES = TRUE,
             h = 12 # forecast horizon
             # other arguments passed to forecast
             )

这个简单的函数将帮助我们在列表中获得均值、下限或上限预测

get_value <- function(x, type = c("mean", "lower", "upper"), 
                      level = c(80, 95))
  if(type == "mean")
    out <- x[["mean"]]
  
  if(type == "lower")
    if(level == 80)
      out <- x[["lower"]][,1]
    
    if(level == 95)
      out <- x[["lower"]][,2]
    
  
  if(type == "upper")
    if(level == 80)
      out <- x[["upper"]][,1]
    
    if(level == 95)
      out <- x[["upper"]][,2]
    
  
  return(out)

获取平均预测

point_forecast <- sapply(fc, FUN = get_value, simplify = TRUE, 
                         USE.NAMES = TRUE,
                         type = "mean")

获得具有 95 % 置信区间的上限值

fc_upper_95 <- sapply(fc, FUN = get_value, simplify = TRUE, 
                         USE.NAMES = TRUE,
                         type = "upper", level = 95)

用 80 % 的置信区间获取上限值

fc_upper_80 <- sapply(fc, FUN = get_value, simplify = TRUE, 
                      USE.NAMES = TRUE,
                      type = "upper", level = 80)
                     

由于您有许多时间序列，最好并行拟合模型以有效利用计算资源

library(parallel)

n_cores <- parallel::detectCores()-1 # number of cores in your machine -1 core

cl <- makeCluster(n_cores)

fit_par <- parallel::parSapply(cl, bts, FUN = auto.arima, 
                               simplify = FALSE, USE.NAMES = TRUE,
                               # auto.arima arguments
                               max.p = 5,
                               max.q = 5,
                               max.P = 2,
                               max.Q = 2)

fc_par <- parallel::parSapply(cl, fit_par, FUN = forecast, simplify = FALSE, 
                              USE.NAMES = TRUE,
                              h = 12
                              # other arguments passed to forecast
                              )

point_forecast <- parallel::parSapply(cl, fc_par, FUN = get_value, 
                                      simplify = TRUE, USE.NAMES = TRUE, 
                                      type = "mean")