Demo Week: Tidy Forecasting with sweep” is an excellent article that uses tidy methods with time series. This article uses their analysis with rsample to get performance estimates for future observations using rolling forecast origin resampling.

The data are sales of alcoholic beverages originally from the Federal Reserve Bank of St. Louis website.

library(tidymodels)
## ── Attaching packages ───────────────────────────────── tidymodels 0.0.1 ──
## ✔ tibble    1.4.2     ✔ yardstick 0.0.2
## ✔ recipes   0.1.4     ✔ infer     0.3.1
## ✔ broom     0.5.0
## ── Conflicts ──────────────────────────────────── tidymodels_conflicts() ──
## ✖ yardstick::accuracy() masks forecast::accuracy()
## ✖ rsample::fill()       masks tidyr::fill()
## ✖ dplyr::filter()       masks stats::filter()
## ✖ dplyr::lag()          masks stats::lag()
## ✖ recipes::step()       masks stats::step()
## ✖ yardstick::tidy()     masks broom::tidy(), recipes::tidy(), rsample::tidy()
data("drinks")
str(drinks, give.att = FALSE)
## Classes 'tbl_df', 'tbl' and 'data.frame':    309 obs. of  2 variables:
##  $ date          : Date, format: "1992-01-01" "1992-02-01" ...
##  $ S4248SM144NCEN: num  3459 3458 4002 4564 4221 ...

Each row is a month of sales (in millions of US dollars).

Suppose that predictions for one year ahead were needed and the model should use the most recent data from the last 20 years. To setup this resampling scheme:

## [1] 58
## # Rolling origin forecast resampling 
## # A tibble: 58 x 2
##    splits           id     
##    <list>           <chr>  
##  1 <split [240/12]> Slice01
##  2 <split [240/12]> Slice02
##  3 <split [240/12]> Slice03
##  4 <split [240/12]> Slice04
##  5 <split [240/12]> Slice05
##  6 <split [240/12]> Slice06
##  7 <split [240/12]> Slice07
##  8 <split [240/12]> Slice08
##  9 <split [240/12]> Slice09
## 10 <split [240/12]> Slice10
## # ... with 48 more rows

Each split element contains the information about that resample:

## <240/12/309>

For plotting, let’s index each split by the first day of the assessment set:

## [1] "2012-01-01" "2012-02-01" "2012-03-01" "2012-04-01" "2012-05-01"
## [6] "2012-06-01"

This resampling scheme has 58 splits of the data so that there will be 58 ARIMA models that are fit. To create the models, the auto.arima function from the forecast package is used. The functions analysis and assessment return the data frame, so another step converts the data in to a ts object called mod_dat using a function in the timetk package.

Each model is saved in a new column:

## Series: . 
## ARIMA(2,1,2)(2,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2    ma1     ma2   sar1    sar2    sma1
##       -0.994  -0.502  0.024  -0.414  0.404  -0.334  -0.554
## s.e.   0.106   0.081  0.111   0.125  0.100   0.076   0.083
## 
## sigma^2 estimated as 71958:  log likelihood=-1591
## AIC=3199   AICc=3199   BIC=3226

(There are some warnings produced by these first regarding extra columns in the data that can be ignored)

Using the model fits, performance will be measured in two ways:

  • interpolation error will measure how well the model fits to the data that were used to create the model. This is most likely optimistic since no holdout method is used.
  • extrapolation or forecast error evaluates the efficacy of the model on the data from the following year (that were not used in the model fit).

In each case, the mean absolute percent error (MAPE) is the statistic used to characterize the model fits. The interpolation error can be computed from the Arima object. to make things easy, the sweep package’s sw_glance function is used:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.84    2.89    2.92    2.93    2.95    3.13

For the extrapolation error, the model and split objects are required. Using these:

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    2.37    3.19    3.56    3.65    4.16    5.45

What do these error estimates look like over time?

It is likely that the interpolation error is an underestimate to some degree.

It is also worth noting that rolling_origin() can be used over calendar periods, rather than just over a fixed window size. This is especially useful for irregular series where a fixed window size might not make sense because of missing data points, or because of calendar features like different months having a different number of days.

The example below demonstrates this idea by splitting drinks into a nested set of 26 years, and rolling over years rather than months. Note that the end result accomplishes a different task than the original example, in this case, each slice moves forward an entire year, rather than just one month.

## # A tibble: 20 x 2
##     year data             
##    <dbl> <list>           
##  1  1992 <tibble [12 × 2]>
##  2  1993 <tibble [12 × 2]>
##  3  1994 <tibble [12 × 2]>
##  4  1995 <tibble [12 × 2]>
##  5  1996 <tibble [12 × 2]>
##  6  1997 <tibble [12 × 2]>
##  7  1998 <tibble [12 × 2]>
##  8  1999 <tibble [12 × 2]>
##  9  2000 <tibble [12 × 2]>
## 10  2001 <tibble [12 × 2]>
## 11  2002 <tibble [12 × 2]>
## 12  2003 <tibble [12 × 2]>
## 13  2004 <tibble [12 × 2]>
## 14  2005 <tibble [12 × 2]>
## 15  2006 <tibble [12 × 2]>
## 16  2007 <tibble [12 × 2]>
## 17  2008 <tibble [12 × 2]>
## 18  2009 <tibble [12 × 2]>
## 19  2010 <tibble [12 × 2]>
## 20  2011 <tibble [12 × 2]>

The workflow to access these calendar slices is to use bind_rows() to join each analysis set together.

## # Rolling origin forecast resampling 
## # A tibble: 6 x 3
##   splits         id     extracted_slice   
## * <list>         <chr>  <list>            
## 1 <split [20/1]> Slice1 <tibble [240 × 2]>
## 2 <split [20/1]> Slice2 <tibble [240 × 2]>
## 3 <split [20/1]> Slice3 <tibble [240 × 2]>
## 4 <split [20/1]> Slice4 <tibble [240 × 2]>
## 5 <split [20/1]> Slice5 <tibble [240 × 2]>
## 6 <split [20/1]> Slice6 <tibble [240 × 2]>