Searching for the optimal hyper-parameters of an ARIMA model in parallel： the tidy gridsearch approach

转载自：http://feedproxy.google.com/~r/RBloggers/~3/AFUW2oM19ms/

Econometrics and Free Software
发表于

Introduction In this blog post, I’ll use the data that I cleaned in a previousblog post, which you can downloadhere. 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 canbe all you need in a lot of situations, but perhaps you’d like to have a better model. I will showhere 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 to search for these hyper-parameters is the so-called Box-Jenkins method. You can read about ithere. This method was describedin a 1970 book, Time series analysis: Forecasting and control by Box and Jenkins. The methodrequires that you first prepare the data by logging it and differencing it, in order to make thetime series stationary. You then need to analyze ACF and PACF plots, in order to determine theright amount of lags… It take some time, but this method made sense in a time were computingpower 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 busydata 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 RMSEfor each model. I will then choose the model that has the lowest RMSE.

Setup

Let’s first load some libraries:

1
2
3
4
5
6
7
8
9
10
library(tidyverse)
library(forecast)
library(lubridate)
library(furrr)
library(tsibble)
library(brotools)

ihs <- function(x){
 log(x + sqrt(x**2 + 1))
}

Now, let’s load the data:

1
avia_clean_monthly <- read_csv("https://raw.githubusercontent.com/b-rodrigues/avia_par_lu/master/avia_clean_monthy.csv")

1
2
3
4
5
6
## Parsed with column specification:
## cols(
## destination = col_character(),
## date = col_date(format = ""),
## passengers = col_integer()
## )

Let’s split the data into a training set and into a testing set:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
avia_clean_train <- 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(2005, 1))

avia_clean_test <- 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 <- ihs(avia_clean_train)

logged_test_data <- ihs(avia_clean_test)

I also define a helper function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
to_tibble <- function(forecast_object){
 point_estimate <- forecast_object$mean %>%
 as_tsibble() %>%
 rename(point_estimate = value,
 date = index)

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

 lower <- forecast_object$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 usefullater, and is based on the code I already used in my previousblog post.

Now, let’s take a closer look at the arima() function:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
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("CSS-ML", "ML", "CSS"), n.cond,
 SSinit = c("Gardner1980", "Rossignol2011"),
 optim.method = "BFGS",
 optim.control = list(), kappa = 1e6)

The user is supposed to enter the hyper-parameters as two lists, one called order for p, d, qand one called seasonal for P, D, Q, S. So what we need is to define these lists:

1
2
3
4
5
order_list <- list("p" = seq(0, 3),
 "d" = seq(0, 2),
 "q" = seq(0, 3)) %>%
 cross() %>%
 map(lift(c))

I first start with order_list. This list has 3 elements, “p”, “d” and “q”. Each element is asequence from 0 to 3 (2 in the case of “d”). When I pass this list to purrr::cross() I get theproduct set of the starting list, so in this case a list of 434 = 48 elements. However, thislist looks pretty bad:

1
2
3
4
5
list("p" = seq(0, 3),
 "d" = seq(0, 2),
 "q" = seq(0, 3)) %>%
 cross() %>%
 head(3)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
## [[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
2
3
4
5
6
7
8
9
10
11
[[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 smallline 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 ofdefinition of a function from one type of input to another. The function whose input I am liftingis c(). So now, c() can take a list instead of a vector. Compare the following:

1
2
3
# The usual

c("a", "b")

1
1
## [1] "a" "b"

1
2
# Nothing happens
c(list("a", "b"))

1
2
3
4
5
## [[1]]
## [1] "a"
## 
## [[2]]
## [1] "b"

1
2
# Magic happens
lift(c)(list("a", "b"))

1
1
## [1] "a" "b"

So order_list is exactly what I wanted:

1
head(order_list)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
## [[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:

1
2
3
4
5
6
season_list <- list("P" = seq(0, 3),
 "D" = seq(0, 2),
 "Q" = seq(0, 3),
 "period" = 12) %>%
 cross() %>%
 map(lift(c))

I now coerce these two lists of vectors to tibbles:

1
2
3
orderdf <- tibble("order" = order_list)

seasondf <- tibble("season" = season_list)

And I can now finally create the grid of hyper-parameters:

1
2
3
4
5
hyper_parameters_df <- crossing(orderdf, seasondf)

nrows <- nrow(hyper_parameters_df)

head(hyper_parameters_df)

1
2
3
4
5
6
7
8
9
## # A tibble: 6 x 2
## order season 
## 
## 1 
## 2 
## 3 
## 4 
## 5 
## 6

The hyper_parameters_df data frame has 2304 rows, meaning, I will now estimate 2304models, and will do so in parallel. Let’s just take a quick look at the internals of hyper_parameters_df:

1
glimpse(hyper_parameters_df)

1
2
3
4
## Observations: 2,304
## Variables: 2
## $ order [<0, 0, 0>, <0, 0, 0>, <0, 0, 0>, <0, 0, 0>, <0, 0, 0>...
## $ season [<0, 0, 0, 12>, <1, 0, 0, 12>, <2, 0, 0, 12>, <3, 0, 0...

So in the order column, the vector 0, 0, 0 is repeated as many times as there are combinationsof P, D, Q, S for season. Same for all the other vectors of the order column.

Training the models

Because training these models might take some time, I will use the fantastic {furrr} packageby Davis Vaughan to train the arima() function in parallel.For this, I first define 8 workers:

1
plan(multiprocess, workers = 8)

And then I run the code:

1
2
3
4
5
6
7
tic <- Sys.time()
models_df <- hyper_parameters_df %>%
 mutate(models = future_map2(.x = order,
 .y = season,
 ~possibly(arima, otherwise = NULL)(x = logged_train_data,
 order = .x, seasonal = .y)))
running_time <- Sys.time() - tic

I use future_map2(), which is just like map2() but running in parallel.I add a new column to the data called models, which will contain the models trained over all thedifferent combinations of order and season. The models are trained on the logged_train_data.

Training the 2304 models took 18 minutes, which isplenty of time to browse the latest memes, but still quick enough that it justifies the whole approach.Let’s take a look at the models_df object:

1
1
head(models_df)

1
2
3
4
5
6
7
8
9
## # A tibble: 6 x 3
## order season models 
## 
## 1 
## 2 
## 3 
## 4 
## 5 
## 6

As you can see, the models column contains all the trained models. The model on the first row,was trained with the hyperparameters of row 1, and so on. But, our work is not over! We now needto find the best model. First, I add a new column to the tibble, which contains the forecast. Fromthe forecast, I extract the point estimate:

1
2
3
4
models_df %>%
 mutate(forecast = map(models, ~possibly(forecast, otherwise = NULL)(., h = 39))) %>%
 mutate(point_forecast = map(forecast, ~.$`mean`)) %>%
 ....

You have to be familiar with a forecast object to understand the last line: a forecast objectis a list with certain elements, the point estimates, the confidence intervals, and so on. To getthe point estimates, I have to extract the “mean” element from the list. Hence the weird ~.$mean.Then I need to add a new list-column, where each element is the vector of true values, meaning the datafrom 2015 to 2018. Because I have to add it as a list of size 2304, I do that with purrr::rerun():

1
rerun(5, c("a", "b", "c"))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
## [[1]]
## [1] "a" "b" "c"
## 
## [[2]]
## [1] "a" "b" "c"
## 
## [[3]]
## [1] "a" "b" "c"
## 
## [[4]]
## [1] "a" "b" "c"
## 
## [[5]]
## [1] "a" "b" "c"

It is then easy to compute the RMSE, which I add as a column to the original data:

1
2
3
4
... %>%
 mutate(true_value = rerun(nrows, logged_test_data)) %>%
 mutate(rmse = map2_dbl(point_forecast, true_value,
 ~sqrt(mean((.x - .y) ** 2))))

The whole workflow is here:

1
2
3
4
5
6
models_df <- models_df %>%
 mutate(forecast = map(models, ~possibly(forecast, otherwise = NULL)(., h = 39))) %>%
 mutate(point_forecast = map(forecast, ~.$`mean`)) %>%
 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:

1
1
head(models_df)

1
2
3
4
5
6
7
8
9
## # A tibble: 6 x 7
## order season models forecast point_forecast true_value rmse
## 
## 1 0.525
## 2 0.236
## 3 0.235
## 4 0.217
## 5 0.190
## 6 0.174

Now, I can finally select the best performing model. I select the model with minimum RMSE:

1
2
best_model <- models_df %>%
 filter(rmse == min(rmse, na.rm = TRUE))

And save the forecast into a new variable, as a tibble, using my to_tibble() function:

1
(best_model_forecast <- to_tibble(best_model$forecast[[1]]))

1
2
## Joining, by = "date"
## Joining, by = "date"

1
2
3
4
5
6
7
8
9
10
11
12
13
14
## # A tsibble: 39 x 6 [1M]
## date point_estimate upper80 upper95 lower80 lower95
## 
## 1 2015 Jan 11.9 12.1 12.1 11.8 11.7
## 2 2015 Feb 11.9 12.0 12.1 11.7 11.6
## 3 2015 Mar 12.1 12.3 12.3 11.9 11.9
## 4 2015 Apr 12.2 12.3 12.4 12.0 11.9
## 5 2015 May 12.2 12.4 12.5 12.1 12.0
## 6 2015 Jun 12.3 12.4 12.5 12.1 12.0
## 7 2015 Jul 12.2 12.3 12.4 12.0 11.9
## 8 2015 Aug 12.3 12.5 12.6 12.2 12.1
## 9 2015 Sep 12.3 12.5 12.6 12.2 12.1
## 10 2015 Oct 12.2 12.4 12.5 12.1 12.0
## # ... with 29 more rows

And now, I can plot it:

1
2
3
4
5
6
7
8
9
10
11
12
avia_clean_monthly %>%
 group_by(date) %>%
 summarise(total = sum(passengers)) %>%
 mutate(total_ihs = ihs(total)) %>%
 ggplot() +
 ggtitle("Logged data") +
 geom_line(aes(y = total_ihs, x = date), colour = "#82518c") +
 scale_x_date(date_breaks = "1 year", date_labels = "%m-%Y") +
 geom_ribbon(data = best_model_forecast, aes(x = date, ymin = lower95, ymax = upper95), 
 fill = "#666018", alpha = 0.2) +
 geom_line(data = best_model_forecast, aes(x = date, y = point_estimate), linetype = 2, colour = "#8e9d98") +
 theme_blog()

Compared to the previous blog post, thedotted line now seems to follow the true line even better!Now, I am not saying that you should always do a gridsearch whenever you have a problem like thisone. In the case of univariate time series, I am still doubtful that a gridsearch like this is reallynecessary. However, it makes for a good exercise and a good illustration of the power of {furrr}.

Hope you enjoyed! If you found this blog post useful, you might want to followme on twitter for blog post updates.

Related