Introduction
In this blog post, I’ll use the data that I cleaned in a previous
blog post, which you can download
here. If you want to follow along,
download the monthly data.
In the previous blog post, I used the auto.arima() function to very quickly get a “good-enough”
model to predict future monthly total passengers flying from LuxAirport. “Good-enough” models can
be all you need in a lot of situations, but perhaps you’d like to have a better model. I will show
here how you can get a better model by searching through a grid of hyper-parameters.
This blog post was partially inspired by: https://drsimonj.svbtle.com/grid-search-in-the-tidyverse
The problem
SARIMA models have a lot of hyper-parameters, 7 in total! Three trend hyper-parameters, p, d, q,
same as for an ARIMA model, and four seasonal hyper-parameters, P, D, Q, S. The traditional way t
o search for these hyper-parameters is the so-called Box-Jenkins method. You can read about it
here. This method was described
in a 1970 book, Time series analysis: Forecasting and control by Box and Jenkins. The method
requires that you first prepare the data by logging it and differencing it, in order to make the
time series stationary. You then need to analyze ACF and PACF plots, in order to determine the
right amount of lags… It take some time, but this method made sense in a time were computing
power was very expensive. Today, we can simply let our computer search through thousands of models,
check memes on the internet, and come back to the best fit. This blog post is for you, the busy
data scientist meme connoisseurs who cannot waste time with theory and other such useless time drains,
when there are literally thousands of new memes being created and shared every day. Every second counts.
To determine what model is best, I will do pseudo out-of-sample forecasting and compute the RMSE
for each model. I will then choose the model that has the lowest RMSE.
Setup
Let’s first load some libraries:
library(tidyverse)
library(forecast)
library(lubridate)
library(furrr)
library(tsibble)
library(brotools)
ihs %
group_by(date) %__%
summarise(total_passengers = sum(passengers)) %__%
pull(total_passengers) %__%
ts(., frequency = 12, start = c(2005, 1))
avia_clean_test %
select(date, passengers) %__%
filter(year(date) __= 2015) %__%
group_by(date) %__%
summarise(total_passengers = sum(passengers)) %__%
pull(total_passengers) %__%
ts(., frequency = 12, start = c(2015, 1))
logged_train_data %
rename(date = index,
upper80 = `80%`,
upper95 = `95%`)
lower %
as_tsibble() %__%
spread(key, value) %__%
rename(date = index,
lower80 = `80%`,
lower95 = `95%`)
reduce(list(point_estimate, upper, lower), full_join)
}
This function takes a forecast object as argument, and returns a nice tibble. This will be useful
later, and is based on the code I already used in my previous
blog post.
Now, let’s take a closer look at the arima() function:
ARIMA Modelling of Time Series
Description
Fit an ARIMA model to a univariate time series.
Usage
arima(x, order = c(0L, 0L, 0L),
seasonal = list(order = c(0L, 0L, 0L), period = NA),
xreg = NULL, include.mean = TRUE,
transform.pars = TRUE,
fixed = NULL, init = NULL,
method = c(&quot;CSS-ML&quot;, &quot;ML&quot;, &quot;CSS&quot;), n.cond,
SSinit = c(&quot;Gardner1980&quot;, &quot;Rossignol2011&quot;),
optim.method = &quot;BFGS&quot;,
optim.control = list(), kappa = 1e6)
The user is supposed to enter the hyper-parameters as two lists, one called order for p, d, q
and one called seasonal for P, D, Q, S. So what we need is to define these lists:
order_list %
cross() %__%
map(lift(c))
I first start with order_list. This list has 3 elements, “p”, “d” and “q”. Each element is a
sequence from 0 to 3 (2 in the case of “d”). When I pass this list to purrr::cross() I get the
product set of the starting list, so in this case a list of 4*3*4 = 48 elements. However, this
list looks pretty bad:
list(&quot;p&quot; = seq(0, 3),
&quot;d&quot; = seq(0, 2),
&quot;q&quot; = seq(0, 3)) %__%
cross() %__%
head(3)
## [[1]]
## [[1]]$p
## [1] 0
## ## [[1]]$d
## [1] 0
## ## [[1]]$q
## [1] 0
## ## ## [[2]]
## [[2]]$p
## [1] 1
## ## [[2]]$d
## [1] 0
## ## [[2]]$q
## [1] 0
## ## ## [[3]]
## [[3]]$p
## [1] 2
## ## [[3]]$d
## [1] 0
## ## [[3]]$q
## [1] 0
I would like to have something like this instead:
[[1]]
p d q 0 0 0 [[2]]
p d q 1 0 0 [[3]]
p d q 2 0 0 This is possible with the last line, map(lift(c)). There’s a lot going on in this very small
line of code. First of all, there’s map(). map() iterates over lists, and applies a function,
in this case lift(c). purrr::lift() is a very interesting function that lifts the domain of
definition of a function from one type of input to another. The function whose input I am lifting
is c(). So now, c() can take a list instead of a vector. Compare the following:
# The usual
c(&quot;a&quot;, &quot;b&quot;)
## [1] &quot;a&quot; &quot;b&quot;
# Nothing happens
c(list(&quot;a&quot;, &quot;b&quot;))
## [[1]]
## [1] &quot;a&quot;
## ## [[2]]
## [1] &quot;b&quot;
# Magic happens
lift(c)(list(&quot;a&quot;, &quot;b&quot;))
## [1] &quot;a&quot; &quot;b&quot;
So order_list is exactly what I wanted:
head(order_list)
## [[1]]
## p d q ## 0 0 0 ## ## [[2]]
## p d q ## 1 0 0 ## ## [[3]]
## p d q ## 2 0 0 ## ## [[4]]
## p d q ## 3 0 0 ## ## [[5]]
## p d q ## 0 1 0 ## ## [[6]]
## p d q ## 1 1 0
I do the same for season_list:
season_list %
cross() %__%
map(lift(c))
I now coerce these two lists of vectors to tibbles:
orderdf %
mutate(true_value = rerun(nrows, logged_test_data)) %__%
mutate(rmse = map2_dbl(point_forecast, true_value,
~sqrt(mean((.x - .y) ** 2))))
This is how models_df looks now:
head(models_df)
## # A tibble: 6 x 7
## order season models forecast point_forecast true_value rmse
## ## 1