<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/using-a-genetic-algorithm-for-the-hyperparameter-optimization-of-a-sarima-model/"></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-16-rgenoud_arima/"> 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://keiwan.itch.io/evolution" target="_blank"><br>
<img width="400" title="Nietzsche's Übermensch" src="https://i2.wp.com/www.brodrigues.co/blog/2018-11-16-rgenoud_arima/img/tap-walker.gif?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. In my <a rel="nofollow" href="https://www.brodrigues.co/blog/2018-11-15-tidy_gridsearch/" target="_blank">last blog post</a><br>
I showed how to perform a grid search the “tidy” way. As an example, I looked for the right<br>
hyperparameters of a SARIMA model. However, the goal of the post was not hyperparameter optimization<br>
per se, so I did not bother with tuning the hyperparameters on a validation set, and used the test<br>
set for both validation of the hyperparameters and testing the forecast. Of course, this is not great<br>
because doing this might lead to overfitting the hyperparameters to the test set. So in this blog post<br>
I split my data into trainig, validation and testing sets and use a genetic algorithm to look<br>
for the hyperparameters. Again, this is not the most optimal way to go about this problem, since<br>
the <code>{forecast}</code> package contains the very useful <code>auto.arima()</code> function. I just wanted to see<br>
what kind of solution a genetic algorithm would return, and also try different cost functions.<br>
If you’re interested, read on!</p>
</div>
<div>
<h2>Setup</h2>
<p>Let’s first load some libraries and define some helper functions (the helper functions were explained<br>
in the previous blog posts):</p>
<pre><code>library(tidyverse)
library(forecast)
library(rgenoud)
library(parallel)
library(lubridate)
library(furrr)
library(tsibble)
library(brotools)

ihs &lt;- function(x){
    log(x + sqrt(x**2 + 1))
}

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>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 train set, a validation set and a test set:</p>
<pre><code>avia_clean_train &lt;- avia_clean_monthly %>%
    select(date, passengers) %>%
    filter(year(date) &lt; 2013) %>%
    group_by(date) %>%
    summarise(total_passengers = sum(passengers)) %>%
    pull(total_passengers) %>%
    ts(., frequency = 12, start = c(2005, 1))

avia_clean_validation &lt;- avia_clean_monthly %>%
    select(date, passengers) %>%
    filter(between(year(date), 2013, 2016)) %>%
    group_by(date) %>%
    summarise(total_passengers = sum(passengers)) %>%
    pull(total_passengers) %>%
    ts(., frequency = 12, start = c(2013, 1))

avia_clean_test &lt;- avia_clean_monthly %>%
    select(date, passengers) %>%
    filter(year(date) >= 2016) %>%
    group_by(date) %>%
    summarise(total_passengers = sum(passengers)) %>%
    pull(total_passengers) %>%
    ts(., frequency = 12, start = c(2016, 1))

logged_test_data &lt;- ihs(avia_clean_test)

logged_validation_data &lt;- ihs(avia_clean_validation)

logged_train_data &lt;- ihs(avia_clean_train)</code></pre>
<p>I will train the models on data from 2005 to 2012, look for the hyperparameters on data from 2013<br>
to 2016 and test the accuracy on data from 2016 to March 2018. For this kind of exercise, the ideal<br>
situation would be to perform cross-validation. Doing this with time-series data is not obvious<br>
because of the autocorrelation between observations, which would be broken by sampling independently<br>
which is required by CV. Also, if for example you do leave-one-out CV,<br>
you would end up trying to predict a point in, say, 2017, with data<br>
from 2018, which does not make sense. So you should be careful about that. <code>{forecast}</code> is able<br>
to perform <a rel="nofollow" href="https://robjhyndman.com/hyndsight/tscv/" target="_blank">CV for time series</a> and <code>scikit-learn</code>, the<br>
Python package, is able to perform<br>
<a rel="nofollow" href="https://scikit-learn.org/stable/modules/cross_validation.html#time-series-split" target="_blank">cross-validation of time series data</a><br>
too. I will not do it in this blog post and simply focus on the genetic algorithm part.</p>
<p>Let’s start by defining the cost function to minimize. I’ll try several, in the first one I will<br>
minimize the RMSE:</p>
<pre><code>cost_function_rmse &lt;- function(param, train_data, validation_data, forecast_periods){
    order &lt;- param[1:3]
    season &lt;- c(param[4:6], 12)
    model &lt;- purrr::possibly(arima, otherwise = NULL)(x = train_data, order = order, 
                                                      seasonal = season,
                                                      method = &quot;ML&quot;)
    if(is.null(model)){
        return(9999999)
    } else {
      forecast_model &lt;- forecast::forecast(model, h = forecast_periods)
      point_forecast &lt;- forecast_model$mean
      sqrt(mean(point_forecast - validation_data) ** 2)
    }
}</code></pre>
<p>If <code>arima()</code> is not able to estimate a model for the given parameters, I force it to return <code>NULL</code>,<br>
and in that case force the cost function to return a very high cost. If a model was successfully estimated,<br>
then I compute the RMSE.</p>
<p>Let’s also take a look at what <code>auto.arima()</code> says:</p>
<pre><code>starting_model &lt;- auto.arima(logged_train_data)
summary(starting_model)</code></pre>
<pre><code>## Series: logged_train_data 
## ARIMA(1,0,2)(2,1,0)[12] 
## 
## Coefficients:
##          ar1      ma1     ma2     sar1     sar2
##       0.9754  -0.7872  0.2091  -0.7285  -0.4413
## s.e.  0.0261   0.1228  0.1213   0.1063   0.1150
## 
## sigma^2 estimated as 0.004514:  log likelihood=105.61
## AIC=-199.22   AICc=-198.13   BIC=-184.64
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE
## Training set 0.008398036 0.06095102 0.03882593 0.07009285 0.3339574
##                   MASE       ACF1
## Training set 0.4425794 0.02073886</code></pre>
<p>Let’s compute the cost at this vector of parameters:</p>
<pre><code>cost_function_rmse(c(1, 0, 2, 2, 1, 0),
              train_data = logged_train_data,
              validation_data = logged_validation_data,
              forecast_periods = 65)</code></pre>
<pre><code>## [1] 0.1731473</code></pre>
<p>Ok, now let’s start with optimizing the hyperparameters. Let’s help the genetic algorithm a little<br>
bit by defining where it should perform the search:</p>
<pre><code>domains &lt;- matrix(c(0, 3, 0, 2, 0, 3, 0, 3, 0, 2, 0, 3), byrow = TRUE, ncol = 2)</code></pre>
<p>This matrix constraints the first parameter to lie between 0 and 3, the second one between 0 and 2,<br>
and so on.</p>
<p>Let’s call the <code>genoud()</code> function from the <code>{rgenoud}</code> package, and use 8 cores:</p>
<pre><code>cl &lt;- makePSOCKcluster(8)
clusterExport(cl, c('logged_train_data', 'logged_validation_data'))

tic &lt;- Sys.time()

auto_arima_rmse &lt;- genoud(cost_function_rmse,
                     nvars = 6,
                     data.type.int = TRUE,
                     starting.values = c(1, 0, 2, 2, 1, 0), # &lt;- from auto.arima
                     Domains = domains,
                     cluster = cl,
                     train_data = logged_train_data,
                     validation_data = logged_validation_data,
                     forecast_periods = length(logged_validation_data),
                     hard.generation.limit = TRUE)
toc_rmse &lt;- Sys.time() - tic</code></pre>
<p><code>makePSOCKcluster()</code> is a function from the <code>{parallel}</code> package. I must also <em>export</em> the global<br>
variables <code>logged_train_data</code> or <code>logged_validation_data</code>. If I don’t do that, the workers called<br>
by <code>genoud()</code> will not <em>know</em> about these variables and an error will be returned. The option<br>
<code>data.type.int = TRUE</code> force the algorithm to look only for integers, and <code>hard.generation.limit = TRUE</code><br>
forces the algorithm to stop after 100 generations.</p>
<p>The process took 7 minutes, which is faster than doing the grid search.<br>
What was the solution found?</p>
<pre><code>auto_arima_rmse</code></pre>
<pre><code>## $value
## [1] 0.0001863039
## 
## $par
## [1] 3 2 1 1 2 1
## 
## $gradients
## [1] NA NA NA NA NA NA
## 
## $generations
## [1] 11
## 
## $peakgeneration
## [1] 1
## 
## $popsize
## [1] 1000
## 
## $operators
## [1] 122 125 125 125 125 126 125 126   0</code></pre>
<p>Let’s train the model using the <code>arima()</code> function at these parameters:</p>
<pre><code>best_model_rmse &lt;- arima(logged_train_data, order = auto_arima_rmse$par[1:3], 
                         season = list(order = auto_arima_rmse$par[4:6], period = 12),
                         method = &quot;ML&quot;)

summary(best_model_rmse)</code></pre>
<pre><code>## 
## Call:
## arima(x = logged_train_data, order = auto_arima_rmse$par[1:3], seasonal = list(order = auto_arima_rmse$par[4:6], 
##     period = 12), method = &quot;ML&quot;)
## 
## Coefficients:
##           ar1      ar2      ar3      ma1     sar1     sma1
##       -0.6999  -0.4541  -0.0476  -0.9454  -0.4996  -0.9846
## s.e.   0.1421   0.1612   0.1405   0.1554   0.1140   0.2193
## 
## sigma^2 estimated as 0.006247:  log likelihood = 57.34,  aic = -100.67
## 
## Training set error measures:
##                         ME       RMSE        MAE          MPE      MAPE
## Training set -0.0006142355 0.06759545 0.04198561 -0.005408262 0.3600483
##                   MASE         ACF1
## Training set 0.4386693 -0.008298546</code></pre>
<p>Let’s extract the forecasts:</p>
<pre><code>best_model_rmse_forecast &lt;- forecast::forecast(best_model_rmse, h = 65)

best_model_rmse_forecast &lt;- to_tibble(best_model_rmse_forecast)</code></pre>
<pre><code>## Joining, by = &quot;date&quot;
## Joining, by = &quot;date&quot;</code></pre>
<pre><code>starting_model_forecast &lt;- forecast(starting_model, h = 65)

starting_model_forecast &lt;- to_tibble(starting_model_forecast)</code></pre>
<pre><code>## Joining, by = &quot;date&quot;
## Joining, by = &quot;date&quot;</code></pre>
<p>and plot the forecast to see how it looks:</p>
<pre><code>avia_clean_monthly %>%
    group_by(date) %>%
    summarise(total = sum(passengers)) %>%
    mutate(total_ihs = ihs(total)) %>%
    ggplot() +
    ggtitle(&quot;Minimization of RMSE&quot;) +
    geom_line(aes(y = total_ihs, x = date), colour = &quot;#82518c&quot;) +
    scale_x_date(date_breaks = &quot;1 year&quot;, date_labels = &quot;%m-%Y&quot;) +
    geom_ribbon(data = best_model_rmse_forecast, aes(x = date, ymin = lower95, ymax = upper95),
                fill = &quot;#666018&quot;, alpha = 0.2) +
    geom_line(data = best_model_rmse_forecast, aes(x = date, y = point_estimate), 
              linetype = 2, colour = &quot;#8e9d98&quot;) +
    geom_ribbon(data = starting_model_forecast, aes(x = date, ymin = lower95, ymax = upper95),
                fill = &quot;#98431e&quot;, alpha = 0.2) +
    geom_line(data = starting_model_forecast, aes(x = date, y = point_estimate), 
              linetype = 2, colour = &quot;#a53031&quot;) +
    theme_blog()</code></pre>
<p><img src="https://i0.wp.com/www.brodrigues.co/blog/2018-11-16-rgenoud_arima/blog/2018-11-16-rgenoud_arima_files/figure-html/unnamed-chunk-14-1.png?w=450" data-recalc-dims="1"></p>
<p>The yellowish line and confidence intervals come from minimizing the genetic algorithm, and the<br>
redish from <code>auto.arima()</code>. Interesting; the point estimate is very precise, but the confidence<br>
intervals are very wide. Low bias, high variance.</p>
<p>Now, let’s try with another cost function, where I minimize the BIC, similar to the <code>auto.arima()</code> function:</p>
<pre><code>cost_function_bic &lt;- function(param, train_data, validation_data, forecast_periods){
    order &lt;- param[1:3]
    season &lt;- c(param[4:6], 12)
    model &lt;- purrr::possibly(arima, otherwise = NULL)(x = train_data, order = order, 
                                                      seasonal = season,
                                                      method = &quot;ML&quot;)
    if(is.null(model)){
        return(9999999)
    } else {
        BIC(model)
    }
}</code></pre>
<p>Let’s take a look at the cost at the parameter values returned by <code>auto.arima()</code>:</p>
<pre><code>cost_function_bic(c(1, 0, 2, 2, 1, 0),
              train_data = logged_train_data,
              validation_data = logged_validation_data,
              forecast_periods = 65)</code></pre>
<pre><code>## [1] -184.6397</code></pre>
<p>Let the genetic algorithm run again:</p>
<pre><code>cl &lt;- makePSOCKcluster(8)
clusterExport(cl, c('logged_train_data', 'logged_validation_data'))

tic &lt;- Sys.time()

auto_arima_bic &lt;- genoud(cost_function_bic,
                     nvars = 6,
                     data.type.int = TRUE,
                     starting.values = c(1, 0, 2, 2, 1, 0), # &lt;- from auto.arima
                     Domains = domains,
                     cluster = cl,
                     train_data = logged_train_data,
                     validation_data = logged_validation_data,
                     forecast_periods = length(logged_validation_data),
                     hard.generation.limit = TRUE)
toc_bic &lt;- Sys.time() - tic</code></pre>
<p>This time, it took 6 minutes, a bit slower than before. Let’s take a look at the solution:</p>
<pre><code>auto_arima_bic</code></pre>
<pre><code>## $value
## [1] -201.0656
## 
## $par
## [1] 0 1 1 1 0 1
## 
## $gradients
## [1] NA NA NA NA NA NA
## 
## $generations
## [1] 12
## 
## $peakgeneration
## [1] 1
## 
## $popsize
## [1] 1000
## 
## $operators
## [1] 122 125 125 125 125 126 125 126   0</code></pre>
<p>Let’s train the model at these parameters:</p>
<pre><code>best_model_bic &lt;- arima(logged_train_data, order = auto_arima_bic$par[1:3], 
                        season = list(order = auto_arima_bic$par[4:6], period = 12),
                        method = &quot;ML&quot;)

summary(best_model_bic)</code></pre>
<pre><code>## 
## Call:
## arima(x = logged_train_data, order = auto_arima_bic$par[1:3], seasonal = list(order = auto_arima_bic$par[4:6], 
##     period = 12), method = &quot;ML&quot;)
## 
## Coefficients:
##           ma1    sar1    sma1
##       -0.6225  0.9968  -0.832
## s.e.   0.0835  0.0075   0.187
## 
## sigma^2 estimated as 0.004145:  log likelihood = 109.64,  aic = -211.28
## 
## Training set error measures:
##                       ME       RMSE        MAE        MPE      MAPE
## Training set 0.003710982 0.06405303 0.04358164 0.02873561 0.3753513
##                   MASE        ACF1
## Training set 0.4553447 -0.03450603</code></pre>
<p>And let’s plot the results:</p>
<pre><code>best_model_bic_forecast &lt;- forecast::forecast(best_model_bic, h = 65)

best_model_bic_forecast &lt;- to_tibble(best_model_bic_forecast)</code></pre>
<pre><code>## Joining, by = &quot;date&quot;
## Joining, by = &quot;date&quot;</code></pre>
<pre><code>avia_clean_monthly %>%
    group_by(date) %>%
    summarise(total = sum(passengers)) %>%
    mutate(total_ihs = ihs(total)) %>%
    ggplot() +
    ggtitle(&quot;Minimization of BIC&quot;) +
    geom_line(aes(y = total_ihs, x = date), colour = &quot;#82518c&quot;) +
    scale_x_date(date_breaks = &quot;1 year&quot;, date_labels = &quot;%m-%Y&quot;) +
    geom_ribbon(data = best_model_bic_forecast, aes(x = date, ymin = lower95, ymax = upper95),
                fill = &quot;#5160a0&quot;, alpha = 0.2) +
    geom_line(data = best_model_bic_forecast, aes(x = date, y = point_estimate), 
              linetype = 2, colour = &quot;#208480&quot;) +
    geom_ribbon(data = starting_model_forecast, aes(x = date, ymin = lower95, ymax = upper95),
                fill = &quot;#98431e&quot;, alpha = 0.2) +
    geom_line(data = starting_model_forecast, aes(x = date, y = point_estimate), 
              linetype = 2, colour = &quot;#a53031&quot;) +
    theme_blog()</code></pre>
<p><img src="https://i1.wp.com/www.brodrigues.co/blog/2018-11-16-rgenoud_arima/blog/2018-11-16-rgenoud_arima_files/figure-html/unnamed-chunk-21-1.png?w=450" data-recalc-dims="1"></p>
<p>The solutions are very close, both in terms of point estimates and confidence intervals. Bias<br>
increased, but variance lowered… This gives me an idea! What if I minimize the RMSE, while<br>
keeping the number of parameters low, as a kind of regularization? This is somewhat what minimising<br>
BIC does, but let’s try to do it a more “naive” approach:</p>
<pre><code>cost_function_rmse_low_k &lt;- function(param, train_data, validation_data, forecast_periods, max.order){
    order &lt;- param[1:3]
    season &lt;- c(param[4:6], 12)
    if(param[1] + param[3] + param[4] + param[6] > max.order){
        return(9999999)
    } else {
        model &lt;- purrr::possibly(arima, otherwise = NULL)(x = train_data, 
                                                          order = order, 
                                                          seasonal = season,
                                                          method = &quot;ML&quot;)
    }
    if(is.null(model)){
        return(9999999)
    } else {
        forecast_model &lt;- forecast::forecast(model, h = forecast_periods)
        point_forecast &lt;- forecast_model$mean
        sqrt(mean(point_forecast - validation_data) ** 2)
    }
}</code></pre>
<p>This is also similar to what <code>auto.arima()</code> does; by default, the <code>max.order</code> argument in <code>auto.arima()</code><br>
is set to 5, and is the sum of <code>p + q + P + Q</code>. So I’ll try something similar.</p>
<p>Let’s take a look at the cost at the parameter values returned by <code>auto.arima()</code>:</p>
<pre><code>cost_function_rmse_low_k(c(1, 0, 2, 2, 1, 0),
              train_data = logged_train_data,
              validation_data = logged_validation_data,
              forecast_periods = 65,
              max.order = 5)</code></pre>
<pre><code>## [1] 0.1731473</code></pre>
<p>Let’s see what will happen:</p>
<pre><code>cl &lt;- makePSOCKcluster(8)
clusterExport(cl, c('logged_train_data', 'logged_validation_data'))

tic &lt;- Sys.time()

auto_arima_rmse_low_k &lt;- genoud(cost_function_rmse_low_k,
                         nvars = 6,
                         data.type.int = TRUE,
                         starting.values = c(1, 0, 2, 2, 1, 0), # &lt;- from auto.arima
                         max.order = 5,
                         Domains = domains,
                         cluster = cl,
                         train_data = logged_train_data,
                         validation_data = logged_validation_data,
                         forecast_periods = length(logged_validation_data),
                         hard.generation.limit = TRUE)
toc_rmse_low_k &lt;- Sys.time() - tic</code></pre>
<p>It took 1 minute to train this one, quite fast! Let’s take a look:</p>
<pre><code>auto_arima_rmse_low_k</code></pre>
<pre><code>## $value
## [1] 0.002503478
## 
## $par
## [1] 1 2 0 3 1 0
## 
## $gradients
## [1] NA NA NA NA NA NA
## 
## $generations
## [1] 11
## 
## $peakgeneration
## [1] 1
## 
## $popsize
## [1] 1000
## 
## $operators
## [1] 122 125 125 125 125 126 125 126   0</code></pre>
<p>And let’s plot it:</p>
<pre><code>best_model_rmse_low_k &lt;- arima(logged_train_data, order = auto_arima_rmse_low_k$par[1:3], 
                               season = list(order = auto_arima_rmse_low_k$par[4:6], period = 12),
                               method = &quot;ML&quot;)

summary(best_model_rmse_low_k)</code></pre>
<pre><code>## 
## Call:
## arima(x = logged_train_data, order = auto_arima_rmse_low_k$par[1:3], seasonal = list(order = auto_arima_rmse_low_k$par[4:6], 
##     period = 12), method = &quot;ML&quot;)
## 
## Coefficients:
##           ar1     sar1     sar2     sar3
##       -0.6468  -0.7478  -0.5263  -0.1143
## s.e.   0.0846   0.1171   0.1473   0.1446
## 
## sigma^2 estimated as 0.01186:  log likelihood = 57.88,  aic = -105.76
## 
## Training set error measures:
##                        ME      RMSE        MAE         MPE      MAPE
## Training set 0.0005953302 0.1006917 0.06165919 0.003720452 0.5291736
##                   MASE       ACF1
## Training set 0.6442205 -0.3706693</code></pre>
<pre><code>best_model_rmse_low_k_forecast &lt;- forecast::forecast(best_model_rmse_low_k, h = 65)

best_model_rmse_low_k_forecast &lt;- to_tibble(best_model_rmse_low_k_forecast)</code></pre>
<pre><code>## Joining, by = &quot;date&quot;
## Joining, by = &quot;date&quot;</code></pre>
<pre><code>avia_clean_monthly %>%
    group_by(date) %>%
    summarise(total = sum(passengers)) %>%
    mutate(total_ihs = ihs(total)) %>%
    ggplot() +
    ggtitle(&quot;Minimization of RMSE + low k&quot;) +
    geom_line(aes(y = total_ihs, x = date), colour = &quot;#82518c&quot;) +
    scale_x_date(date_breaks = &quot;1 year&quot;, date_labels = &quot;%m-%Y&quot;) +
    geom_ribbon(data = best_model_rmse_low_k_forecast, aes(x = date, ymin = lower95, ymax = upper95),
                fill = &quot;#5160a0&quot;, alpha = 0.2) +
    geom_line(data = best_model_rmse_low_k_forecast, aes(x = date, y = point_estimate), 
              linetype = 2, colour = &quot;#208480&quot;) +
    geom_ribbon(data = starting_model_forecast, aes(x = date, ymin = lower95, ymax = upper95),
                fill = &quot;#98431e&quot;, alpha = 0.2) +
    geom_line(data = starting_model_forecast, aes(x = date, y = point_estimate), 
              linetype = 2, colour = &quot;#a53031&quot;) +
    theme_blog()</code></pre>
<p><img src="https://i1.wp.com/www.brodrigues.co/blog/2018-11-16-rgenoud_arima/blog/2018-11-16-rgenoud_arima_files/figure-html/unnamed-chunk-28-1.png?w=450" data-recalc-dims="1"></p>
<p>Looks like this was not the right strategy. There might be a better cost function than what I have<br>
tried, but looks like minimizing the BIC is the way to go.</p>
<p>Hope you enjoyed! If you found this blog post useful, you might want to follow<br>
me on <a rel="nofollow" href="https://www.twitter.com/brodriguesco" target="_blank">twitter</a> for blog post updates or<br>
<a rel="nofollow" href="https://www.buymeacoffee.com/brodriguesco" target="_blank">buy me an espresso</a>.</p>

<a rel="nofollow" href="https://www.buymeacoffee.com/brodriguesco" target="_blank"><img alt="Buy me an Espresso" src="https://www.buymeacoffee.com/assets/img/BMC-btn-logo.svg"><span>Buy me an Espresso</span></a>
</div>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://www.brodrigues.co/blog/2018-11-16-rgenoud_arima/"> Econometrics and Free Software</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a rel="nofollow" href="https://www.r-project.org/" title="The R Project for Statistical Computing">R</a> news and <a rel="nofollow" href="https://www.r-bloggers.com/search/tutorial" title="R tutorials">tutorials</a> on topics such as: <a rel="nofollow" href="https://www.r-bloggers.com/search/data%20science" title="Data science">Data science</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data" title="Big Data">Big Data, <a rel="nofollow" href="https://www.r-users.com/" title="R jobs">R jobs</a>, visualization (<a rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2" title="ggplot and ggplot2 tutorials">ggplot2</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/boxplot" title="Boxplots using lattice and ggplot2 tutorials">Boxplots</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/map" title="Maps and gis">maps</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/animation" title="Animation in R">animation</a>), programming (<a rel="nofollow" href="https://www.r-bloggers.com/search/RStudio" title="RStudio IDE for R">RStudio</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/sweave" title="Sweave and literate programming">Sweave</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX" title="LaTeX in R">LaTeX</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/SQL" title="SQL and databases">SQL</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/eclipse" title="Eclipse IDE for R">Eclipse</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/git" title="git and github, Version Control System">git</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/hadoop" title="Large data in R using Hadoop">hadoop</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping" title="Web Scraping of google, facebook, yahoo, twitter and more using R">Web Scraping</a>) statistics (<a rel="nofollow" href="https://www.r-bloggers.com/search/regression" title="Regressions and ANOVA analysis tutorials">regression</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/PCA" title="principal component analysis tutorial">PCA</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/time+series" title="Time series">time series</a>, <a rel="nofollow" href="https://www.r-bloggers.com/search/trading" title="finance trading">trading</a>) and more...
<img width="1" alt="" src="http://feeds.feedburner.com/~r/RBloggers/~4/MvsLLIe4sko" height="1"></a></div>