cross validation and forecast combination with the forecast package

[newsagent]

Using the fable package we ca easily apply the cross validation on the forecast combination. For example

  google_stock <- gafa_stock |>
  filter(Symbol == "GOOG", year(Date) >= 2015) |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)

google_2015 <- google_stock |> filter(year(Date) == 2015)

google_2015_tr <- google_2015 |>
  stretch_tsibble(.init = 180, .step = 1) |>
  relocate(Date, Symbol, .id)
 
fc<-google_2015_tr |>
  model(auto_arima=ARIMA(Close), arima210 = ARIMA(log(Close) ~ pdq(1,1,0))) |>
  mutate(combination = (auto_arima + arima210) / 2) |>
 forecast(h = 8) |>
  group_by(.id) |>
  mutate(h = row_number()) |>
  ungroup() |>
  as_fable(response = "Close", distribution = Close)

fc |> forecast:: accuracy(google_2015, by = c("h", ".model"))

I am trying to translate the above code with the forecast package, using also the theory in these links

accuracy

combinations

So my attempt is this

m1<-  forecast(auto.arima(google_2015$Close), h=8)
m2 <- forecast(Arima(google_2015$Close, order=c(1,1,0), lambda=0,biasadj=TRUE)),  h=8) 
Combination <- (m1[["mean"]] +m2[["mean"]])/2

But then I am not sure how to proceed.

[robjhyndman]

The simplest approach is to use the forecast::tsCV() function. The following will compute the MSE of the three methods.

library(forecast)
library(fpp3)
google_stock <- gafa_stock |>
  filter(Symbol == "GOOG", year(Date) >= 2015) |>
  mutate(day = row_number()) |>
  update_tsibble(index = day, regular = TRUE)
google_2015 <- google_stock |> filter(year(Date) == 2015)

google_tscv_m1 <- tsCV(google_2015$Close, h = 8, function(x, h) {
  auto.arima(x) |>
    forecast(h = h)
})
google_tscv_m2 <- tsCV(google_2015$Close, h = 8, function(x, h) {
  Arima(x, order = c(1, 1, 0), lambda = 0, biasadj = TRUE) |>
    forecast(h = h)
})
google_tscv_comb <- tsCV(google_2015$Close, h = 8, function(x, h) {
  fc1 <- Arima(x, order = c(1, 1, 0), lambda = 0, biasadj = TRUE) |>
    forecast(h = h)
  fc2 <- auto.arima(x) |> forecast(h = h)
  fc2$mean <- 0.5 * (fc1$mean + fc2$mean)
  return(fc2)
})

common_times <- apply(
  !is.na(google_tscv_m1) &
    !is.na(google_tscv_m2) &
    !is.na(google_tscv_comb),
  1,
  all
)

mse <- rbind(
  m1 = colMeans(google_tscv_m1[common_times, ]^2),
  m2 = colMeans(google_tscv_m2[common_times, ]^2),
  comb = colMeans(google_tscv_comb[common_times, ]^2)
)
mse
#>           h=1      h=2      h=3      h=4      h=5      h=6      h=7      h=8
#> m1   138.8834 300.1304 444.5767 575.8839 691.6277 800.2433 918.4026 1022.969
#> m2   133.9775 297.6269 448.3832 581.2425 698.6947 796.2140 909.7168 1005.752
#> comb 135.7658 297.3055 443.5407 573.8943 688.6614 790.0591 904.4725 1003.522

^{Created on 2026-06-02 with reprex v2.1.1}

This is not the most efficient approach, because the forecasts from each model are computed twice, but it is probably the simplest approach. To do it more efficiently, you would need to implement your own variation of tsCV().

[newsagent]

That was very helpful. Thank you!