调整 R 中多个变量的预测代码

Question

这个问题派生自这个问题 在摘要中按组连接残差 R 中的预测表其中使用 ETS 函数对每组（1 和 2）进行预测唯一一个严重的区别是它适用于一个变量。 如果我有很多变量，我必须一次对所有变量执行预测。让我们举个例子

df=structure(list(Variable = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L), .Label = c("x", "y"), class = "factor"), group = c(1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 2L, 2L), year = c(1973L, 1974L, 1975L, 1976L, 
1977L, 1978L, 1973L, 1974L, 1975L, 1976L, 1977L, 1978L, 1973L, 
1974L, 1975L, 1976L, 1977L, 1978L, 1973L, 1974L, 1975L, 1976L, 
1977L, 1978L), Jan = c(9007L, 7750L, 8162L, 7717L, 7792L, 7836L, 
9007L, 7750L, 8162L, 7717L, 7792L, 7836L, 9007L, 7750L, 8162L, 
7717L, 7792L, 7836L, 9007L, 7750L, 8162L, 7717L, 7792L, 7836L
), Feb = c(8106L, 6981L, 7306L, 7461L, 6957L, 6892L, 8106L, 6981L, 
7306L, 7461L, 6957L, 6892L, 8106L, 6981L, 7306L, 7461L, 6957L, 
6892L, 8106L, 6981L, 7306L, 7461L, 6957L, 6892L), Mar = c(8928L, 
8038L, 8124L, 7767L, 7726L, 7791L, 8928L, 8038L, 8124L, 7767L, 
7726L, 7791L, 8928L, 8038L, 8124L, 7767L, 7726L, 7791L, 8928L, 
8038L, 8124L, 7767L, 7726L, 7791L), Apr = c(9137L, 8422L, 7870L, 
7925L, 8106L, 8192L, 9137L, 8422L, 7870L, 7925L, 8106L, 8192L, 
9137L, 8422L, 7870L, 7925L, 8106L, 8192L, 9137L, 8422L, 7870L, 
7925L, 8106L, 8192L), May = c(10017L, 8714L, 9387L, 8623L, 8890L, 
9115L, 10017L, 8714L, 9387L, 8623L, 8890L, 9115L, 10017L, 8714L, 
9387L, 8623L, 8890L, 9115L, 10017L, 8714L, 9387L, 8623L, 8890L, 
9115L), Jun = c(10826L, 9512L, 9556L, 8945L, 9299L, 9434L, 10826L, 
9512L, 9556L, 8945L, 9299L, 9434L, 10826L, 9512L, 9556L, 8945L, 
9299L, 9434L, 10826L, 9512L, 9556L, 8945L, 9299L, 9434L), Jul = c(11317L, 
10120L, 10093L, 10078L, 10625L, 10484L, 11317L, 10120L, 10093L, 
10078L, 10625L, 10484L, 11317L, 10120L, 10093L, 10078L, 10625L, 
10484L, 11317L, 10120L, 10093L, 10078L, 10625L, 10484L), Aug = c(10744L, 
9823L, 9620L, 9179L, 9302L, 9827L, 10744L, 9823L, 9620L, 9179L, 
9302L, 9827L, 10744L, 9823L, 9620L, 9179L, 9302L, 9827L, 10744L, 
9823L, 9620L, 9179L, 9302L, 9827L), Sep = c(9713L, 8743L, 8285L, 
8037L, 8314L, 9110L, 9713L, 8743L, 8285L, 8037L, 8314L, 9110L, 
9713L, 8743L, 8285L, 8037L, 8314L, 9110L, 9713L, 8743L, 8285L, 
8037L, 8314L, 9110L), Oct = c(9938L, 9129L, 8466L, 8488L, 8850L, 
9070L, 9938L, 9129L, 8466L, 8488L, 8850L, 9070L, 9938L, 9129L, 
8466L, 8488L, 8850L, 9070L, 9938L, 9129L, 8466L, 8488L, 8850L, 
9070L), Nov = c(9161L, 8710L, 8160L, 7874L, 8265L, 8633L, 9161L, 
8710L, 8160L, 7874L, 8265L, 8633L, 9161L, 8710L, 8160L, 7874L, 
8265L, 8633L, 9161L, 8710L, 8160L, 7874L, 8265L, 8633L), Dec = c(8927L, 
8680L, 8034L, 8647L, 8796L, 9240L, 8927L, 8680L, 8034L, 8647L, 
8796L, 9240L, 8927L, 8680L, 8034L, 8647L, 8796L, 9240L, 8927L, 
8680L, 8034L, 8647L, 8796L, 9240L)), .Names = c("Variable", "group", 
"year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", 
"Sep", "Oct", "Nov", "Dec"), class = "data.frame", row.names = c(NA, 
-24L))

变量

列具有变量 x 和 y，每个变量都有组 1 和 2 以执行预测

加载库

load_pkgs <- c("forecast", "zoo", "timetk", "tidyverse") 
sapply(load_pkgs, function(x) suppressPackageStartupMessages(library(x, character.only = T)))

第 1 步：预处理

# perform split by group
ld <- split(df[, -1], df$group)

 # Tidy-up the splits
library(lubridate)
ld <- lapply(ld, function(x) {
  x %>%
    gather(key, value, -year) %>%
    unite(date, year, key, sep = "-") %>%
    mutate(date = paste0(date, "-01")) %>%
    mutate(date =ymd(date))    
})

然后

# Transform time series to ts objects
ld <- lapply(ld, function(x) {
  yr <- lubridate::year(min(x$date))
  mth <- lubridate::month(min(x$date))
  timetk::tk_ts(data = x, select = value, frequency = 12,
                start = c(yr, mth))
})

第 2 步：使用 ets 进行训练和预测

# helping function
make_df <- function(ts_obj) {
  ts_df <- timetk::tk_tbl(preserve_index = TRUE, ts_obj) %>%
    mutate(index = zoo::as.Date(x = .$index, frac = 0)) %>% 
    dplyr::rename(date = index)
  return(ts_df)
}

以下函数训练 ets 并预测未来 12 个月;然后，它使用拟合值和预测值准备表：

lts <- lapply(ld, function(ts_obj) {
# train ets model and get fitted results
res_model <- ets(ts_obj, model = "ZZZ")
res_fit <- ts(as.numeric(res_model$fitted), start = start(ts_obj), frequency = 12)
# add extra metrics you may be interested in
model <- res_model[["method"]]
mse <- res_model[["mse"]]
# get forecasts for the next 12 months
res_fct <- forecast(res_model, h = 12)
res_fcst <- ts(res_fct$mean, start = end(ts_obj) + 1/12, frequency = 12)
# transform results to tbl
# for fitted output we keep the residuals and the 95% CI
res_fit_tbl <- make_df(res_fit) %>%
  mutate(residuals = as.numeric(res_model[["residuals"]])) %>%
  mutate(CI95_upper = value + 1.96*sqrt(res_model$sigma2), 
         CI95_lower = value - 1.96*sqrt(res_model$sigma2))
# the forecast output does not have residuals
res_fcst_tbl <- make_df(res_fcst)
return(list(res_fit_tbl = res_fit_tbl, res_fcst_tbl = res_fcst_tbl, model = model, mse = mse)) # don't forget to pass the extra metrics as output
})

步骤 3：将不同组的拟合输出和预测输出汇集在一起

# add groups back + other metrics of interest
lts_all <- lapply(names(lts), function(x, lts) {
  output_fit <- lts[[x]][["res_fit_tbl"]] %>%
    mutate(group = x,
           model = lts[[x]][["model"]],
           mse = lts[[x]][["mse"]])
  output_fcst <- lts[[x]][["res_fcst_tbl"]] %>%
    mutate(group = x)
  return(list(output_fit=output_fit, output_fcst=output_fcst))
  }, lts)

然后

# bring together the fitted respectively forecasting results
output_fit_all <- lapply(lts_all, function(x) x[[1]])
output_fit_all <- bind_rows(output_fit_all)
output_fcst_all <- lapply(lts_all, function(x) x[[2]])
output_fcst_all <- bind_rows(output_fcst_all)

如何做，此代码对所有变量执行预测，例如可重现的示例

Answer 1

您只需要从"步骤 1"开头的Variable/group对中获取一个新组，代码应该可以工作：

# load libraries
load_pkgs <- c("forecast", "zoo", "timetk", "tidyverse", "lubridate") 
sapply(load_pkgs, function(x) suppressPackageStartupMessages(library(x, character.only = T)))
# get new group from pair variable/group
df <- df %>%
  unite_("group", c("Variable", "group"))

接下来，运行步骤 1 - 步骤 3 中的代码，在步骤 3 中进行小更新：将拟合的预测结果放在一起时，为了提取Variable和group列，您需要更新步骤 3 的第二部分，如下所示：

# bring together the fitted respectively forecasting results
output_fit_all <- lapply(lts_all, function(x) x[[1]])
output_fit_all <- bind_rows(output_fit_all)  %>%
  separate(group, c("Variable", "group"))
output_fcst_all <- lapply(lts_all, function(x) x[[2]])
output_fcst_all <- bind_rows(output_fcst_all)   %>%
  separate(group, c("Variable", "group"))