3.1.1 THE SCENARIO

To help us benchmark multiple models, we will need to split some portion of our data for validation. The strategy here is to cut the last 1 week as our test dataset. Then, we will cut again some (bigger) portion as the train dataset, say, 1 weeks times 8–approximately 2 months. This strategy is just the simplified version of Rolling Forecasting Origin, with only having one pair of train and test sample.

So the first step here is to get the start and end date of the train and test sample. The most straighforward way is to define the train and test size, then recursively get the start and end index for each:

# train-val-test size
test_size <- 24 * 7
train_size <- 24 * 7 * 8

# get the min-max of the time index for each sample
test_end <- max(subarea_train$datetime)
test_start <- test_end - hours(test_size) + hours(1)

train_end <- test_start - hours(1)
train_start <- train_end - hours(train_size) + hours(1)

We combine the start and end indices into a date interval:

# get the interval of each samples
intrain <- interval(train_start, train_end)
intest <- interval(test_start, test_end)

intrain

## [1] 2017-10-01 UTC--2017-11-25 23:00:00 UTC

#> [1] 2017-10-01 UTC--2017-11-25 23:00:00 UTC

intest

## [1] 2017-11-26 UTC--2017-12-02 23:00:00 UTC

#> [1] 2017-11-26 UTC--2017-12-02 23:00:00 UTC

Here’s the data after splitting process: data train & test

# plot the train and test
p <- subarea_train %>%
  mutate(sample = case_when(
    datetime %within% intrain ~ "train",
    datetime %within% intest ~ "test"
  )) %>%
  drop_na() %>%
  mutate(sample = factor(sample, levels = c("train", "test"))) %>%
  ggplot(aes(x = datetime, y = demand, colour = sample)) +
    geom_line() +
    labs(x = NULL, y = NULL, colour = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    tidyquant::theme_tq() +
    tidyquant::scale_colour_tq()

ggplotly(p)

3.1.2 DATA PREPROCESS USING RECIPES

Data preprocessing is a very crucial step in time series model fitting. In this Capstone project, we will use recipes package for data preprocessing and comparing multiple preprocess specifications.

Since recipe package work columnwise, we need to convert our data into a wide format first:

# convert to wide format
subarea_wide <- subarea_train %>% 
  spread(src_sub_area, demand)

subarea_wide

## # A tibble: 1,512 x 4
##    datetime            sxk97 sxk9e sxk9s
##    <dttm>              <dbl> <dbl> <dbl>
##  1 2017-10-01 00:00:00     6     8    10
##  2 2017-10-01 01:00:00     4     2     3
##  3 2017-10-01 02:00:00     9     3     1
##  4 2017-10-01 03:00:00     2     2     0
##  5 2017-10-01 04:00:00     5     1     0
##  6 2017-10-01 05:00:00     4     2     0
##  7 2017-10-01 06:00:00     1     1     0
##  8 2017-10-01 07:00:00     0     0     1
##  9 2017-10-01 08:00:00     3     2     5
## 10 2017-10-01 09:00:00     3    11     5
## # ... with 1,502 more rows

Then we could start to define the preprocess recipe(), and bake() our data based on the defined recipe:

# recipe is a set of steps, obj: make the data additive
# we only do this with data train
rec <- recipe(~ ., filter(subarea_wide, datetime %within% intrain)) %>%
  step_sqrt(all_numeric()) %>%
  step_center(all_numeric()) %>%
  step_scale(all_numeric()) %>%
  prep()

# preview the bake results
subarea_wide <- bake(rec, subarea_wide)

subarea_wide

## # A tibble: 1,512 x 4
##    datetime             sxk97  sxk9e  sxk9s
##    <dttm>               <dbl>  <dbl>  <dbl>
##  1 2017-10-01 00:00:00 -0.594 -0.697 -0.111
##  2 2017-10-01 01:00:00 -0.841 -1.28  -0.780
##  3 2017-10-01 02:00:00 -0.291 -1.15  -1.12 
##  4 2017-10-01 03:00:00 -1.16  -1.28  -1.59 
##  5 2017-10-01 04:00:00 -0.711 -1.44  -1.59 
##  6 2017-10-01 05:00:00 -0.841 -1.28  -1.59 
##  7 2017-10-01 06:00:00 -1.39  -1.44  -1.59 
##  8 2017-10-01 07:00:00 -1.94  -1.85  -1.12 
##  9 2017-10-01 08:00:00 -0.989 -1.28  -0.544
## 10 2017-10-01 09:00:00 -0.989 -0.497 -0.544
## # ... with 1,502 more rows

If we use recipes package, the next steps is to create a revert back function:

# revert back function
rec_revert <- function(vector, rec, varname) {

  # store recipe values
  rec_center <- rec$steps[[2]]$means[varname]
  rec_scale <- rec$steps[[3]]$sds[varname]

  # convert back based on the recipe
  results <- (vector * rec_scale + rec_center) ^ 2

  # add additional adjustment if necessary
  results <- round(results)

  # return the results
  results
}

This revert back function will convert back the data to its original form.

Now we can convert our data into the long format again:

# convert back to long format
subarea_long <- subarea_wide %>% 
  gather(src_sub_area, demand, -datetime)
  
subarea_long

## # A tibble: 4,536 x 3
##    datetime            src_sub_area demand
##    <dttm>              <chr>         <dbl>
##  1 2017-10-01 00:00:00 sxk97        -0.594
##  2 2017-10-01 01:00:00 sxk97        -0.841
##  3 2017-10-01 02:00:00 sxk97        -0.291
##  4 2017-10-01 03:00:00 sxk97        -1.16 
##  5 2017-10-01 04:00:00 sxk97        -0.711
##  6 2017-10-01 05:00:00 sxk97        -0.841
##  7 2017-10-01 06:00:00 sxk97        -1.39 
##  8 2017-10-01 07:00:00 sxk97        -1.94 
##  9 2017-10-01 08:00:00 sxk97        -0.989
## 10 2017-10-01 09:00:00 sxk97        -0.989
## # ... with 4,526 more rows

3.1.3 NESTED MODEL FITTING AND FORECASTING

In functional programming using purrr, we need to convert our tbl into a nested tbl. Nested data is like a table inside a table, which could be controlled using a key indicator; in other words, we can have tbl for each subarea and samples. Using this format, the fitting and forecasting process would be very versatile, yet we can still convert the results as long as we have a proper key like provider.

Let’s start by converting our tbl into a nested tbl. First, we need to add sample indicator so it could be recognized as a key when we nest() the data:

# adjust by sample
subarea_nest <- subarea_long %>%
  mutate(sample = case_when(
    datetime %within% intrain ~ "train",
    datetime %within% intest ~ "test"
  )) %>%
  drop_na()

# nest the data train
subarea_train_nest <- subarea_nest %>% 
  group_by(src_sub_area, sample) %>%
  nest(.key = "data") %>%
  spread(sample, data)

subarea_train_nest

## # A tibble: 3 x 3
##   src_sub_area test               train               
##   <chr>        <list>             <list>              
## 1 sxk97        <tibble [168 x 2]> <tibble [1,344 x 2]>
## 2 sxk9e        <tibble [168 x 2]> <tibble [1,344 x 2]>
## 3 sxk9s        <tibble [168 x 2]> <tibble [1,344 x 2]>

# nest the data submission too
submit_nest <- data_submission %>% nest(datetime)
submit_nest 

## # A tibble: 3 x 3
##   src_sub_area demand data              
##   <chr>        <lgl>  <list>            
## 1 sxk97        NA     <tibble [168 x 1]>
## 2 sxk9e        NA     <tibble [168 x 1]>
## 3 sxk9s        NA     <tibble [168 x 1]>

3.1.4 PREPARING THE DATA MODEL LIST

For data and model combination, we could start by defining some options for data representation. Recall that our series have a relatively high frequency, so we could consider two option of data representation here:
- a ts object with daily seasonality (24 hours), and
- an msts with daily and weekly seasonality (24 * 7).

To incorporate them into our nested data, we need to create another nested data frame containing the data daily and weekly, after that we will do a function for converting the data into the specified data. In this case, we would like to compare multiple seasonality specifications.

Let’s start with making a named list containing the transformation functions:

# function: data funs list
data_funs <- list(
  ts = function(x) ts(x$demand, frequency = 24),
  msts = function(x) msts(x$demand, seasonal.periods = c(24, 24 * 7))
)

data_funs

## $ts
## function (x) 
## ts(x$demand, frequency = 24)
## 
## $msts
## function (x) 
## msts(x$demand, seasonal.periods = c(24, 24 * 7))

Then we convert the list into a tbl using enframe().
Note that we should also give a key – which is the src_sub_area in our case – so we could use left_join() later.

Use rep() function:

# convert to nested
data_funs <- data_funs %<>%
  rep(length(unique(subarea_nest$src_sub_area))) %>%
  enframe("data_fun_name", "data_fun") %>%
  mutate(src_sub_area =
    sort(rep(unique(subarea_nest$src_sub_area), length(unique(.$data_fun_name))))
  )

data_funs

## # A tibble: 6 x 3
##   data_fun_name data_fun src_sub_area
##   <chr>         <list>   <chr>       
## 1 ts            <fn>     sxk97       
## 2 msts          <fn>     sxk97       
## 3 ts            <fn>     sxk9e       
## 4 msts          <fn>     sxk9e       
## 5 ts            <fn>     sxk9s       
## 6 msts          <fn>     sxk9s

Then the last steps here is to join the nested function with our Nested Model Fitting:

# combine with models
subarea_nest_fit <- subarea_train_nest %>% 
  left_join(data_funs)

## Joining, by = "src_sub_area"

subarea_nest_fit

## # A tibble: 6 x 5
##   src_sub_area test               train              data_fun_name data_fun
##   <chr>        <list>             <list>             <chr>         <list>  
## 1 sxk97        <tibble [168 x 2]> <tibble [1,344 x ~ ts            <fn>    
## 2 sxk97        <tibble [168 x 2]> <tibble [1,344 x ~ msts          <fn>    
## 3 sxk9e        <tibble [168 x 2]> <tibble [1,344 x ~ ts            <fn>    
## 4 sxk9e        <tibble [168 x 2]> <tibble [1,344 x ~ msts          <fn>    
## 5 sxk9s        <tibble [168 x 2]> <tibble [1,344 x ~ ts            <fn>    
## 6 sxk9s        <tibble [168 x 2]> <tibble [1,344 x ~ msts          <fn>

3.1.5 PREPARING THE TIME SERIES MODEL LIST

Similar to when we create the data representation (daily and weekly) list, we also make some time series models as a nested list.

Again, let’s start by making a list of models. We would like to compare multiple model specifications. For this Capstone Project, we will use auto.arima(), ets(), stlm(), and tbats(). We need to make some functions to call those model, and store them inside a list:

# time series models list
models <- list(
  auto.arima = function(x) auto.arima(x),
  ets = function(x) ets(x),
  stlm = function(x) stlm(x),
  tbats = function(x) tbats(x, use.box.cox = FALSE)
)

models

## $auto.arima
## function (x) 
## auto.arima(x)
## 
## $ets
## function (x) 
## ets(x)
## 
## $stlm
## function (x) 
## stlm(x)
## 
## $tbats
## function (x) 
## tbats(x, use.box.cox = FALSE)

Then we can convert it into a nested format like previous example:

# convert to nested
models <- models %>%
  rep(length(unique(subarea_nest$src_sub_area))) %>%
  enframe("model_name", "model") %>%
  mutate(src_sub_area =
    sort(rep(unique(subarea_nest$src_sub_area), length(unique(.$model_name))))
  )

models

## # A tibble: 12 x 3
##    model_name model  src_sub_area
##    <chr>      <list> <chr>       
##  1 auto.arima <fn>   sxk97       
##  2 ets        <fn>   sxk97       
##  3 stlm       <fn>   sxk97       
##  4 tbats      <fn>   sxk97       
##  5 auto.arima <fn>   sxk9e       
##  6 ets        <fn>   sxk9e       
##  7 stlm       <fn>   sxk9e       
##  8 tbats      <fn>   sxk9e       
##  9 auto.arima <fn>   sxk9s       
## 10 ets        <fn>   sxk9s       
## 11 stlm       <fn>   sxk9s       
## 12 tbats      <fn>   sxk9s

And finally, we can join the result into our nested data. We also apply some rule here:
since ets() and auto.arima() are not suitable for data with multiple seasonality time series, we use filter to remove them out:

# combine with models
subarea_nest_model <- subarea_nest_fit %>%
  left_join(models) %>%
  filter(
    !(model_name == "ets" & data_fun_name == "msts"),
    !(model_name == "auto.arima" & data_fun_name == "msts")
  )

## Joining, by = "src_sub_area"

subarea_nest_model

## # A tibble: 18 x 7
##    src_sub_area test      train     data_fun_name data_fun model_name model
##    <chr>        <list>    <list>    <chr>         <list>   <chr>      <lis>
##  1 sxk97        <tibble ~ <tibble ~ ts            <fn>     auto.arima <fn> 
##  2 sxk97        <tibble ~ <tibble ~ ts            <fn>     ets        <fn> 
##  3 sxk97        <tibble ~ <tibble ~ ts            <fn>     stlm       <fn> 
##  4 sxk97        <tibble ~ <tibble ~ ts            <fn>     tbats      <fn> 
##  5 sxk97        <tibble ~ <tibble ~ msts          <fn>     stlm       <fn> 
##  6 sxk97        <tibble ~ <tibble ~ msts          <fn>     tbats      <fn> 
##  7 sxk9e        <tibble ~ <tibble ~ ts            <fn>     auto.arima <fn> 
##  8 sxk9e        <tibble ~ <tibble ~ ts            <fn>     ets        <fn> 
##  9 sxk9e        <tibble ~ <tibble ~ ts            <fn>     stlm       <fn> 
## 10 sxk9e        <tibble ~ <tibble ~ ts            <fn>     tbats      <fn> 
## 11 sxk9e        <tibble ~ <tibble ~ msts          <fn>     stlm       <fn> 
## 12 sxk9e        <tibble ~ <tibble ~ msts          <fn>     tbats      <fn> 
## 13 sxk9s        <tibble ~ <tibble ~ ts            <fn>     auto.arima <fn> 
## 14 sxk9s        <tibble ~ <tibble ~ ts            <fn>     ets        <fn> 
## 15 sxk9s        <tibble ~ <tibble ~ ts            <fn>     stlm       <fn> 
## 16 sxk9s        <tibble ~ <tibble ~ ts            <fn>     tbats      <fn> 
## 17 sxk9s        <tibble ~ <tibble ~ msts          <fn>     stlm       <fn> 
## 18 sxk9s        <tibble ~ <tibble ~ msts          <fn>     tbats      <fn>

3.1.6 EXECUTE THE NESTED FITTING

To execute the model fitting, we need to wrap up the needed arguments as a list using map() function. Then, we could call the function using invoke_map().

We need to do this for data transformation using the function inside data_fun, then continue to fit the model with the same process using the function inside model.

Note: this function will take some time, so we already save it into RDS file, we can only readRDS and get the result faster.

# invoke nested fitting
# subarea_nest_model <- subarea_nest_model  %>% 
#   mutate(
#     params = map(train, ~ list(x = .x)),
#     data = invoke_map(data_fun, params),
#     params = map(data, ~ list(x = .x)),
#     fitted = invoke_map(model, params)
#     # error = map_dbl(fitted, mae_vec(.x$demand, .x$estimate))
#   ) %>%
#   select(-data, -params)
# 
# saveRDS(subarea_nest_model, "subarea_nest_model.RDS")

subarea_nest_model <- readRDS("subarea_nest_model.RDS")
subarea_nest_model

## # A tibble: 18 x 8
##    src_sub_area test   train data_fun_name data_fun model_name model fitted
##    <chr>        <list> <lis> <chr>         <list>   <chr>      <lis> <list>
##  1 sxk97        <tibb~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
##  2 sxk97        <tibb~ <tib~ ts            <fn>     ets        <fn>  <ets> 
##  3 sxk97        <tibb~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
##  4 sxk97        <tibb~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
##  5 sxk97        <tibb~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
##  6 sxk97        <tibb~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
##  7 sxk9e        <tibb~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
##  8 sxk9e        <tibb~ <tib~ ts            <fn>     ets        <fn>  <ets> 
##  9 sxk9e        <tibb~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
## 10 sxk9e        <tibb~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
## 11 sxk9e        <tibb~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
## 12 sxk9e        <tibb~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
## 13 sxk9s        <tibb~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
## 14 sxk9s        <tibb~ <tib~ ts            <fn>     ets        <fn>  <ets> 
## 15 sxk9s        <tibb~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
## 16 sxk9s        <tibb~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
## 17 sxk9s        <tibb~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
## 18 sxk9s        <tibb~ <tib~ msts          <fn>     tbats      <fn>  <tbat~

# calculate test errors
subarea_error <- subarea_nest_model %>%
  mutate(error =
    map(fitted, ~ forecast(.x, h = 24 * 7)) %>%
    map2_dbl(test, ~ mae_vec(truth = .y$demand, estimate = .x$mean))
  ) %>%
  arrange(src_sub_area, error) #arrange based on smallest error for each area

subarea_error

## # A tibble: 18 x 9
##    src_sub_area test  train data_fun_name data_fun model_name model fitted
##    <chr>        <lis> <lis> <chr>         <list>   <chr>      <lis> <list>
##  1 sxk97        <tib~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
##  2 sxk97        <tib~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
##  3 sxk97        <tib~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
##  4 sxk97        <tib~ <tib~ ts            <fn>     ets        <fn>  <ets> 
##  5 sxk97        <tib~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
##  6 sxk97        <tib~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
##  7 sxk9e        <tib~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
##  8 sxk9e        <tib~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
##  9 sxk9e        <tib~ <tib~ ts            <fn>     ets        <fn>  <ets> 
## 10 sxk9e        <tib~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
## 11 sxk9e        <tib~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
## 12 sxk9e        <tib~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
## 13 sxk9s        <tib~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
## 14 sxk9s        <tib~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
## 15 sxk9s        <tib~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
## 16 sxk9s        <tib~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
## 17 sxk9s        <tib~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
## 18 sxk9s        <tib~ <tib~ ts            <fn>     ets        <fn>  <ets> 
## # ... with 1 more variable: error <dbl>

# show the area, data fun name, model & error
subarea_error %>%
  select(src_sub_area, ends_with("_name"), error)

## # A tibble: 18 x 4
##    src_sub_area data_fun_name model_name error
##    <chr>        <chr>         <chr>      <dbl>
##  1 sxk97        msts          tbats      0.488
##  2 sxk97        msts          stlm       0.543
##  3 sxk97        ts            tbats      0.573
##  4 sxk97        ts            ets        0.574
##  5 sxk97        ts            stlm       0.600
##  6 sxk97        ts            auto.arima 0.702
##  7 sxk9e        msts          tbats      0.359
##  8 sxk9e        msts          stlm       0.366
##  9 sxk9e        ts            ets        0.419
## 10 sxk9e        ts            tbats      0.422
## 11 sxk9e        ts            stlm       0.438
## 12 sxk9e        ts            auto.arima 0.502
## 13 sxk9s        msts          tbats      0.421
## 14 sxk9s        msts          stlm       0.428
## 15 sxk9s        ts            tbats      0.431
## 16 sxk9s        ts            auto.arima 0.503
## 17 sxk9s        ts            stlm       0.540
## 18 sxk9s        ts            ets        0.550

Get the best model by choosing the smallest error:

subarea_best_model <- subarea_error %>%
  group_by(src_sub_area) %>% 
  arrange(error) %>% 
  slice(1) %>%  #choose the smallest error & best model for each area directly
  ungroup() %>% 
  select(src_sub_area, ends_with("_name"),error)

Do left join in order to get the whole information regarding the Best model:

subarea_best_model <- subarea_best_model %>% 
  select(-error) %>% 
  left_join(subarea_error)%>% 
  select(-error)

## Joining, by = c("src_sub_area", "data_fun_name", "model_name")

subarea_best_model

## # A tibble: 3 x 8
##   src_sub_area data_fun_name model_name test   train  data_fun model fitted
##   <chr>        <chr>         <chr>      <list> <list> <list>   <lis> <list>
## 1 sxk97        msts          tbats      <tibb~ <tibb~ <fn>     <fn>  <tbat~
## 2 sxk9e        msts          tbats      <tibb~ <tibb~ <fn>     <fn>  <tbat~
## 3 sxk9s        msts          tbats      <tibb~ <tibb~ <fn>     <fn>  <tbat~

3.1.7 UNNESTING THE RESULT

Beside measuring the error, we can also compare the forecast results to the real test series through graphical analysis. But to do that, we need to make a tbl containing our forecast to the test dataset, then do some spread-gather tricks to make a set of keys that unique for each data representations, models, and also one for the forecast itself. If we get to that format, we could conveniently unnest() the data into a proper long format

Let’s start with creating the forecast:

subarea_test <- subarea_error %>% 
  mutate(
    forecast =
      map(fitted, ~ forecast(.x, h = 24 * 7)) %>%
      map2(test, ~ tibble(
        datetime = .y$datetime,
        demand = as.vector(.x$mean)
      )),
    key = paste(data_fun_name, model_name, sep = "-")
  )
  
subarea_test

## # A tibble: 18 x 11
##    src_sub_area test  train data_fun_name data_fun model_name model fitted
##    <chr>        <lis> <lis> <chr>         <list>   <chr>      <lis> <list>
##  1 sxk97        <tib~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
##  2 sxk97        <tib~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
##  3 sxk97        <tib~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
##  4 sxk97        <tib~ <tib~ ts            <fn>     ets        <fn>  <ets> 
##  5 sxk97        <tib~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
##  6 sxk97        <tib~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
##  7 sxk9e        <tib~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
##  8 sxk9e        <tib~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
##  9 sxk9e        <tib~ <tib~ ts            <fn>     ets        <fn>  <ets> 
## 10 sxk9e        <tib~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
## 11 sxk9e        <tib~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
## 12 sxk9e        <tib~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
## 13 sxk9s        <tib~ <tib~ msts          <fn>     tbats      <fn>  <tbat~
## 14 sxk9s        <tib~ <tib~ msts          <fn>     stlm       <fn>  <stlm>
## 15 sxk9s        <tib~ <tib~ ts            <fn>     tbats      <fn>  <tbat~
## 16 sxk9s        <tib~ <tib~ ts            <fn>     auto.arima <fn>  <ARIM~
## 17 sxk9s        <tib~ <tib~ ts            <fn>     stlm       <fn>  <stlm>
## 18 sxk9s        <tib~ <tib~ ts            <fn>     ets        <fn>  <ets> 
## # ... with 3 more variables: error <dbl>, forecast <list>, key <chr>

Then do some spread-gather to create a proper key:

subarea_test_key <- subarea_test %>%
  select(src_sub_area, key, actual = test, forecast) %>%
  spread(key, forecast) %>%
  gather(key, value, -src_sub_area)

subarea_test_key

## # A tibble: 21 x 3
##    src_sub_area key           value             
##    <chr>        <chr>         <list>            
##  1 sxk97        actual        <tibble [168 x 2]>
##  2 sxk9e        actual        <tibble [168 x 2]>
##  3 sxk9s        actual        <tibble [168 x 2]>
##  4 sxk97        msts-stlm     <tibble [168 x 2]>
##  5 sxk9e        msts-stlm     <tibble [168 x 2]>
##  6 sxk9s        msts-stlm     <tibble [168 x 2]>
##  7 sxk97        msts-tbats    <tibble [168 x 2]>
##  8 sxk9e        msts-tbats    <tibble [168 x 2]>
##  9 sxk9s        msts-tbats    <tibble [168 x 2]>
## 10 sxk97        ts-auto.arima <tibble [168 x 2]>
## # ... with 11 more rows

The last but not least, unnest() the series data, and apply the revert back function:

subarea_test_unnest <- subarea_test_key %>%
  unnest(value) %>%
  mutate(demand = rec_revert(demand, rec, src_sub_area))
  
subarea_test_unnest

## # A tibble: 3,528 x 4
##    src_sub_area key    datetime            demand
##    <chr>        <chr>  <dttm>               <dbl>
##  1 sxk97        actual 2017-11-26 00:00:00     38
##  2 sxk97        actual 2017-11-26 01:00:00     21
##  3 sxk97        actual 2017-11-26 02:00:00     10
##  4 sxk97        actual 2017-11-26 03:00:00     10
##  5 sxk97        actual 2017-11-26 04:00:00      5
##  6 sxk97        actual 2017-11-26 05:00:00      2
##  7 sxk97        actual 2017-11-26 06:00:00      0
##  8 sxk97        actual 2017-11-26 07:00:00     11
##  9 sxk97        actual 2017-11-26 08:00:00      4
## 10 sxk97        actual 2017-11-26 09:00:00      4
## # ... with 3,518 more rows

With the resulting tbl, we can compare the forecast and actual data on test like this:

# plot forecast on test
p <- subarea_test_unnest %>%
  ggplot(aes(x = datetime, y = demand, colour = key)) +
    geom_line() +
    labs(x = NULL, y = NULL, colour = NULL) +
    facet_wrap(~ src_sub_area, scale = "free", ncol = 1) +
    tidyquant::theme_tq() +
    tidyquant::scale_colour_tq()

ggplotly(p)