Using ProbReco with the fable package

Anastasios Panagiotelis

2023-04-05

The ProbReco package assumes that base probabilistic forecasts are available. This vignette describes how these can be obtained using the fable package. Note that the fable package does currently allow for probabilistic forecast reconciliation, but only under Gaussianity and not using score optimisation.

This vignette considers the case of training reconciliation weights using a rolling window of probabilistic forecasts. A simpler method is to simply use predicted values from a single window of data using the function inscoreopt().

The data

The data sim_hierarchy refer to a simulated 7-variable hierarchy. The bottom level series are all simulated from stationary ARMA models. Noise terms are added so that the residual terms on the bottom levels have higher variance than the middle level residuals, which in turn have higher variance than the top level. For details see [@wp].

library(magrittr)
library(dplyr)
library(tidyr)
library(purrr)
library(fable)
library(ProbReco)
set.seed(1983)
sim_hierarchy
#> # A tibble: 1,506 × 8
#>     Time    Tot     A      B      AA     AB     BA    BB
#>    <dbl>  <dbl> <dbl>  <dbl>   <dbl>  <dbl>  <dbl> <dbl>
#>  1     1  4.75   7.25 -2.50    0.242  7.01   -5.47  2.97
#>  2     2  2.18   2.96 -0.780   6.91  -3.95    3.95 -4.73
#>  3     3 -7.22  -9.18  1.96   -6.92  -2.26   -2.89  4.85
#>  4     4 -0.781  3.95 -4.73    3.26   0.691  -2.46 -2.28
#>  5     5 -4.93   2.31 -7.25    3.08  -0.766  -2.95 -4.29
#>  6     6 -9.30  -5.76 -3.53  -10.1    4.39  -10.8   7.23
#>  7     7 -5.71  -5.47 -0.239   1.65  -7.13    3.80 -4.03
#>  8     8 -5.55  -4.87 -0.686   3.14  -8.00    5.21 -5.90
#>  9     9 -7.26  -4.89 -2.37    4.59  -9.47    6.53 -8.90
#> 10    10 -1.37  -3.41  2.05    1.15  -4.57    4.36 -2.31
#> # … with 1,496 more rows

To ensure the results are stable we have also set a seed.

Set up rolling window

To set up, first we should break the data into a series of rolling windows. This can be done using the map function in the purrr package.

#Length of window
N<-500 

#Make data windows

data_windows<-purrr::map(1:(nrow(sim_hierarchy)-N+1),
                  function(i){return(sim_hierarchy[i:(i+N-1),])})

This creates a list, the first element of which is the data from \(t=1\) to \(t=500\), the second element is the data from \(t=2\) to \(t=501\), etc…

data_windows[[1]]%>%head(3)
#> # A tibble: 3 × 8
#>    Time   Tot     A      B     AA    AB    BA    BB
#>   <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl>
#> 1     1  4.75  7.25 -2.50   0.242  7.01 -5.47  2.97
#> 2     2  2.18  2.96 -0.780  6.91  -3.95  3.95 -4.73
#> 3     3 -7.22 -9.18  1.96  -6.92  -2.26 -2.89  4.85
data_windows[[1]]%>%tail(3)
#> # A tibble: 3 × 8
#>    Time    Tot     A     B    AA     AB      BA     BB
#>   <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>  <dbl>
#> 1   498 -0.776 -1.56 0.787  5.48  -7.04  5.69   -4.91 
#> 2   499  1.14  -1.76 2.90  10.9  -12.6  11.9    -8.98 
#> 3   500  4.10   3.91 0.189  1.80   2.11 -0.0485  0.237
data_windows[[2]]%>%head(3)
#> # A tibble: 3 × 8
#>    Time    Tot     A      B    AA     AB    BA    BB
#>   <dbl>  <dbl> <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl>
#> 1     2  2.18   2.96 -0.780  6.91 -3.95   3.95 -4.73
#> 2     3 -7.22  -9.18  1.96  -6.92 -2.26  -2.89  4.85
#> 3     4 -0.781  3.95 -4.73   3.26  0.691 -2.46 -2.28
data_windows[[2]]%>%tail(3)
#> # A tibble: 3 × 8
#>    Time   Tot     A     B    AA      AB      BA     BB
#>   <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl>  <dbl>
#> 1   499  1.14 -1.76 2.90  10.9  -12.6   11.9    -8.98 
#> 2   500  4.10  3.91 0.189  1.80   2.11  -0.0485  0.237
#> 3   501 11.6   8.70 2.85   8.90  -0.195  6.22   -3.37
data_windows[[500]]%>%head(3)
#> # A tibble: 3 × 8
#>    Time   Tot     A     B     AA     AB      BA     BB
#>   <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>   <dbl>  <dbl>
#> 1   500  4.10  3.91 0.189  1.80   2.11  -0.0485  0.237
#> 2   501 11.6   8.70 2.85   8.90  -0.195  6.22   -3.37 
#> 3   502  6.33  3.86 2.47  -0.515  4.37  -0.911   3.38
data_windows[[500]]%>%tail(3)
#> # A tibble: 3 × 8
#>    Time   Tot     A     B    AA    AB    BA    BB
#>   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1   997  4.57 -2.53  7.11 -2.70 0.169  1.42 5.69 
#> 2   998  7.29 -4.07 11.4  -4.45 0.381  4.71 6.66 
#> 3   999 10.8   7.17  3.63  3.16 4.01   2.69 0.934

Modelling and forecasting

A function for modelling and then obtaining the forecast mean and variance using data from a single window is written below. This can be used with the map family of functions from the purrr package. This function is given as


forecast_window <- function(data_w){
  data_w%>%
    tidyr::pivot_longer(cols = -Time,
               names_to = 'var',
               values_to = 'value')%>%
    as_tsibble(index = 'Time',key='var')%>%
    model(arma11=ARIMA(value~1+pdq(1,0,1)+PDQ(0,0,0)))%>%
    forecast(h=1)%>%
    dplyr::arrange(match(var,c("Tot","A","B","AA","AB","BA","BB")))->f
    mean<-map_dbl(f$.distribution,use_series,mean)
    sd<-map_dbl(f$.distribution,use_series,sd)
 return(list(mean=mean,sd=sd))
}

Using the fable package requires some manipulation of the data. So first, the data frame is converted to long format using pivot_longer. The data must then be converted to a tsibble. Finally, the model function can be used from fable. Here an ARMA(1,1) is fit to each data set. Note that this is only for illustrative purposes, and there will be some misspecification. In practice the order of the ARMA can be chosen automatically. The forecast function is used to obtain the forecast mean and variance for each variable. Only 1 step ahead forecasts are obtained for this example. The variables are arranged in the correct order and then the forecast mean and variance are extracted. Here, dependence is completely ignored, i.e. the base forecasts are indepenedent.

To fit the model and obtain the base forecast mean and variance for each window, the map function can be used. Here we will only fit models to the first 10 windows meaning that t=501 to t=510 will constitute the training data for learning the reconciliation weights.

#Number of windows for training
W<-10
all_fc<-purrr::map(data_windows[1:W],forecast_window)
#> Warning: Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.
#> Unknown or uninitialised column: `.distribution`.

Setting up arguments for reconciliation

The S matrix

The hierarchy has the following \(\boldsymbol{S}\) matrix

S<-matrix(c(1,1,1,1,
            1,1,0,0,
            0,0,1,1,
            1,0,0,0,
            0,1,0,0,
            0,0,1,0,
            0,0,0,1),7,4,byrow = TRUE)

The realisations

The following code obtains the realisations in the form required.


obs_i<-function(i){
  sim_hierarchy%>%
  dplyr::filter(Time==i)%>%
  tidyr::pivot_longer(-Time,names_to = 'var')%>%
  dplyr::arrange(match(var,c("Tot","A","B","AA","AB","BA","BB")))%>%
  dplyr::pull(value)
}

all_y<-purrr::map((N+1):(N+W),obs_i)

This list of length 10 has the vector of true realisations from t=501 as the first element, the vector of true realisations from t=502 as the second element, etc. Note that the arrange and match functions are used to preserve the ordering of the variables from top to bottom. Although any ordering is acceptable, the order must agree with the ordering of the rows in the \(\boldsymbol{S}\) matrix.

The list of probabilistic forecast distributions

The next step is to create a list of functions where the first element generates from the probabilistic forecast distribution at time \(t=501\), the second element generates from the probabilistic forecast distribution at time \(t=502\), etc. To do this write a function that returns a function as follows

make_genfunc<-function(input){
  f<-function(){
    fc_mean<-input$mean
    fc_sd<-input$sd
    out<-matrix(rnorm(350,mean=fc_mean,sd=fc_sd),7,50)
    return(out)
  }
  return(f)
}

The ‘inner’ function f generates fifty, vector-valued, one-step ahead forecasts from a independent Gaussian distributions with mean and standard deviation extracted from input. The ‘outer’ function make_genfunc is required so that the entire list can be created using map

all_prob<-purrr::map(all_fc,make_genfunc)

The object all_prob is in the form required.

Using the ProbReco functions

The total score for bottom up can be found using

  G_bu<-as.vector(rbind(matrix(0,4,4),diag(rep(1,4))))
  es_bu<-total_score(all_y,all_prob,S,G_bu)

The optimal score can now be found using

  opt<-scoreopt(all_y,all_prob,S)

Reference