使用auto.arima生成30天预测失败,预测周期过长问题排查
问题:ARIMA预测周期超出设定的30天
我有2016年5月至11月三个分组的计数数据(本次仅选三组测试),用auto.arima模型尝试生成30天的计数预测,但实际预测周期延伸到了1月,其中一个分组甚至到了3月,请问问题出在哪?
实现代码
library(tidyverse) library(tidyquant) library(timetk) library(sweep) library(forecast) sub <- structure(list(group = c("group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_1", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_2", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3", "group_3"), date = structure(c(16934, 16947, 16952, 16955, 16959, 16962, 16965, 16968, 16971, 16974, 16977, 16980, 16983, 16986, 16989, 16992, 16995, 16998, 17001, 17004, 17007, 17010, 17013, 17016, 17019, 17022, 17025, 17028, 17031, 17034, 17037, 17040, 17043, 17046, 17049, 17052, 17055, 17058, 17061, 17064, 17067, 17070, 17073, 17076, 17079, 17082, 17085, 17088, 17091, 17094, 17097, 17101, 16963, 16968, 16974, 16977, 16983, 16986, 16989, 16992, 16995, 16998, 17001, 17004, 17007, 17010, 17013, 17016, 17019, 17022, 17025, 17028, 17031, 17034, 17037, 17040, 17043, 17046, 17049, 17052, 17055, 17061, 17066, 17071, 17074, 17079, 17082, 17088, 17093, 17099, 17103, 17108, 17113, 16994, 17001, 17004, 17008, 17012, 17016, 17019, 17022, 17025, 17029, 17032, 17035, 17038, 17042, 17045, 17049, 17052, 17056, 17059, 17062, 17067, 17071, 17075, 17080, 17086, 17092, 17099, 17104, 17108), class = "Date"), count = c(65, 12, 46, 33, 19, 18, 56, 21, 50, 13, 80, 70, 56, 59, 78, 96, 111, 140, 147, 132, 86, 96, 186, 169, 153, 106, 94, 80, 134, 172, 217, 148, 106, 94, 102, 74, 132, 75, 108, 50, 81, 78, 38, 91, 109, 44, 101, 82, 102, 28, 44, 48, 56, 82, 64, 74, 16, 69, 87, 11, 97, 144, 41, 95, 99, 83, 54, 62, 131, 92, 90, 104, 113, 51, 74, 72, 84, 36, 25, 94, 100, 58, 32, 62, 41, 70, 17, 80, 37, 53, 63, 67, 73, 63, 27, 36, 17, 55, 16, 38, 48, 97, 88, 84, 39, 34, 24, 60, 61, 10, 25, 20, 85, 21, 78, 85, 16, 16, 82, 81, 53, 25)), row.names = c(NA, -122L), class = c("tbl_df", "tbl", "data.frame" )) dta <- sub %>% mutate(order = as_date((date))) %>% select(-date) dta_nest <- dta %>% group_by(group) %>% nest() ## Create a daily Date object inds <- seq(min(sub$date, na.rm=T), max(sub$date, na.rm=T), by = "day") # Create the time series data dta_ts <- dta_nest %>% mutate(data.ts = map(.x = data, .f = tk_ts, select = -order, start = c(2022, as.numeric(format(inds[1], "%j")))), freq = 365) # Fit ARIMA dta_fit <- dta_ts %>% mutate(fit.arima = map(data.ts, auto.arima)) # Obtain the augmented fitted and residual values augment_fit_arima <- dta_fit %>% mutate(augment = map(fit.arima, sw_augment, timetk_idx = TRUE, rename_index = "date")) %>% unnest(augment) # Forecast dta_fcast <- dta_fit %>% mutate(fcast.arima = map(fit.arima, forecast, h = 30)) # 30 day forecast dta_fcast_tidy <- dta_fcast %>% mutate(sweep = map(fcast.arima, sw_sweep, fitted = FALSE, timetk_idx = TRUE)) %>% unnest(sweep) # Plot the forecast dta_fcast_tidy %>% ggplot(aes(x = index, y = count, color = key, group = group)) + geom_ribbon(aes(ymin = lo.95, ymax = hi.95), fill = "#D5DBFF", color = NA, size = 0) + geom_ribbon(aes(ymin = lo.80, ymax = hi.80, fill = key), fill = "#596DD5", color = NA, size = 0, alpha = 0.8) + geom_line() + labs(title = "Counts by Group", subtitle = "ARIMA Model Forecasts", x = "", y = "Units") + scale_x_date(date_breaks = "2 weeks", date_labels = "%b %d") + scale_color_tq() + scale_fill_tq() + facet_wrap(~ group, scales = "free_y", ncol = 1) + theme_tq() + theme(axis.text.x = element_text(angle = 45, hjust = 1))
当前预测图

问题根源
核心问题出在时间序列的频率设置和起始时间定义错误:
- 你使用全局的
inds[1](所有分组的最早日期)作为每个分组时间序列的起始点,同时设置freq=365,但每个分组的实际数据起始日期远晚于这个全局起始点,导致时间序列中存在大量未观测的“隐含日期”。 auto.arima会把这些缺失日期视为时间序列的一部分,forecast(h=30)的h参数指的是时间序列的观测步长数,不是实际天数。由于原始数据是间隔观测,30个步长对应的实际天数会远大于30天。
修复方案
需要为每个分组单独构建连续的时间序列,确保观测点与实际日期一一对应,具体修改如下:
修改后的代码
library(tidyverse) library(tidyquant) library(timetk) library(sweep) library(forecast) # 原始数据(保持不变) sub <- structure(list(group = c("group_1", "group_1", ...)), row.names = c(NA, -122L), class = c("tbl_df", "tbl", "data.frame")) # 1. 为每个分组补全每日日期,缺失计数填充为0(可根据业务调整为NA) dta_full <- sub %>% group_by(group) %>% complete(date = seq(min(date), max(date), by = "day"), fill = list(count = 0)) %>% ungroup() # 2. 嵌套分组数据 dta_nest <- dta_full %>% group_by(group) %>% nest() # 3. 为每个分组创建正确的时间序列(基于自身最早日期) dta_ts <- dta_nest %>% mutate( first_date = map(data, ~ min(.$date)), data.ts = map2(.x = data, .y = first_date, .f = function(df, start_date) { tk_ts(df, select = count, start = c(year(start_date), yday(start_date)), freq = 365) }) ) # 4. 拟合ARIMA模型 dta_fit <- dta_ts %>% mutate(fit.arima = map(data.ts, auto.arima)) # 5. 生成30天预测(此时h=30对应30个连续日) dta_fcast <- dta_fit %>% mutate(fcast.arima = map(fit.arima, forecast, h = 30)) # 6. 整理预测结果并绘图 dta_fcast_tidy <- dta_fcast %>% mutate(sweep = map(fcast.arima, sw_sweep, fitted = FALSE, timetk_idx = TRUE)) %>% unnest(sweep) dta_fcast_tidy %>% ggplot(aes(x = index, y = count, color = key, group = group)) + geom_ribbon(aes(ymin = lo.95, ymax = hi.95), fill = "#D5DBFF", color = NA, size = 0) + geom_ribbon(aes(ymin = lo.80, ymax = hi.80, fill = key), fill = "#596DD5", color = NA, size = 0, alpha = 0.8) + geom_line() + labs(title = "分组计数预测", subtitle = "ARIMA模型30天预测结果", x = "", y = "计数") + scale_x_date(date_breaks = "2 weeks", date_labels = "%b %d") + scale_color_tq() + scale_fill_tq() + facet_wrap(~ group, scales = "free_y", ncol = 1) + theme_tq() + theme(axis.text.x = element_text(angle = 45, hjust = 1))
额外说明
- 如果你的数据仅记录有计数的日期,必须补全所有连续日期,否则
tk_ts会将每个观测视为时间序列的一个“步长”,而非对应实际日期。 forecast的h参数仅当时间序列为连续日数据时,才会对应实际天数。
内容的提问来源于stack exchange,提问作者JeniFav
相关产品推荐
相关产品推荐

