Time Series Forecasting Lab (Part 2) - Feature Engineering with Recipes

Time Series Forecasting Lab (Part 2) - Feature Engineering with Recipes

Cover photo by Kristine Tumanyan on Unsplash

Go to R-bloggers for R news and tutorials contributed by hundreds of R bloggers.

This is the second of a series of 6 articles about time series forecasting with panel data and ensemble stacking with R. Click here for Part 1.

Through these articles I will be putting into practice what I have learned from the Business Science University training course 1 DS4B 203-R: High-Performance Time Series Forecasting", delivered by Matt Dancho. If you are looking to gain a high level of expertise in time series with R I strongly recommend this course.

The objective of this article is to explain how to perform time series preprocessing and feature enginneering with timetk . timetk extends recipes for time series.

I will explain what can and cannot be done with recipes for time series panel data.

Introduction

Since the ultimate goal of this series of articles is to perform time series forecast with machine learning, we all know that maximing model performance requires a lot of exprimentation with different feature engineered datasets.

Feature engineering is the most critical part of time series analysis and with recipes you can "use dplyr-like pipeable sequences of feature engineering steps to get your data ready for modeling".

Prerequisites

I assume you are already familiar with the following topics, packages and terms:

  • time series data format: timestamp, measure, key(s), features

  • time series lags, rolling lags, Fourier terms

  • dplyr or tidyverse R packages

  • normalization and standardization

  • dummy variables (one-hot encoding)

  • dataset split into training and test for modeling

Packages

The following packages must be loaded:

# PART 1 - FEATURE ENGINEERING
library(tidyverse)  # loading dplyr, tibble, ggplot2, .. dependencies
library(timetk)  # using timetk plotting, diagnostics and augment operations
library(tsibble)  # for month to Date conversion
library(tsibbledata)  # for aus_retail dataset
library(fastDummies)  # for dummyfying categorical variables
library(skimr) # for quick statistics

# PART 2 -  FEATURE ENGINEERING WITH RECIPES
library(tidymodels)

End-to-end FE process with recipe

Below is a summary diagram of the end-to-end FE process that will be implemented in this article. Notice that the recipe creation will be the last step of the process. The whole process is described hereafter:

  1. Perform measure transformations (here, Turnover) .

  2. Add a future frame with prediction length

  3. Add Fourier, lag, and rolling lag features. Add a rowid column.

  4. Perform a first split into prepared and future datasets.

  5. Apply the training and test split on the prepared dataset.

  6. Create the recipe based upon the training dataset of the split.

  7. Add steps to the recipe such as adding calendar-based features and other tranformations and cleaning.

Recipe limitation for panel data

Generally, we can perform 1., 2., and 3. and preprocessing within a recipe by simply adding the correspong step to it.

The problem is that we are working with panel data , thus it is not possible to include 1., 2., and 3. as different steps within a recipe. We must first perfom these "steps" per group (by "Industry"), and then bind all groups into a single tibble again to be able to do 4., 5., 6. and 7.

p2_01_FE.png

Note: the "NA" values in the prepared dataset correspond the (unknown) Turnover predictions within the future frame.

Our data

Let us recall our timeseries dataset from Part 1: the Australian retail trade turnover dataset tsibble::aus_retail filtered by the “Australian Capital Territory” State value."

> aus_retail_tbl <- aus_retail %>%
     tk_tbl()

> monthly_retail_tbl <- aus_retail_tbl %>%
     filter(State == "Australian Capital Territory") %>%
     mutate(Month = as.Date(Month)) %>%
     mutate(Industry = as_factor(Industry)) %>%
     select(Month, Industry, Turnover)

> monthly_retail_tbl %>% glimpse()
Rows: 8,820
Columns: 3
$ Month    <date> 1982-04-01, 1982-05-01, 1982-06-01, 1982-07-01, 1982-08-01, 1982-09-01, 1982-10-01, 1982-11-01, 1982-12-01, 1983-01-01, 1983-02-01, 1983-03-01, 1983-04-01, 1983-05-01,~
$ Industry <fct> "Cafes, restaurants and catering services", "Cafes, restaurants and catering services", "Cafes, restaurants and catering services", "Cafes, restaurants and catering ser~
$ Turnover <dbl> 4.4, 3.4, 3.6, 4.0, 3.6, 4.2, 4.8, 5.4, 6.9, 3.8, 4.2, 4.0, 4.4, 4.3, 4.3, 4.6, 4.9, 4.5, 5.5, 6.7, 9.6, 3.8, 3.8, 4.4, 3.8, 3.9, 4.1, 4.3, 4.8, 5.0, 5.3, 5.6, 6.8, 4.6~

The code snippet below will steps (not recipe steps!) 1., 2., and 3. as explained above. You can always refer to Part 1 for function arguments explanations.

We will generate a list (groups) of 20 tibbles, each corresponding to a specific industry. In addition to steps 1., 2. and 3. I have also arranged in ascending order the timestamp (Month) for each Industry group for the lags, rolling lags and Fourier terms.

We add a future frame, starting from 2019-01-01 to 2019-12-01 (12-month length), which will hold NA values for the measure (Turnover) but to which the Fourier, lag and rolling lag features will also be added.

Finally, we bind them together to produce the first feature engineered tibble groups_fe_tbl.

groups <- lapply(X = 1:length(Industries), FUN = function(x){

  monthly_retail_tbl %>%
    filter(Industry == Industries[x]) %>%
    arrange(Month) %>%
    mutate(Turnover =  log1p(x = Turnover)) %>%
    mutate(Turnover =  standardize_vec(Turnover)) %>%
    future_frame(Month, .length_out = "12 months", .bind_data = TRUE) %>%
    mutate(Industry = Industries[x]) %>%
    tk_augment_fourier(.date_var = Month, .periods = 12, .K = 1) %>%
    tk_augment_lags(.value = Turnover, .lags = 12) %>%
    tk_augment_slidify(.value   = Turnover_lag12,
                       .f       = ~ mean(.x, na.rm = TRUE), 
                       .period  = c(3, 6, 9, 12),
                       .partial = TRUE,
                       .align   = "center")
})

groups_fe_tbl <- bind_rows(groups) %>%
  rowid_to_column(var = "rowid")

Below I will display the last 13 months of groups_fe_tbl, you should see the last historical row (2018-12-01) and the future frame values for the last industry (obviously, we have the same situation for other industries within the same tibble):

> groups_fe_tbl %>% tail(n=13) %>% glimpse()
Rows: 13
Columns: 16
$ rowid                  <int> 9048, 9049, 9050, 9051, 9052, 9053, 9054, 9055, 9056, 9057, 9058, 9059, 9060
$ Month                  <date> 2018-12-01, 2019-01-01, 2019-02-01, 2019-03-01, 2019-04-01, 2019-05-01, 2019-06-01, 2019-07-01, 20~
$ Industry               <chr> "Takeaway food services", "Takeaway food services", "Takeaway food services", "Takeaway food serv~
$ Turnover               <dbl> 2.078015, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA
$ Month_sin12_K1         <dbl> 0.16810089, 0.63846454, 0.93775213, 0.99301874, 0.80100109, 0.40982032, -0.10116832, -0.57126822, ~
$ Month_cos12_K1         <dbl> 0.98576980, 0.76965124, 0.34730525, -0.11795672, -0.59866288, -0.91216627, -0.99486932, -0.8207634~
$ Turnover_lag12         <dbl> 1.668977, 1.397316, 1.542003, 1.894938, 1.894938, 1.894938, 1.878972, 1.886969, 1.957669, 1.926527~
$ Turnover_lag13         <dbl> 1.755391, 1.668977, 1.397316, 1.542003, 1.894938, 1.894938, 1.894938, 1.878972, 1.886969, 1.957669~
$ Turnover_lag12_roll_3  <dbl> 1.607228, 1.536098, 1.611419, 1.777293, 1.894938, 1.889616, 1.886960, 1.907870, 1.923722, 1.972616~
$ Turnover_lag13_roll_3  <dbl> 1.723756, 1.607228, 1.536098, 1.611419, 1.777293, 1.894938, 1.889616, 1.886960, 1.907870, 1.923722~
$ Turnover_lag12_roll_6  <dbl> 1.667587, 1.692260, 1.715518, 1.750517, 1.832126, 1.901404, 1.906669, 1.929788, 1.941529, 1.974703~
$ Turnover_lag13_roll_6  <dbl> 1.668911, 1.667587, 1.692260, 1.715518, 1.750517, 1.832126, 1.901404, 1.906669, 1.929788, 1.941529~
$ Turnover_lag12_roll_9  <dbl> 1.750363, 1.744253, 1.741597, 1.757160, 1.779636, 1.808252, 1.878956, 1.925999, 1.946341, 1.952766~
$ Turnover_lag13_roll_9  <dbl> 1.748589, 1.750363, 1.744253, 1.741597, 1.757160, 1.779636, 1.808252, 1.878956, 1.925999, 1.929881~
$ Turnover_lag12_roll_12 <dbl> 1.783846, 1.784512, 1.785157, 1.787128, 1.811024, 1.828524, 1.862610, 1.904910, 1.941200, 1.946341~
$ Turnover_lag13_roll_12 <dbl> 1.766346, 1.783846, 1.784512, 1.785157, 1.787128, 1.811024, 1.828524, 1.843028, 1.887599, 1.925999~

Save standarization parameters

Since we are working with panel data (there are 20 time series, one for each Industry), consequently there will be 20 Standardization Parameters (mean, standard deviation) generated by standardize_vec(Turnover).

The values of these parameters must be saved, Ajay Mehta proposed an elegant solution with group_map(), I'm including his code snippet based on his comment at the end of this article. The mean and standard deviation values will be saved in std_mean and std_sd vectors respectively.

tmp <- monthly_retail_tbl %>%
  group_by(Industry) %>% 
  arrange(Month) %>%
  mutate(Turnover = log1p(x = Turnover)) %>% 
  group_map(~ c(mean = mean(.x$Turnover, na.rm = TRUE),
                sd = sd(.x$Turnover, na.rm = TRUE))) %>% 
  bind_rows()

std_mean <- tmp$mean
std_sd <- tmp$sd
rm('tmp')

Prepared and future datasets

We perform the split of groups_fe_tbl into prepared and future datasets .

The prepared dataset is the one with historic values of the measure (here Turnover) which must not have NA values. Although we know Turnover does have missing values we still apply filter(!is.na(Turnover)). We will also drop all rows with NA values elsewhere e.g., those introduced by lag and rolling lag features.

Finally, we obtain an 8,560-row tibble (data_prepared_tbl ) from a 9,060-row one (groups_fe_tbl). You can run a groups_fe_tbl %>% group_by(Industry) %>% skim() to check the total number of missing values (NA and NaN ) that were removed.

Prepared dataset starts from 1983-04-01.

> data_prepared_tbl <- groups_fe_tbl %>%
    filter(!is.na(Turnover)) %>%
    drop_na()

> data_prepared_tbl %>% glimpse()
Rows: 8,560
Columns: 16
$ rowid                  <int> 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38~
$ Month                  <date> 1983-05-01, 1983-06-01, 1983-07-01, 1983-08-01, 1983-09-01, 1983-10-01, 1983-11-01, 1983-12-01, 1~
$ Industry               <chr> "Cafes, restaurants and catering services", "Cafes, restaurants and catering services", "Cafes, re~
$ Turnover               <dbl> -2.0426107, -2.0426107, -1.9499231, -1.8620736, -1.9802554, -1.6990366, -1.4138382, -0.8757670, -2~
$ Month_sin12_K1         <dbl> 0.51455540, 0.87434662, 1.00000000, 0.86602540, 0.50000000, 0.01688948, -0.48530196, -0.84864426, ~
$ Month_cos12_K1         <dbl> 8.574571e-01, 4.853020e-01, 1.567970e-14, -5.000000e-01, -8.660254e-01, -9.998574e-01, -8.743466e-~
$ Turnover_lag12         <dbl> -2.355895, -2.281065, -2.140701, -2.281065, -2.074676, -1.890850, -1.725136, -1.370672, -2.209420,~
$ Turnover_lag13         <dbl> -2.011144, -2.355895, -2.281065, -2.140701, -2.281065, -2.074676, -1.890850, -1.725136, -1.370672,~
$ Turnover_lag12_roll_3  <dbl> -2.216035, -2.259220, -2.234277, -2.165481, -2.082197, -1.896888, -1.662219, -1.768409, -1.884923,~
$ Turnover_lag13_roll_3  <dbl> -2.183520, -2.216035, -2.259220, -2.234277, -2.165481, -2.082197, -1.896888, -1.662219, -1.768409,~
$ Turnover_lag12_roll_6  <dbl> -2.213974, -2.190758, -2.170709, -2.065582, -1.913850, -1.925303, -1.890905, -1.901909, -1.921958,~
$ Turnover_lag13_roll_6  <dbl> -2.197201, -2.213974, -2.190758, -2.170709, -2.065582, -1.913850, -1.925303, -1.890905, -1.901909,~
$ Turnover_lag12_roll_9  <dbl> -2.190758, -2.147914, -2.095067, -2.014578, -2.036609, -2.005362, -1.989766, -1.975371, -1.948876,~
$ Turnover_lag13_roll_9  <dbl> -2.213974, -2.190758, -2.147914, -2.095067, -2.014578, -2.036609, -2.005362, -1.989766, -1.975371,~
$ Turnover_lag12_roll_12 <dbl> -2.095067, -2.014578, -2.034063, -2.037755, -2.046334, -2.046334, -2.020226, -2.000355, -1.984457,~
$ Turnover_lag13_roll_12 <dbl> -2.147914, -2.095067, -2.014578, -2.034063, -2.037755, -2.046334, -2.046334, -2.020226, -2.000355,~

The future dataset hold the 12 future (and unknown) values of Turnover e.g., the NA values introduced by the future frame. We will predict and fill in the values in Part 6 with the final stacking ensemble model.

> future_tbl <- groups_fe_tbl %>%
    filter(is.na(Turnover))

> future_tbl %>% head(n=12) %>% glimpse()
Rows: 12
Columns: 16
$ rowid                  <int> 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453
$ Month                  <date> 2019-01-01, 2019-02-01, 2019-03-01, 2019-04-01, 2019-05-01, 2019-06-01, 2019-07-01, 2019-08-01, 20~
$ Industry               <chr> "Cafes, restaurants and catering services", "Cafes, restaurants and catering services", "Cafes, r~
$ Turnover               <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA
$ Month_sin12_K1         <dbl> 0.63846454, 0.93775213, 0.99301874, 0.80100109, 0.40982032, -0.10116832, -0.57126822, -0.90511451,~
$ Month_cos12_K1         <dbl> 0.76965124, 0.34730525, -0.11795672, -0.59866288, -0.91216627, -0.99486932, -0.82076344, -0.425167~
$ Turnover_lag12         <dbl> 1.125397, 1.282320, 1.569292, 1.425855, 1.437951, 1.429897, 1.550608, 1.469789, 1.457920, 1.469789~
$ Turnover_lag13         <dbl> 1.384894, 1.125397, 1.282320, 1.569292, 1.425855, 1.437951, 1.429897, 1.550608, 1.469789, 1.457920~
$ Turnover_lag12_roll_3  <dbl> 1.264204, 1.325670, 1.425822, 1.477699, 1.431234, 1.472819, 1.483431, 1.492773, 1.465833, 1.495272~
$ Turnover_lag13_roll_3  <dbl> 1.316081, 1.264204, 1.325670, 1.425822, 1.477699, 1.431234, 1.472819, 1.483431, 1.492773, 1.465833~
$ Turnover_lag12_roll_6  <dbl> 1.370952, 1.370952, 1.378452, 1.449320, 1.480565, 1.462003, 1.469326, 1.489352, 1.490694, 1.478711~
$ Turnover_lag13_roll_6  <dbl> 1.378274, 1.370952, 1.370952, 1.378452, 1.449320, 1.480565, 1.462003, 1.469326, 1.489352, 1.501243~
$ Turnover_lag12_roll_9  <dbl> 1.411416, 1.395927, 1.404907, 1.408445, 1.416559, 1.454825, 1.485468, 1.470874, 1.476502, 1.482009~
$ Turnover_lag13_roll_9  <dbl> 1.420563, 1.411416, 1.395927, 1.404907, 1.408445, 1.416559, 1.454825, 1.485468, 1.474990, 1.482009~
$ Turnover_lag12_roll_12 <dbl> 1.433627, 1.429420, 1.420139, 1.420139, 1.430152, 1.434573, 1.462680, 1.480716, 1.470874, 1.476502~
$ Turnover_lag13_roll_12 <dbl> 1.429496, 1.433627, 1.429420, 1.420139, 1.420139, 1.430152, 1.434266, 1.465153, 1.485468, 1.474990~

Train and test datasets split

We perform the split with timetk::time_series_split() function. The latter returns a split object (an R list) which can be further used for plotting both the traininig and the test datasets.

One important argument is assess, it determines the length of the test dataset. As per 2 , "The size of the test set is typically about 20% of the total sample, although this value depends on how long the sample is and how far ahead you want to forecast. The test dataset should ideally be at least as large as the maximum forecast horizon required." Since there are 428 observations (months) per industry, then 20% of 428 is 86 months (rounded up).

cumulative = TRUE means that we consider the training dataset from the begining (first row) of the prepared dataset.

> splits <- data_prepared_tbl %>%
    time_series_split(Month, 
                      assess = "86 months", 
                      cumulative = TRUE)

# split object
> splits
<Analysis/Assess/Total>
<6840/1720/8560>

Which means

  • 6840 rows for training (all 20 industries). Here, we only care about the rows since the number of columns will be updated by the recipe (explained soon).

  • 1720 rows for test dataset (86 months x 20 industries)

  • 8560 rows of data_prepared_tbl (all 20 industries)

You can confirm the size of the training and test datasets with the following code snippet using the tk_time_series_cv_plan() function which creates a new tibble with two new columns .id and .key

> splits %>%
   tk_time_series_cv_plan() %>% glimpse()
Rows: 8,560
Columns: 18
$ .id                    <chr> "Slice1", "Slice1", "Slice1", "Slice1", "Slice1", "Slice1", "Slice1", "Slice1", "Slice1", "Slic~
$ .key                   <fct> training, training, training, training, training, training, training, training, training, train~
$ rowid                  <int> 14, 467, 920, 1373, 1826, 2279, 2732, 3185, 3638, 4091, 4544, 4997, 5450, 5903, 6356, 6809, 726~
$ Month                  <date> 1983-05-01, 1983-05-01, 1983-05-01, 1983-05-01, 1983-05-01, 1983-05-01, 1983-05-01, 1983-05-01~
$ Industry               <chr> "Cafes, restaurants and catering services", "Cafes, restaurants and takeaway food services", "C~
$ Turnover               <dbl> -2.042611, -1.866127, -1.383977, -1.533674, -1.772910, -1.575604, -1.941481, -1.635071, -1.9995~
$ Month_sin12_K1         <dbl> 0.5145554, 0.5145554, 0.5145554, 0.5145554, 0.5145554, 0.5145554, 0.5145554, 0.5145554, 0.51455~
$ Month_cos12_K1         <dbl> 8.574571e-01, 8.574571e-01, 8.574571e-01, 8.574571e-01, 8.574571e-01, 8.574571e-01, 8.574571e-0~
$ Turnover_lag12         <dbl> -2.355895, -2.397538, -1.674143, -1.891997, -1.973832, -1.550576, -2.113506, -2.109525, -2.2924~
$ Turnover_lag13         <dbl> -2.011144, -2.195619, -1.713772, -1.891997, -2.044287, -1.680229, -2.080545, -2.035540, -2.4340~
$ Turnover_lag12_roll_3  <dbl> -2.216035, -2.299396, -1.771137, -1.971674, -2.053104, -1.602136, -2.091532, -2.110334, -2.4204~
$ Turnover_lag13_roll_3  <dbl> -2.183520, -2.296579, -1.693957, -1.891997, -2.009059, -1.615402, -2.097026, -2.072532, -2.3632~
$ Turnover_lag12_roll_6  <dbl> -2.213974, -2.275130, -1.824439, -2.023204, -2.202758, -1.587032, -2.072703, -2.125292, -2.4065~
$ Turnover_lag13_roll_6  <dbl> -2.197201, -2.278793, -1.787835, -1.988232, -2.146638, -1.577044, -2.076792, -2.110131, -2.4117~
$ Turnover_lag12_roll_9  <dbl> -2.190758, -2.244673, -1.826687, -2.015706, -2.222243, -1.607126, -2.067327, -2.110334, -2.4193~
$ Turnover_lag13_roll_9  <dbl> -2.213974, -2.275130, -1.824439, -2.023204, -2.202758, -1.587032, -2.072703, -2.125292, -2.4065~
$ Turnover_lag12_roll_12 <dbl> -2.095067, -2.152525, -1.819038, -1.978000, -2.167591, -1.593199, -2.043333, -2.048362, -2.3049~
$ Turnover_lag13_roll_12 <dbl> -2.147914, -2.200931, -1.828293, -2.002079, -2.217778, -1.606260, -2.056832, -2.089405, -2.3594~

> splits %>%
   tk_time_series_cv_plan() %>% select(.id) %>% table()

Slice1 
  8560 

> splits %>%
   tk_time_series_cv_plan() %>% select(.key) %>% table()

training  testing 
    6840     1720

Let us plot the training and test datasets for one Industry value with the plot_time_series_cv_plan() function:

> Industries <- unique(aus_retail_tbl$Industry)
> splits %>%
  tk_time_series_cv_plan() %>%
  filter(Industry == Industries[1]) %>%
  plot_time_series_cv_plan(.date_var = Month, 
                           .value = Turnover, 
                           .title = paste0("Split for ", industry))

p2_splits_plot.png

You can extract the training and test datasets by applying the training() and testing() functions respectively on the splits object.

> training(splits) %>% 
     filter(Industry == Industries[industry]) %>% 
     glimpse()
Rows: 342
Columns: 16
$ rowid                  <int> 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,~
$ Month                  <date> 1983-05-01, 1983-06-01, 1983-07-01, 1983-08-01, 1983-09-01, 1983-10-01, 1983-11-01, 1983-12-01~
$ Industry               <chr> "Cafes, restaurants and catering services", "Cafes, restaurants and catering services", "Cafes,~
$ Turnover               <dbl> -2.0426107, -2.0426107, -1.9499231, -1.8620736, -1.9802554, -1.6990366, -1.4138382, -0.8757670,~
$ Month_sin12_K1         <dbl> 0.51455540, 0.87434662, 1.00000000, 0.86602540, 0.50000000, 0.01688948, -0.48530196, -0.8486442~
$ Month_cos12_K1         <dbl> 8.574571e-01, 4.853020e-01, 1.567970e-14, -5.000000e-01, -8.660254e-01, -9.998574e-01, -8.74346~
$ Turnover_lag12         <dbl> -2.355895, -2.281065, -2.140701, -2.281065, -2.074676, -1.890850, -1.725136, -1.370672, -2.2094~
$ Turnover_lag13         <dbl> -2.011144, -2.355895, -2.281065, -2.140701, -2.281065, -2.074676, -1.890850, -1.725136, -1.3706~
$ Turnover_lag12_roll_3  <dbl> -2.216035, -2.259220, -2.234277, -2.165481, -2.082197, -1.896888, -1.662219, -1.768409, -1.8849~
$ Turnover_lag13_roll_3  <dbl> -2.183520, -2.216035, -2.259220, -2.234277, -2.165481, -2.082197, -1.896888, -1.662219, -1.7684~
$ Turnover_lag12_roll_6  <dbl> -2.213974, -2.190758, -2.170709, -2.065582, -1.913850, -1.925303, -1.890905, -1.901909, -1.9219~
$ Turnover_lag13_roll_6  <dbl> -2.197201, -2.213974, -2.190758, -2.170709, -2.065582, -1.913850, -1.925303, -1.890905, -1.9019~
$ Turnover_lag12_roll_9  <dbl> -2.190758, -2.147914, -2.095067, -2.014578, -2.036609, -2.005362, -1.989766, -1.975371, -1.9488~
$ Turnover_lag13_roll_9  <dbl> -2.213974, -2.190758, -2.147914, -2.095067, -2.014578, -2.036609, -2.005362, -1.989766, -1.9753~
$ Turnover_lag12_roll_12 <dbl> -2.095067, -2.014578, -2.034063, -2.037755, -2.046334, -2.046334, -2.020226, -2.000355, -1.9844~
$ Turnover_lag13_roll_12 <dbl> -2.147914, -2.095067, -2.014578, -2.034063, -2.037755, -2.046334, -2.046334, -2.020226, -2.0003~

It is important to understand that the actual training dataset to be used to learn the model will be the output of the recipe , as explained hereafter.

Create recipe

As stated in Get Started, before training the model, we can use a recipe to create (new) predictors and conduct the preprocessing required by the model.

recipe(Turnover ~ ., data = training(splits))

"The recipe() function as we used it here has two arguments:

  • A formula. Any variable on the left-hand side of the tilde (~) is considered the model outcome (here, Turnover). On the right-hand side of the tilde are the predictors. Variables may be listed by name, or you can use the dot (.) to indicate all other variables as predictors.

  • The data. A recipe is associated with the data set used to create the model. This will typically be the training set, so data = training(splits) here. Naming a data set doesn’t actually change the data itself; it is only used to catalog the names of the variables and their types, like factors, integers, dates, etc."

For each preprocessing and feature engineering steps explained in Part 1 there is a specific step_<> :

  • add calendar-based features: step_timeseries_signature()

  • remove columns: step_rm()

  • perform one-hot encoding on categorical variables: step_dummy()

  • normalize numeric variables: step_normalize()

recipe_spec <- recipe(Turnover ~ ., data = training(splits)) %>%
  update_role(rowid, new_role = "indicator") %>%  
  step_other(Industry) %>%
  step_timeseries_signature(Month) %>%
  step_rm(matches("(.xts$)|(.iso$)|(hour)|(minute)|(second)|(day)|(week)|(am.pm)")) %>%
  step_dummy(all_nominal(), one_hot = TRUE) %>%
  step_normalize(Month_index.num, Month_year)

Notice the two very useful functions explained hereafter:

update_role(rowid, new_role = "indicator"): it lets recipes know that rowid variable has a custom role called "indicator" (a role can have any character value). This tells the recipe to keep the rowid variable without using it neither as an outcome nor a predictor.

step_other(Industry): in case the training dataset does not include all Industry (categorical variable) values, it will pool the infrequently occurring values into an "other" category.

Below we display the recipe's spec, its summary, and finally the (real) training dataset:

> recipe_spec
Recipe

Inputs:

      role #variables
 indicator          1
   outcome          1
 predictor         14

Operations:

Collapsing factor levels for Industry
Timeseries signature features from Month
Delete terms matches("(.xts$)|(.iso$)|(hour)|(minute)|(second)|(day)|(week)|(am.pm)")
Dummy variables from all_nominal()
Centering and scaling for Month_index.num, Month_year

# Recipe summary: 50 predictors (including all dummy variables) + 1 indicator + 1 outcome
> recipe_spec %>% 
+   prep() %>% 
+   summary()
# A tibble: 52 x 4
   variable              type    role      source  
   <chr>                 <chr>   <chr>     <chr>   
 1 rowid                 numeric indicator original
 2 Month                 date    predictor original
 3 Month_sin12_K1        numeric predictor original
 4 Month_cos12_K1        numeric predictor original
 5 Turnover_lag12        numeric predictor original
 6 Turnover_lag13        numeric predictor original
 7 Turnover_lag12_roll_3 numeric predictor original
 8 Turnover_lag13_roll_3 numeric predictor original
 9 Turnover_lag12_roll_6 numeric predictor original
10 Turnover_lag13_roll_6 numeric predictor original
# ... with 42 more rows

# Extracted training dataset
> recipe_spec %>% 
   prep() %>%
   juice() %>% 
   glimpse()
Rows: 6,840
Columns: 52
$ rowid                                                                      <int> 14, 467, 920, 1373, 1826, 2279, 2732, 3185,~
$ Month                                                                      <date> 1983-05-01, 1983-05-01, 1983-05-01, 1983-0~
$ Month_sin12_K1                                                             <dbl> 0.5145554, 0.5145554, 0.5145554, 0.5145554,~
$ Month_cos12_K1                                                             <dbl> 8.574571e-01, 8.574571e-01, 8.574571e-01, 8~
$ Turnover_lag12                                                             <dbl> -2.355895, -2.397538, -1.674143, -1.891997,~
$ Turnover_lag13                                                             <dbl> -2.011144, -2.195619, -1.713772, -1.891997,~
$ Turnover_lag12_roll_3                                                      <dbl> -2.216035, -2.299396, -1.771137, -1.971674,~
$ Turnover_lag13_roll_3                                                      <dbl> -2.183520, -2.296579, -1.693957, -1.891997,~
$ Turnover_lag12_roll_6                                                      <dbl> -2.213974, -2.275130, -1.824439, -2.023204,~
$ Turnover_lag13_roll_6                                                      <dbl> -2.197201, -2.278793, -1.787835, -1.988232,~
$ Turnover_lag12_roll_9                                                      <dbl> -2.190758, -2.244673, -1.826687, -2.015706,~
$ Turnover_lag13_roll_9                                                      <dbl> -2.213974, -2.275130, -1.824439, -2.023204,~
$ Turnover_lag12_roll_12                                                     <dbl> -2.095067, -2.152525, -1.819038, -1.978000,~
$ Turnover_lag13_roll_12                                                     <dbl> -2.147914, -2.200931, -1.828293, -2.002079,~
$ Turnover                                                                   <dbl> -2.042611, -1.866127, -1.383977, -1.533674,~
$ Month_index.num                                                            <dbl> -1.727262, -1.727262, -1.727262, -1.727262,~
$ Month_year                                                                 <dbl> -1.710287, -1.710287, -1.710287, -1.710287,~
$ Month_half                                                                 <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
$ Month_quarter                                                              <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2~
$ Month_month                                                                <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5~
$ Industry_Cafes..restaurants.and.catering.services                          <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Cafes..restaurants.and.takeaway.food.services                     <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Clothing.retailing                                                <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Clothing..footwear.and.personal.accessory.retailing               <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Department.stores                                                 <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Electrical.and.electronic.goods.retailing                         <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Food.retailing                                                    <dbl> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Footwear.and.other.personal.accessory.retailing                   <dbl> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Furniture..floor.coverings..houseware.and.textile.goods.retailing <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0~
$ Industry_Hardware..building.and.garden.supplies.retailing                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0~
$ Industry_Household.goods.retailing                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0~
$ Industry_Liquor.retailing                                                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0~
$ Industry_Newspaper.and.book.retailing                                      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0~
$ Industry_Other.recreational.goods.retailing                                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0~
$ Industry_Other.retailing                                                   <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1~
$ Industry_Other.retailing.n.e.c.                                            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Other.specialised.food.retailing                                  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Pharmaceutical..cosmetic.and.toiletry.goods.retailing             <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Supermarket.and.grocery.stores                                    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Industry_Takeaway.food.services                                            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_01                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_02                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_03                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_04                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_05                                                         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1~
$ Month_month.lbl_06                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_07                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_08                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_09                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_10                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_11                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~
$ Month_month.lbl_12                                                         <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0~

As for the standardization parameters, we save the normalization parameters generated by the step_normalize() function applied on Month_index.num and Month_year columns. Min-max normalization was used for these columns.

To get the original min and max values prior normalization you must customize the skim() function as I did in Part 1 with my myskim() functionand apply it to training(splits).

recipe(Turnover ~ ., data = training(splits)) %>%
  step_timeseries_signature(Month) %>%
  prep() %>%
  juice() %>% 
  myskim()

# Extract of the output for the sake of clarity:
... ...

-- Variable type: numeric ------------------------------------------------------------------------------------------------------------------------------------------------------------------------
   skim_variable          n_missing complete_rate     mean            sd           p0           p25      p50     p75          p100           max          min

 ... ...

15 Month_index.num                0             1  8.69e+8 259648888.    420595200    644198400      8.69e+8 1.09e+9 1317427200    1317427200    420595200   
16 Month_year                     0             1  2.00e+3         8.23       1983         1990      2.00e+3 2.00e+3       2011          2011         1983  

... ...

# Month_index.num normalization Parameters
Month_index.num_limit_lower = 420595200
Month_index.num_limit_upper = 1317427200

# Month_year normalization Parameters
Month_year_limit_lower = 1983
Month_year_limit_upper = 2011

Save your work

We will keep all data, recipes, splits, parameters and any other useful stuff in a list and save it as an RDS file.

feature_engineering_artifacts_list <- list(
  # Data
  data = list(
    data_prepared_tbl = data_prepared_tbl,
    future_tbl      = future_tbl,
    industries = Industries
  ),

  # Recipes
  recipes = list(
    recipe_spec = recipe_spec
  ),

  # Splits
  splits = splits,

  # Inversion Parameters
  standardize = list(
    std_mean = std_mean,
    std_sd   = std_sd
  ),

  normalize = list(
    Month_index.num_limit_lower = Month_index.num_limit_lower, 
    Month_index.num_limit_upper = Month_index.num_limit_upper,
    Month_year_limit_lower = Month_year_limit_lower,
    Month_year_limit_upper = Month_year_limit_upper
  )  
)

feature_engineering_artifacts_list %>% 
  write_rds("feature_engineering_artifacts_list.rds")

Conclusion

In this article I have explained how to perform time series feature engineering with recipes thanks to the timetk package.

The recipe performed preprocessing, transformations and feature engineering applicable to all 20 time series (panel data). For operations such as measure standardization, adding lags, rolling lags and Fourier terms we needed to code a specific by group (here, by Industry) feature engineering.

The recipe is an important step in the data science process with modeltime since it will be added to a workflow along with a model, this will be explained in my next article (Part 3).

Note: I haven't shown all the power of recipes e.g., model-specific recipes can be built from a "base" recipe which can be modified in order to adapt to the model requirements. Remember, it is important to perform experimentation with different datasets and creating interdependent recipes can be very useful.

References

1 Dancho, Matt, "DS4B 203-R: High-Performance Time Series Forecasting", Business Science University

2 Hyndman R., Athanasopoulos G., "Forecasting: Principles and Practice", 3rd edition