Thanks for providing a reproducible example. It makes answering much easier. Here is some code to do what you want.
suppressWarnings(suppressMessages(library(fpp3)))
# Data
toy_data <- PBS %>%
filter(ATC1 == "A", ATC2 == "A01", Concession == "General") %>%
select(- Scripts)
train_data <- toy_data %>%
filter_index(~ "2005 Dec")
# Fit all models
fit <- train_data %>%
model(
# Model 1
`STL + ARIMA` = decomposition_model(
STL(Cost ~ trend(window = 21) + season(window=13), robust = TRUE),
ARIMA(season_adjust)),
# Model 2
ARIMA = ARIMA(Cost)
)
# Forecasts from all models
forecasts <- fit %>%
forecast(h = 36)
# Find best models using RMSE
bestrmse <- accuracy(forecasts, toy_data) %>%
group_by(Concession, Type, ATC1, ATC2) %>%
filter(RMSE == min(RMSE)) %>%
select(.model:ATC2)
# Keep best forecasts
bestfc <- forecasts %>%
right_join(bestrmse)
#> Joining, by = c("Concession", "Type", "ATC1", "ATC2", ".model")
# Modify mable to only keep the best models
bestfits <- fit %>%
pivot_longer(cols=`STL + ARIMA`:ARIMA, names_to = ".model", values_to = "fit") %>%
right_join(bestrmse) %>%
mutate(.model = "best") %>%
pivot_wider(Concession:ATC2, names_from = ".model", values_from = "fit") %>%
as_mable(key = c(Concession, Type, ATC1, ATC2), model=best)
#> Joining, by = c("Concession", "Type", "ATC1", "ATC2", ".model")
Created on 2020-10-28 by the reprex package (v0.3.0)