#### James Young

## Packages Used

**library**(tidyverse)
**library**(patchwork)
**library**(tidyquant)
**library**(timetk)
**library**(sweep)
**library**(forecast)
**library**(knitr)

## Purpose

In the last article we completed an EDA of the sales data and determined it was most likely in units sold rather than dollars and that the data had a hierarchical and grouped nature. Here we will forecast at the dept_id crossed with state_id level. That will allow us to roll up to higher levels including dept_id, state_id, and total sales. We want to try a couple different models, including exponential smoothing and tbats from the forecast package and compare them against “control” models including naive, seasonal naive, and mean based forecasts. Ideally we will find a model type that performs best across all levels and can serve as our base model in production. Once we compare these models and find a base model, we can also look deeper into the results to try to build ideas for further improvement upon our base model.

## Set-Up the Data

Since we will be forecasting multiple series that are part of the hierarchy, we want to organize the data in a way that can store information in a way that we can reference what the training data was, what the forecasts are, what the model type was and its parameters, and other important identifying information for the time-series we are forecasting.

I will start by pivoting the data from wide to long by state_id, dept_id, and date with the corresponding sales in a single column. Next I will compress the sales data from its current daily format to a pseudo-monthly format by summing in 30 day increments. After these operations are done, I apply a glimpse to the new data format and we can see it is now 4 columns and 1,344 rows.

```
long_data <- data %>%
group_by(dept_id,state_id) %>%
summarise(across(where(is.numeric), sum, na.rm = T)) %>%
pivot_longer(cols = !ends_with("id"), names_to = "day", values_to = "sales") %>%
mutate(day = as.numeric(gsub("d_", "", day))) %>%
ungroup() %>%
group_by(dept_id, state_id, month=cut(day, breaks= seq(1, 1940, by = 30), labels = F) ) %>%
summarise(sales= sum(sales, na.rm = T)) %>%
arrange(as.numeric(month)) %>%
filter(!is.na(month))
glimpse(long_data)
```

```
## Rows: 1,344
## Columns: 4
## Groups: dept_id, state_id [21]
## $ dept_id <fct> FOODS_1, FOODS_1, FOODS_1, FOODS_2, FOODS_2, FOODS_2, FOOD...
## $ state_id <fct> CA, TX, WI, CA, TX, WI, CA, TX, WI, CA, TX, WI, CA, TX, WI...
## $ month <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ sales <int> 29383, 12677, 15968, 44610, 35033, 25659, 165886, 119593, ...
```

Now I will create a new column of train data for each row, with all known sales up to that date. To contain the required train data in a single cell for a given row I will use a list as columns structure. I want to make sure there is enough train data so I will only start making forecasts on the 45th month of data and will make 15 forecasts from there out. These 15 forecasts will further be compared against the actual sales to give an idea of what level of accuracy could be expected in future forecasts and this will also allow the comparison of our different model types.

`train <- NULL`

for(iin45:60) { long_data2 <- long_data %>% filter(month < i) long_data2$month <- long_data2$month*30 long_data2$month <- as.Date(long_data2$month) ts <- long_data2 %>% group_by(dept_id, state_id) %>% nest() train <- rbind(train, ts) }`actuals <- NULL for (i in 45:60) { long_data2 <- long_data long_data2$month <- long_data2$month*30 long_data2$month <- as.Date(long_data2$month) ts <- long_data2 %>% group_by(dept_id, state_id) %>% nest() actuals <- rbind(actuals, ts) } train$train <- train$data train$data <- NULL train$actuals <- actuals$data`

Let’s look at where the data is at now. We have a column that signifies the state_id and another column for dept_id. We also have our lists as columns for train data and actual sales.

`glimpse(train)`

```
## Rows: 336
## Columns: 4
## Groups: dept_id, state_id [21]
## $ dept_id <fct> FOODS_1, FOODS_1, FOODS_1, FOODS_2, FOODS_2, FOODS_2, FOOD...
## $ state_id <fct> CA, TX, WI, CA, TX, WI, CA, TX, WI, CA, TX, WI, CA, TX, WI...
## $ train <list> [<tbl_df[44 x 2]>, <tbl_df[44 x 2]>, <tbl_df[44 x 2]>, <t...
## $ actuals <list> [<tbl_df[64 x 2]>, <tbl_df[64 x 2]>, <tbl_df[64 x 2]>, <t...
```

Looking at an example of the train column and actual sales column, we can see they have the corresponding dates. We can also see that the actual sales has observations going out further into the future which is necessary to compare with the forecasts for those corresponding times. Since we compressed the data every 30 days, notice how we fluctuate between the first of the month to the last of the month periodically.

`train[["train"]][[1]]`

```
## # A tibble: 44 x 2
## month sales
## <date> <int>
## 1 1970-01-31 29383
## 2 1970-03-02 26431
## 3 1970-04-01 25321
## 4 1970-05-01 24078
## 5 1970-05-31 24839
## 6 1970-06-30 27553
## 7 1970-07-30 23682
## 8 1970-08-29 26562
## 9 1970-09-28 30354
## 10 1970-10-28 28668
## # ... with 34 more rows
```

`train[["actuals"]][[1]]`

```
## # A tibble: 64 x 2
## month sales
## <date> <int>
## 1 1970-01-31 29383
## 2 1970-03-02 26431
## 3 1970-04-01 25321
## 4 1970-05-01 24078
## 5 1970-05-31 24839
## 6 1970-06-30 27553
## 7 1970-07-30 23682
## 8 1970-08-29 26562
## 9 1970-09-28 30354
## 10 1970-10-28 28668
## # ... with 54 more rows
```

Now I will convert the lists as columns to time series format. We will say the series started in 1970 as an arbitrary starting point for this example. We will also report this data as having a frequency of 12 meaning monthly data. Again, when we made the breaks above it was every 30 days so it is really pseudo-monthly, but will be very close.

```
train_tbl_ts <- train %>%
mutate(train.ts = map(.x = train,
.f = tk_ts,
select = -month,
start = 1970,
freq = 12),
actuals.ts = map(.x = actuals,
.f = tk_ts,
select = -month,
start = 1970,
freq = 12))
```

Again, we’ll check out the updated table we are building. We have added two new columns, train.ts and actual sales.ts. Looking at examples of these new columns we see that the organization of the data has changed now that it is formatted as time series relative to the format we saw above.

`head(train_tbl_ts)`

```
## # A tibble: 6 x 6
## # Groups: dept_id, state_id [6]
## dept_id state_id train actuals train.ts actuals.ts
## <fct> <fct> <list> <list> <list> <list>
## 1 FOODS_1 CA <tibble [44 x 2]> <tibble [64 x 2]> <ts [44 x 1~ <ts [64 x 1~
## 2 FOODS_1 TX <tibble [44 x 2]> <tibble [64 x 2]> <ts [44 x 1~ <ts [64 x 1~
## 3 FOODS_1 WI <tibble [44 x 2]> <tibble [64 x 2]> <ts [44 x 1~ <ts [64 x 1~
## 4 FOODS_2 CA <tibble [44 x 2]> <tibble [64 x 2]> <ts [44 x 1~ <ts [64 x 1~
## 5 FOODS_2 TX <tibble [44 x 2]> <tibble [64 x 2]> <ts [44 x 1~ <ts [64 x 1~
## 6 FOODS_2 WI <tibble [44 x 2]> <tibble [64 x 2]> <ts [44 x 1~ <ts [64 x 1~
```

`train_tbl_ts[["train.ts"]][[1]]`

```
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1970 29383 26431 25321 24078 24839 27553 23682 26562 30354 28668 34134 28916
## 1971 37051 36916 39814 34074 38221 38615 41239 43200 37438 35174 34089 35964
## 1972 38025 38081 37044 36384 44320 43742 46269 45282 39207 34923 36288 39130
## 1973 39292 39368 39065 36508 38213 37984 38002 40269
```

`train_tbl_ts[["actuals.ts"]][[1]]`

```
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 1970 29383 26431 25321 24078 24839 27553 23682 26562 30354 28668 34134 28916
## 1971 37051 36916 39814 34074 38221 38615 41239 43200 37438 35174 34089 35964
## 1972 38025 38081 37044 36384 44320 43742 46269 45282 39207 34923 36288 39130
## 1973 39292 39368 39065 36508 38213 37984 38002 40269 38349 34771 34718 40430
## 1974 41020 38059 37743 36241 41732 44178 43870 44190 41671 38848 35935 45646
## 1975 43603 42202 42884 41968
```

Now we have arrived at the fun part!! I have completed the set-up of the data for building models for different time series at different points in time. Below I apply 5 different model types, exponential smoothing (ets), tbats which is handy for data with complex seasonality, and our control models including mean, naive, and seasonal naive. We can apply all of these below in a few lines of code and save the models as new columns. I will pre-preemptively tell our control models to forecast 3 months into the future, approximately one quarter. We will make sure our ets and tbats model on forecast 3 months ahead as well in a step below.

```
start <- Sys.time()
train_tbl_ts <- train_tbl_ts %>%
mutate(fit.ets = map(train.ts, ets),
fit.tbats = map(train.ts, tbats),
fit.snaive = map(train.ts, snaive, h =3),
fit.naive = map(train.ts, naive, h =3),
fit.mean = map(train.ts, meanf, h =3))
end <- Sys.time()
cat("How Long Did It Take?", start - end, "minutes")
```

`## How Long Did It Take? -8.83988 minutes`

That took a decent amount of time to make forecast models for the few time series we consider here over 15 unique points in time. The more granular we go in the hierarchy, the longer it will take to build the increased model count. Building more model types or trying external regressors with models like arima or nnetar can add even more time but may be worth it for an adequate accuracy increase.

Taking a look at our updated table we see it now has the models saved as columns for each observation.

`head(train_tbl_ts)`

```
## # A tibble: 6 x 11
## # Groups: dept_id, state_id [6]
## dept_id state_id train actuals train.ts actuals.ts fit.ets fit.tbats
## <fct> <fct> <lis> <list> <list> <list> <list> <list>
## 1 FOODS_1 CA <tib~ <tibbl~ <ts [44~ <ts [64 x~ <ets> <bats>
## 2 FOODS_1 TX <tib~ <tibbl~ <ts [44~ <ts [64 x~ <ets> <bats>
## 3 FOODS_1 WI <tib~ <tibbl~ <ts [44~ <ts [64 x~ <ets> <bats>
## 4 FOODS_2 CA <tib~ <tibbl~ <ts [44~ <ts [64 x~ <ets> <tbats>
## 5 FOODS_2 TX <tib~ <tibbl~ <ts [44~ <ts [64 x~ <ets> <tbats>
## 6 FOODS_2 WI <tib~ <tibbl~ <ts [44~ <ts [64 x~ <ets> <bats>
## # ... with 3 more variables: fit.snaive <list>, fit.naive <list>,
## # fit.mean <list>
```

Below, I use the trained models to make forecasts that represent the coming 3 months of sales for the respective time series we are forecasting. I also take the corresponding actual sales and date that the forecast was made. We will use this information to evaluate performance of our models.

#This function will take the actual sales that line up with the forecasted sales for a proper comparison and sum them into into one value representing the whole 3 month forecasted period.actuals_sum <-function(x, y) { sum(subset(x, x$month > max(y$month))[1:3,2]) }#This function sums the mean, or point forecast, of the forecasts from 3 months into one value representing the whole 3 month forecasted period.forecasts_sum <-function(x) { sum(x[["mean"]]) }#This function will strip out the last date in the train data so we can keep track of performance at different points in time.str_date <-function(x) { max(x$month) }`#I use map and map2 functions below to apply a through a column. train_tbl_ts2 <- train_tbl_ts %>% mutate(fcast.ets = map(fit.ets, forecast, h = 3), fcast.tbats = map(fit.tbats, forecast, h = 3), fcst.ets_sum = map_dbl(fcast.ets, forecasts_sum), fcst.tbats_sum = map_dbl(fcast.tbats, forecasts_sum), fcst.mean_sum = map_dbl(fit.mean, forecasts_sum), fcst.snaive_sum = map_dbl(fit.snaive, forecasts_sum), fcst.naive_sum = map_dbl(fit.naive, forecasts_sum), fcst.meanf_sum = map(fcast.tbats, forecasts_sum), actuals_sum = map2_dbl(actuals, train, actuals_sum), date = map_dbl(train, str_date))`

If we glimpse the data we can see all the new columns we just added to the data frame. This will make it convenient to save as a RDS to look back on and redistribute without having to rerun all the code above.

`glimpse(train_tbl_ts2)`

```
## Rows: 336
## Columns: 21
## Groups: dept_id, state_id [21]
## $ dept_id <fct> FOODS_1, FOODS_1, FOODS_1, FOODS_2, FOODS_2, FOODS_...
## $ state_id <fct> CA, TX, WI, CA, TX, WI, CA, TX, WI, CA, TX, WI, CA,...
## $ train <list> [<tbl_df[44 x 2]>, <tbl_df[44 x 2]>, <tbl_df[44 x ...
## $ actuals <list> [<tbl_df[64 x 2]>, <tbl_df[64 x 2]>, <tbl_df[64 x ...
## $ train.ts <list> [<ts[44 x 1]>, <ts[44 x 1]>, <ts[44 x 1]>, <ts[44 ...
## $ actuals.ts <list> [<ts[64 x 1]>, <ts[64 x 1]>, <ts[64 x 1]>, <ts[64 ...
## $ fit.ets <list> [<-437.2786, 880.5573, 885.9098, 881.1573, 9743508...
## $ fit.tbats <list> [<NULL, 0.7851275, NULL, NULL, NULL, NULL, NULL, 8...
## $ fit.snaive <list> [<Seasonal naive method, 29383, 26431, 25321, 2407...
## $ fit.naive <list> [<Naive method, 29383, 26431, 25321, 24078, 24839,...
## $ fit.mean <list> [<Mean, 80, 95, 29383, 26431, 25321, 24078, 24839,...
## $ fcast.ets <list> [<-437.2786, 880.5573, 885.9098, 881.1573, 9743508...
## $ fcast.tbats <list> [<0.7851275, 875.0101, 0, 9844305, 879.0101, 0.785...
## $ fcst.ets_sum <dbl> 119187.240, 53978.199, 65650.080, 133621.939, 96223...
## $ fcst.tbats_sum <dbl> 119342.658, 53352.993, 65662.145, 127589.767, 94219...
## $ fcst.mean_sum <dbl> 106984.909, 57611.250, 63413.523, 130268.864, 99066...
## $ fcst.snaive_sum <dbl> 110418, 62274, 64834, 118661, 89780, 109986, 710521...
## $ fcst.naive_sum <dbl> 120807, 53559, 66063, 134037, 96150, 131637, 727974...
## $ fcst.meanf_sum <list> [119342.7, 53352.99, 65662.15, 127589.8, 94219.05,...
## $ actuals_sum <dbl> 107838, 57818, 64674, 143109, 101642, 141617, 67387...
## $ date <dbl> 1320, 1320, 1320, 1320, 1320, 1320, 1320, 1320, 132...
```

Below I make a function to get a quick and dirty approximation of the mean absolute percent error (MAPE) of the different model types. A lower MAPE represents a better model. These MAPE scores can be interpreted as a MAPE of 0.1 equating to 10 % mean absolute percent error, 0.2 equates to 20 % mean absolute percent error, and so on.

`MAPE_performance <- `**function**(tb_ts22, **...**) {
tb_ts22 %>%
group_by(!!!enquos(**...**), date) %>%
summarise(across(where(is.numeric), sum)) %>%
ungroup() %>%
mutate(ets_AE = abs(actuals_sum - fcst.ets_sum),
snaive_AE = abs(actuals_sum - fcst.snaive_sum),
naive_AE = abs(actuals_sum - fcst.naive_sum),
mean_AE = abs(actuals_sum - fcst.mean_sum),
tbats_AE = abs(actuals_sum - fcst.tbats_sum) ) %>%
select( ends_with("id"), date, ends_with("_AE"), actuals_sum) %>%
summarise(ets_err = sum(ets_AE)/sum(actuals_sum),
tbats_MAPE = sum(tbats_AE)/sum(actuals_sum),
naive_MAPE = sum(naive_AE)/sum(actuals_sum),
snaive_MAPE = sum(snaive_AE)/sum(actuals_sum),
mean_MAPE = sum(mean_AE)/sum(actuals_sum))
}

Now I apply the function made above at the different levels and groupings we can make in this hierarchy. Looking at MAPE approximations below we can see tbats has the best results *in aggregate* at all levels. It looks like the naive model, which forecasts that the next value will be the same as the last, does the best out of the control models but still not as good as the tbats model or the ets model. We can also notice how MAPE generally increases as we get into more granular levels below total sales. This is because of the nature of MAPE as a performance metric and is a normal trend, which I may add another article for later. If we were to go more granular and forecast at the SKU level, for example, you could expect the MAPE to become even higher. However, a higher MAPE does not mean the forecast is not useful and also presents an opportunity to use advanced methods to cut the MAPE down.

`results <- cbind(rbind(cbind(MAPE_performance(train_tbl_ts2)), cbind(MAPE_performance(train_tbl_ts2, state_id)), cbind(MAPE_performance(train_tbl_ts2, dept_id)), cbind(MAPE_performance(train_tbl_ts2, dept_id, state_id))), c("Total", "state_id", "dept_id", "dept_id and state_id")) colnames(results)[6] <- "Level"`

`kable(results, digits = 3)`

ets_err | tbats_MAPE | naive_MAPE | snaive_MAPE | mean_MAPE | Level |
---|---|---|---|---|---|

0.042 | 0.024 | 0.045 | 0.049 | 0.110 | Total |

0.046 | 0.031 | 0.050 | 0.056 | 0.110 | state_id |

0.058 | 0.045 | 0.062 | 0.100 | 0.121 | dept_id |

0.064 | 0.054 | 0.067 | 0.103 | 0.135 | dept_id and state_id |

Above, we looked at MAPE in aggregate at different levels. However, just because a model performs the best in aggregate does not mean it performs best for every time series at that level. To explore what this means we can look at some of the different time series below.

**NOTE** In all of the plots below, the actual sales are *green*, tbats results are *dark red*, ets is *lighter red*, and naive forecasts are *blue*. Sometimes ets results are covered by naive model results because they forecasted the same or very similar results at that point in time.

Looking at some of the sub-series in dept_id’s, we can look at the error metrics in absolute error and mean absolute percent error that sometimes the naive model is better than the tbats model, like HOBBIES_1 for example. So we are already finding that we could potentially improve the overall results even further if we could use more advanced methods to do model selection at individual time series within the hierarchy. However, we do see that most of the time tbats, the darker red line, does provider lower error metrics, supporting our MAPE table above.

`PLOT_level <-`

function(tb_ts22,...) { tb_ts22 %>% group_by(!!!enquos(...), date) %>% summarise(across(where(is.numeric), sum)) %>% ungroup() %>% mutate(ets_AE = abs(actuals_sum - fcst.ets_sum), snaive_AE = abs(actuals_sum - fcst.snaive_sum), naive_AE = abs(actuals_sum - fcst.naive_sum), mean_AE = abs(actuals_sum - fcst.mean_sum), tbats_AE = abs(actuals_sum - fcst.tbats_sum), ets_err = (ets_AE)/(actuals_sum), ets_MAPE = (ets_AE)/(actuals_sum), tbats_MAPE = (tbats_AE)/(actuals_sum), naive_MAPE = (naive_AE)/(actuals_sum), snaive_MAPE = (snaive_AE)/(actuals_sum), mean_MAPE = (mean_AE)/(actuals_sum), date = as.Date(date)) } plot1 <- PLOT_level(train_tbl_ts2, dept_id ) %>% filter(str_detect(dept_id ,'HO')) %>% ggplot(aes(x=date)) + geom_line(aes(y = fcst.tbats_sum), color = "darkred") + geom_line(aes(y = fcst.ets_sum), color = "red") + geom_line(aes(y = actuals_sum), color = "green") + geom_line(aes(y = fcst.naive_sum), color = "blue") +facet_grid(dept_id~.)+ labs(y = "Actual Sales and Forecasts") plot2 <- PLOT_level(train_tbl_ts2, dept_id ) %>% filter(str_detect(dept_id ,'HO')) %>% ggplot(aes(x=date)) + geom_line(aes(y = tbats_AE), color = "darkred") + geom_line(aes(y = ets_AE), color = "red") + geom_line(aes(y = naive_AE), color = "blue") +facet_grid(dept_id~.) + labs(y = "Absolute Error") plot3 <- PLOT_level(train_tbl_ts2, dept_id ) %>% filter(str_detect(dept_id ,'HO')) %>% ggplot(aes(x=date)) + geom_line(aes(y = tbats_MAPE), color = "darkred") + geom_line(aes(y = ets_MAPE), color = "red") + geom_line(aes(y = naive_MAPE), color = "blue") +facet_grid(dept_id~.) + labs(y = "Mean Absolute Percent Error")`plot1+plot2+plot3`

Looking in the food dept_id’s we again see that tbats on average appears to perform best, although there are instances where it does not.

`plot1 <- PLOT_level(train_tbl_ts2, dept_id ) %>% filter(str_detect(dept_id ,'FO')) %>% ggplot(aes(x=date)) + geom_line(aes(y = fcst.tbats_sum), color = "darkred") + geom_line(aes(y = fcst.ets_sum), color = "red") + geom_line(aes(y = actuals_sum), color = "green") + geom_line(aes(y = fcst.naive_sum), color = "blue") +facet_grid(dept_id~.)+ labs(y = "Actual Sales and Forecasts") plot2 <- PLOT_level(train_tbl_ts2, dept_id ) %>% filter(str_detect(dept_id ,'FO')) %>% ggplot(aes(x=date)) + geom_line(aes(y = tbats_AE), color = "darkred") + geom_line(aes(y = ets_AE), color = "red") + geom_line(aes(y = naive_AE), color = "blue") +facet_grid(dept_id~.) + labs(y = "Absolute Error") plot3 <- PLOT_level(train_tbl_ts2, dept_id ) %>% filter(str_detect(dept_id ,'FO')) %>% ggplot(aes(x=date)) + geom_line(aes(y = tbats_MAPE), color = "darkred") + geom_line(aes(y = ets_MAPE), color = "red") + geom_line(aes(y = naive_MAPE), color = "blue") +facet_grid(dept_id~.) + labs(y = "Mean Absolute Percent Error")`

`plot1+plot2+plot3`

Looking at state_id’s we see tbats outperforms in TX and CA but not really in WI. This could be another opportunity for further model improvement.

`plot1 <- PLOT_level(train_tbl_ts2, state_id ) %>%`

#filter(str_detect(dept_id ,'FO')) %>%ggplot(aes(x=date)) + geom_line(aes(y = fcst.tbats_sum), color = "darkred") + geom_line(aes(y = fcst.ets_sum), color = "red") + geom_line(aes(y = actuals_sum), color = "green") + geom_line(aes(y = fcst.naive_sum), color = "blue") +facet_grid(state_id~.)+ labs(y = "Actual Sales and Forecasts") plot2 <- PLOT_level(train_tbl_ts2, state_id ) %>%#filter(str_detect(dept_id ,'FO')) %>%ggplot(aes(x=date)) + geom_line(aes(y = tbats_AE), color = "darkred") + geom_line(aes(y = ets_AE), color = "red") + geom_line(aes(y = naive_AE), color = "blue") +facet_grid(state_id~.) + labs(y = "Absolute Error") plot3 <- PLOT_level(train_tbl_ts2, state_id ) %>%#filter(str_detect(dept_id ,'HO')) %>%ggplot(aes(x=date)) + geom_line(aes(y = tbats_MAPE), color = "darkred") + geom_line(aes(y = ets_MAPE), color = "red") + geom_line(aes(y = naive_MAPE), color = "blue") +facet_grid(state_id~.) + labs(y = "Mean Absolute Percent Error")`plot1+plot2+plot3`

Now as an example of looking at the most granular level that I have forecasted, we look at just a few examples of state_id crossed with dept_id. Again, we see tbats mostly outperforms, but that more sophisticated methods such as a more dynamic model selection, ensembling, and some other tricks, which will be discussed in a later blog, could create a better model.

`plot1 <- PLOT_level(train_tbl_ts2, dept_id, state_id ) %>% filter(str_detect(dept_id ,'FO')) %>% filter(str_detect(state_id ,'CA')) %>% ggplot(aes(x=date)) + geom_line(aes(y = fcst.tbats_sum), color = "darkred") + geom_line(aes(y = fcst.ets_sum), color = "red") + geom_line(aes(y = actuals_sum), color = "green") + geom_line(aes(y = fcst.naive_sum), color = "blue") +facet_grid(state_id+dept_id~.)+ labs(y = "Actual Sales and Forecasts") plot2 <- PLOT_level(train_tbl_ts2, dept_id, state_id ) %>% filter(str_detect(dept_id ,'FO')) %>% filter(str_detect(state_id ,'CA')) %>% ggplot(aes(x=date)) + geom_line(aes(y = tbats_AE), color = "darkred") + geom_line(aes(y = ets_AE), color = "red") + geom_line(aes(y = naive_AE), color = "blue") +facet_grid(state_id+dept_id~.) + labs(y = "Absolute Error") plot3 <- PLOT_level(train_tbl_ts2, dept_id, state_id ) %>% filter(str_detect(dept_id ,'FO')) %>% filter(str_detect(state_id ,'CA')) %>% ggplot(aes(x=date)) + geom_line(aes(y = tbats_MAPE), color = "darkred") + geom_line(aes(y = ets_MAPE), color = "red") + geom_line(aes(y = naive_MAPE), color = "blue") +facet_grid(state_id+dept_id~.) + labs(y = "Mean Absolute Percent Error")`

`plot1+plot2+plot3`