VAR and ARIMA and their forecast combination with respect to the same variable

[nickcox896]

Suppose we want to forecast the Consumption variable, using two models (VAR and ARIMA). So,

library(fpp3)
h <- 10
training <- us_change |> filter(year(Quarter) <= 2010)

fc1 <-  training |>
  model(
    m1 = VAR(vars(Consumption, Income)))

fc1|>
  forecast(h = h)

 A fable: 10 x 4 [1Q]
# Key:     .model [1]
   .model Quarter .distribution .mean[,"Consumption"] [,"Income"]
   <chr>    <qtr>        <dist>                 <dbl>       <dbl>
 1 m1     2019 Q3        MVN[2]                 0.581       0.651
 2 m1     2019 Q4        MVN[2]                 0.619       0.320
 3 m1     2020 Q1        MVN[2]                 0.768       0.806
 4 m1     2020 Q2        MVN[2]                 0.698       0.753
 5 m1     2020 Q3        MVN[2]                 0.712       0.644
 6 m1     2020 Q4        MVN[2]                 0.768       0.824
 7 m1     2021 Q1        MVN[2]                 0.755       0.751
 8 m1     2021 Q2        MVN[2]                 0.737       0.669
 9 m1     2021 Q3        MVN[2]                 0.748       0.738
10 m1     2021 Q4        MVN[2]                 0.742       0.717


 
fc2 <-  training |>
  model(m2=ARIMA(Consumption))

fc2 |>
  forecast(h = h)

# A fable: 10 x 4 [1Q]
# Key:     .model [1]
   .model Quarter   Consumption .mean
   <chr>    <qtr>        <dist> <dbl>
 1 m3     2019 Q3 N(0.67, 0.33) 0.673
 2 m3     2019 Q4 N(0.64, 0.35) 0.635
 3 m3     2020 Q1 N(0.87, 0.37) 0.867
 4 m3     2020 Q2  N(0.73, 0.4) 0.730
 5 m3     2020 Q3 N(0.73, 0.41) 0.734
 6 m3     2020 Q4 N(0.76, 0.41) 0.756
 7 m3     2021 Q1 N(0.81, 0.41) 0.807
 8 m3     2021 Q2 N(0.71, 0.41) 0.714
 9 m3     2021 Q3 N(0.75, 0.42) 0.752
10 m3     2021 Q4 N(0.74, 0.42) 0.736

The goal is to take the simple average of forecasts from models m1 and m2 (forecast combination) with
respect to the common variable Consumption and compute the measures of accuracy for this combination.

Is that possible with the fable package, using something like

 mutate(combination = ( x1+x2 ) / 2)|>
  forecast(h = h)

[robjhyndman]

library(fpp3)
h <- 10
training <- us_change |> filter(year(Quarter) <= 2010)

var <- training |>
  model(m1 = VAR(vars(Consumption, Income)))
arima <- training |>
  model(m2 = ARIMA(Consumption))

fc <- forecast(arima, h = h) |>
  left_join(
    forecast(var, h = h) |>
      as_tibble() |>
      transmute(Quarter, VAR_fc = mean(.distribution)[, 1]),
    by = "Quarter"
  ) |>
  mutate(ave = (.mean + VAR_fc) / 2)

fc
#> # A fable: 10 x 6 [1Q]
#> # Key:     .model [1]
#>    .model Quarter   Consumption .mean VAR_fc   ave
#>    <chr>    <qtr>        <dist> <dbl>  <dbl> <dbl>
#>  1 m2     2011 Q1 N(0.79, 0.38) 0.791  0.813 0.802
#>  2 m2     2011 Q2  N(0.83, 0.4) 0.827  0.692 0.760
#>  3 m2     2011 Q3 N(0.63, 0.42) 0.631  0.683 0.657
#>  4 m2     2011 Q4 N(0.81, 0.46) 0.814  0.753 0.783
#>  5 m2     2012 Q1 N(0.78, 0.46) 0.780  0.742 0.761
#>  6 m2     2012 Q2 N(0.75, 0.46) 0.747  0.748 0.747
#>  7 m2     2012 Q3 N(0.74, 0.47) 0.745  0.778 0.762
#>  8 m2     2012 Q4 N(0.78, 0.47) 0.779  0.772 0.776
#>  9 m2     2013 Q1 N(0.77, 0.47) 0.767  0.764 0.766
#> 10 m2     2013 Q2 N(0.77, 0.47) 0.767  0.771 0.769

^{Created on 2026-04-15 with reprex v2.1.1}

[nickcox896]

Thank you! To obtain all the foreast measures for the combination I tried

fc1 <-fc  |>
pivot_wider(names_from = .model, values_from = .mean) |>
  pivot_longer(!c(Consumption, Quarter), names_to = ".model", values_to = ".mean") |>
relocate(.model)  %>% select(-Consumption) |> arrange(.model) |>
as_tsibble(index = Quarter, key = .model )

# A tsibble: 30 x 3 [1Q]
# Key:       .model [3]
   .model Quarter .mean
   <chr>    <qtr> <dbl>
 1 VAR_fc 2011 Q1 0.813
 2 VAR_fc 2011 Q2 0.692
 3 VAR_fc 2011 Q3 0.683
 4 VAR_fc 2011 Q4 0.753
 5 VAR_fc 2012 Q1 0.742
 6 VAR_fc 2012 Q2 0.748
 7 VAR_fc 2012 Q3 0.778
 8 VAR_fc 2012 Q4 0.772
 9 VAR_fc 2013 Q1 0.764
10 VAR_fc 2013 Q2 0.771
# ℹ 20 more rows
# ℹ Use `print(n = ...)` to see more rows

fc11<-fc1  |> 
left_join( forecast(arima, h = h) |>  select(.model, Quarter, Consumption), by=c(".model", "Quarter") ) |> 
mutate( Consumption = if_else(is.na( Consumption), .mean,  Consumption))

 A tsibble: 30 x 4 [1Q]
# Key:       .model [3]
   .model Quarter .mean Consumption
   <chr>    <qtr> <dbl>      <dist>
 1 VAR_fc 2011 Q1 0.813   0.8125254
 2 VAR_fc 2011 Q2 0.692   0.6921688
 3 VAR_fc 2011 Q3 0.683   0.6830871
 4 VAR_fc 2011 Q4 0.753    0.752544
 5 VAR_fc 2012 Q1 0.742   0.7419287
 6 VAR_fc 2012 Q2 0.748   0.7478853
 7 VAR_fc 2012 Q3 0.778   0.7784147
 8 VAR_fc 2012 Q4 0.772   0.7718377
 9 VAR_fc 2013 Q1 0.764   0.7639837
10 VAR_fc 2013 Q2 0.771   0.7713127
# ℹ 20 more rows
# ℹ Use `print(n = ...)` to see more rows

 fbl <- as_fable(fc11, response = "Consumption", distribution = "Consumption")
Warning message:
The dimnames of the fable's distribution are missing and have been set to match the response variables. 

fbl %>% forecast:: accuracy(us_change)

# A tibble: 3 × 10
  .model .type     ME  RMSE   MAE   MPE  MAPE  MASE RMSSE   ACF1
  <chr>  <chr>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>
1 VAR_fc Test  -0.415 0.467 0.426 -256.  257. 0.621 0.515 -0.440
2 ave    Test  -0.422 0.473 0.429 -259.  260. 0.625 0.522 -0.503
3 m2     Test  -0.428 0.481 0.432 -261.  262. 0.629 0.531 -0.563

It seems to be working but I am not sure if it is the correct approach.