<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/searching-for-the-optimal-hyper-parameters-of-an-arima-model-in-parallel-the-tidy-gridsearch-approach/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://www.brodrigues.co/blog/2018-11-15-tidy_gridsearch/"> Econometrics and Free Software</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<div>
<p><a rel="nofollow" href="https://youtu.be/3NxM-AL18lU?t=33s" target="_blank"><br>
<img width="400" title="What a time to be alive" src="https://i1.wp.com/www.brodrigues.co/blog/2018-11-15-tidy_gridsearch/img/dank_memes.jpg?w=456" data-recalc-dims="1"></a></p>
</div>
<div>
<h2>Introduction</h2>
<p>In this blog post, I’ll use the data that I cleaned in a previous<br>
<a rel="nofollow" href="https://www.brodrigues.co/blog/2018-11-14-luxairport/" target="_blank">blog post</a>, which you can download<br>
<a rel="nofollow" href="https://github.com/b-rodrigues/avia_par_lu/tree/master" target="_blank">here</a>. If you want to follow along,<br>
download the monthly data.</p>
<p>In the previous blog post, I used the <code>auto.arima()</code> function to very quickly get a “good-enough”<br>
model to predict future monthly total passengers flying from LuxAirport. “Good-enough” models can<br>
be all you need in a lot of situations, but perhaps you’d like to have a better model. I will show<br>
here how you can get a better model by searching through a grid of hyper-parameters.</p>
<p>This blog post was partially inspired by: <a rel="nofollow" href="https://drsimonj.svbtle.com/grid-search-in-the-tidyverse" target="_blank">https://drsimonj.svbtle.com/grid-search-in-the-tidyverse</a></p>
</div>
<div>
<h2>The problem</h2>
<p>SARIMA models have a lot of hyper-parameters, 7 in total! Three trend hyper-parameters, <em>p, d, q</em>,<br>
same as for an ARIMA model, and four seasonal hyper-parameters, <em>P, D, Q, S</em>. The traditional way t<br>
o search for these hyper-parameters is the so-called Box-Jenkins method. You can read about it<br>
<a rel="nofollow" href="https://www.itl.nist.gov/div898/handbook/pmc/section4/pmc44a.htm" target="_blank">here</a>. This method was described<br>
in a 1970 book, <em>Time series analysis: Forecasting and control</em> by Box and Jenkins. The method<br>
requires that you first prepare the data by logging it and differencing it, in order to make the<br>
time series stationary. You then need to analyze ACF and PACF plots, in order to determine the<br>
right amount of lags… It take some time, but this method made sense in a time were computing<br>
power was very expensive. Today, we can simply let our computer search through thousands of models,<br>
check memes on the internet, and come back to the best fit. This blog post is for you, the busy<br>
data scientist meme connoisseurs who cannot waste time with theory and other such useless time drains,<br>
when there are literally thousands of new memes being created and shared every day. Every second counts.<br>
To determine what model is best, I will do pseudo out-of-sample forecasting and compute the RMSE<br>
for each model. I will then choose the model that has the lowest RMSE.</p>
</div>
<div>
<h2>Setup</h2>
<p>Let’s first load some libraries:</p>
<pre><code>library(tidyverse)
library(forecast)
library(lubridate)
library(furrr)
library(tsibble)
library(brotools)

ihs &lt;- function(x){
    log(x + sqrt(x**2 + 1))
}</code></pre>
<p>Now, let’s load the data:</p>
<pre><code>avia_clean_monthly &lt;- read_csv(&quot;https://raw.githubusercontent.com/b-rodrigues/avia_par_lu/master/avia_clean_monthy.csv&quot;)</code></pre>
<pre><code>## Parsed with column specification:
## cols(
##   destination = col_character(),
##   date = col_date(format = &quot;&quot;),
##   passengers = col_integer()
## )</code></pre>
<p>Let’s split the data into a training set and into a testing set:</p>
<pre><code>avia_clean_train &lt;- avia_clean_monthly %>%
    select(date, passengers) %>%
    filter(year(date) &lt; 2015) %>%
    group_by(date) %>%
    summarise(total_passengers = sum(passengers)) %>%
    pull(total_passengers) %>%
    ts(., frequency = 12, start = c(2005, 1))

avia_clean_test &lt;- avia_clean_monthly %>%
    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 &lt;- ihs(avia_clean_train)

logged_test_data &lt;- ihs(avia_clean_test)</code></pre>
<p>I also define a helper function:</p>
<pre><code>to_tibble &lt;- function(forecast_object){
    point_estimate &lt;- forecast_object$mean %>%
        as_tsibble() %>%
        rename(point_estimate = value,
               date = index)

    upper &lt;- forecast_object$upper %>%
        as_tsibble() %>%
        spread(key, value) %>%
        rename(date = index,
               upper80 = `80%`,
               upper95 = `95%`)

    lower &lt;- forecast_object$lower %>%
        as_tsibble() %>%
        spread(key, value) %>%
        rename(date = index,
               lower80 = `80%`,
               lower95 = `95%`)

    reduce(list(point_estimate, upper, lower), full_join)
}</code></pre>
<p>This function takes a <code>forecast</code> object as argument, and returns a nice tibble. This will be useful<br>
later, and is based on the code I already used in my previous<br>
<a rel="nofollow" href="https://www.brodrigues.co/blog/2018-11-14-luxairport/" target="_blank">blog post</a>.</p>
<p>Now, let’s take a closer look at the <code>arima()</code> function:</p>
<pre><code>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)</code></pre>
<p>The user is supposed to enter the hyper-parameters as two lists, one called <code>order</code> for <em>p, d, q</em><br>
and one called <code>seasonal</code> for <em>P, D, Q, S</em>. So what we need is to define these lists:</p>
<pre><code>order_list &lt;- list(&quot;p&quot; = seq(0, 3),
                   &quot;d&quot; = seq(0, 2),
                   &quot;q&quot; = seq(0, 3)) %>%
    cross() %>%
    map(lift(c))</code></pre>
<p>I first start with <code>order_list</code>. This list has 3 elements, “p”, “d” and “q”. Each element is a<br>
sequence from 0 to 3 (2 in the case of “d”). When I pass this list to <code>purrr::cross()</code> I get the<br>
product set of the starting list, so in this case a list of 4*3*4 = 48 elements. However, this<br>
list looks pretty bad:</p>
<pre><code>list(&quot;p&quot; = seq(0, 3),
     &quot;d&quot; = seq(0, 2),
     &quot;q&quot; = seq(0, 3)) %>%
    cross() %>%
    head(3)</code></pre>
<pre><code>## [[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</code></pre>
<p>I would like to have something like this instead:</p>
<pre><code>[[1]]
p d q 
0 0 0 

[[2]]
p d q 
1 0 0 

[[3]]
p d q 
2 0 0 </code></pre>
<p>This is possible with the last line, <code>map(lift(c))</code>. There’s a lot going on in this very small<br>
line of code. First of all, there’s <code>map()</code>. <code>map()</code> iterates over lists, and applies a function,<br>
in this case <code>lift(c)</code>. <code>purrr::lift()</code> is a very interesting function that lifts the domain of<br>
definition of a function from one type of input to another. The function whose input I am lifting<br>
is <code>c()</code>. So now, <code>c()</code> can take a list instead of a vector. Compare the following:</p>
<pre><code># The usual

c(&quot;a&quot;, &quot;b&quot;)</code></pre>
<pre><code>## [1] &quot;a&quot; &quot;b&quot;</code></pre>
<pre><code># Nothing happens
c(list(&quot;a&quot;, &quot;b&quot;))</code></pre>
<pre><code>## [[1]]
## [1] &quot;a&quot;
## 
## [[2]]
## [1] &quot;b&quot;</code></pre>
<pre><code># Magic happens
lift(c)(list(&quot;a&quot;, &quot;b&quot;))</code></pre>
<pre><code>## [1] &quot;a&quot; &quot;b&quot;</code></pre>
<p>So <code>order_list</code> is exactly what I wanted:</p>
<pre><code>head(order_list)</code></pre>
<pre><code>## [[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</code></pre>
<p>I do the same for <code>season_list</code>:</p>
<pre><code>season_list &lt;- list(&quot;P&quot; = seq(0, 3),
                    &quot;D&quot; = seq(0, 2),
                    &quot;Q&quot; = seq(0, 3),
                    &quot;period&quot; = 12)  %>%
    cross() %>%
    map(lift(c))</code></pre>
<p>I now coerce these two lists of vectors to tibbles:</p>
<pre><code>orderdf &lt;- tibble(&quot;order&quot; = order_list)

seasondf &lt;- tibble(&quot;season&quot; = season_list)</code></pre>
<p>And I can now finally create the grid of hyper-parameters:</p>
<pre><code>hyper_parameters_df &lt;- crossing(orderdf, seasondf)

nrows &lt;- nrow(hyper_parameters_df)

head(hyper_parameters_df)</code></pre>
<pre><code>## # A tibble: 6 x 2
##   order     season   
##   </code></pre></div>