Forecasting Time Series: Auto ARIMA, Ensemble Method and Artificial Neural Networks

R for Fun

My Code

# Clear workspace: 
rm(list = ls())

# Load some packages: 
library(nnfor)
library(fpp2)
library(tidyverse)
library(magrittr)

# Visualize AirPassengers Data: 

autoplot(AirPassengers) + 
  geom_line(color = "cyan") + 
  geom_point(color = "cyan") + 
  labs(x = NULL, y = NULL, 
       title = "Monthly Airline Passenger Numbers 1949-1960", 
       subtitle = "Data Source: The classic Box & Jenkins airline data.") + 
  theme(
    axis.line = element_blank(),  
    axis.text.x = element_text(color = "white", lineheight = 0.9),  
    axis.text.y = element_text(color = "white", lineheight = 0.9),  
    axis.ticks = element_line(color = "white", size  =  0.2),  
    axis.title.x = element_text(color = "white", margin = margin(0, 10, 0, 0)),  
    axis.title.y = element_text(color = "white", angle = 90, margin = margin(0, 10, 0, 0)),  
    axis.ticks.length = unit(0.3, "lines"),   
    legend.background = element_rect(color = NA, fill = " gray10"),  
    legend.key = element_rect(color = "white",  fill = " gray10"),  
    legend.key.size = unit(1.2, "lines"),  
    legend.key.height = NULL,  
    legend.key.width = NULL,      
    legend.text = element_text(color = "white"),  
    legend.title = element_text(face = "bold", hjust = 0, color = "white"),  
    legend.text.align = NULL,  
    legend.title.align = NULL,  
    legend.direction = "vertical",  
    legend.box = NULL, 
    panel.background = element_rect(fill = "gray10", color  =  NA),  
    panel.border = element_blank(),
    panel.grid.major = element_line(color = "grey35"),  
    panel.grid.minor = element_line(color = "grey20"),  
    panel.spacing = unit(0.5, "lines"),   
    strip.background = element_rect(fill = "grey30", color = "grey10"),  
    strip.text.x = element_text(color = "white"),  
    strip.text.y = element_text(color = "white", angle = -90),  
    plot.background = element_rect(color = "gray10", fill = "gray10"),  
    plot.title = element_text(color = "white", hjust = 0, lineheight = 1.25,
                              margin = margin(2, 2, 2, 2)),  
    plot.subtitle = element_text(color = "white", hjust = 0, margin = margin(2, 2, 2, 2)),  
    plot.caption = element_text(color = "white", hjust = 0),  
    plot.margin = unit(rep(1, 4), "lines"))

# Split data: 

train <- AirPassengers[1:132] %>% ts(frequency = 12)
test <- AirPassengers[133:144] %>% ts(frequency = 12)
h <- length(test)

#======================================================
#  MLP neural networks for Forecasting Time Series
#======================================================

# Prepare a series for train and forecasting: 
tt <- cbind(c(1:(length(train) + h), rep(0, h)))

# Fit a network with no differencing, no univariate lags, and fixed deterministic trend: 
my_mlp <- mlp(train, 
              hd = 10, 
              difforder = 0, 
              lags = 0, 
              xreg = tt, 
              xreg.lags = list(0), 
              xreg.keep = TRUE)

# Make predictions: 
pred <- forecast(my_mlp, h = h, xreg = tt)

# Data Frame for comparing: 

df_mlp <- data_frame(Actual = test %>% as.vector(), 
                     Predicted = pred$mean %>% round(0), 
                     Error = Predicted - Actual,
                     Error_Percent = round(Error / Actual, 2))

df_mlp %>% knitr::kable()

Actual	Predicted	Error	Error_Percent
417	392	-25	-0.06
391	370	-21	-0.05
419	434	15	0.04
461	414	-47	-0.10
472	447	-25	-0.05
535	506	-29	-0.05
622	578	-44	-0.07
606	590	-16	-0.03
508	487	-21	-0.04
461	427	-34	-0.07
390	383	-7	-0.02
432	426	-6	-0.01

# Function calculate some accuracy measures: 
get_accuracy_measures <- function(your_result_df) {
  
  act <- your_result_df %>% pull(Actual)
  pred <- your_result_df %>% pull(Predicted)
  err <- act - pred 
  per_err <- abs(err / act)
  
  # Mean Absolute Error (MAE): 
  mae <- err %>% abs() %>% mean()
  
  # Mean Squared Error (MSE): 
  mse <- mean(err^2)
  
  # Mean Absolute Percentage Error (MAPE): 
  mape <- 100*mean(per_err)
  
  # Return results: 
  return(data_frame(MAE = mae, MSE = mse, MAPE = mape, N = length(act)))
}


#==================================================================
#  Auto ARIMA Approach as proposed by Hyndman and Khandakar (2008)
#==================================================================

my_arima <- auto.arima(train)

# Use the model for forecasting: 
predicted_arima <- forecast(my_arima, h = h)$mean %>% as.vector()


# Data Frame for comparing: 
df_arima <- data_frame(Actual = test %>% as.vector(), 
                       Predicted = predicted_arima %>% round(0), 
                       Error = Predicted - Actual,
                       Error_Percent = round(Error / Actual, 2))

# Model Performance: 
df_arima %>% knitr::kable()

Actual	Predicted	Error	Error_Percent
417	424	7	0.02
391	407	16	0.04
419	471	52	0.12
461	461	0	0.00
472	485	13	0.03
535	537	2	0.00
622	613	-9	-0.01
606	624	18	0.03
508	528	20	0.04
461	472	11	0.02
390	427	37	0.09
432	470	38	0.09

#=======================================
#    My Solution: Ensemble Method 
#=======================================

my_ensemble_predict <- function(your_df, n_ahead) {
  
  N <- length(your_df)
  end1 <- N - n_ahead
  
  train <- your_df[1:end1] %>% ts(frequency = 12)
  test <- your_df[(end1 + 1):N] %>% ts(frequency = 12)
  h <- length(test)
  
  tt <- cbind(c(1:(length(train) + h), rep(0, h)))
  
  my_mlp <- mlp(train, 
                hd = 10, 
                difforder = 0, 
                lags = 0, 
                xreg = tt, 
                xreg.lags = list(0), 
                xreg.keep = TRUE)
  predicted_mlp <- forecast(my_mlp, h = h, xreg = tt)
  
  my_arima <- auto.arima(train)
  predicted_arima <- forecast(my_arima, h = h)$mean %>% as.vector()
  
  ensemble_predicted <- 0.5*(predicted_mlp$mean + predicted_arima)
  
  return(data_frame(Actual = test %>% as.vector(), 
                    Predicted = ensemble_predicted %>% round(0), 
                    Error = Predicted - Actual,
                    Error_Percent = round(Error / Actual, 2)))
} 


# Use the function: 
my_ensemble_predict(AirPassengers, 12) -> df_ensemble
df_ensemble %>% knitr::kable()

Actual	Predicted	Error	Error_Percent
417	408	-9	-0.02
391	390	-1	0.00
419	450	31	0.07
461	440	-21	-0.05
472	466	-6	-0.01
535	521	-14	-0.03
622	596	-26	-0.04
606	609	3	0.00
508	506	-2	0.00
461	453	-8	-0.02
390	407	17	0.04
432	450	18	0.04

# Compare three approaches: 

bind_rows(df_mlp %>% get_accuracy_measures(), 
          df_arima %>% get_accuracy_measures(), 
          df_ensemble %>% get_accuracy_measures()) %>% 
  mutate(Approach = c("MLP", "ARIMA", "Ensemble")) %>% 
  select(Approach, everything()) %>% 
  mutate_if(is.numeric, function(x) {round(x, 2)}) %>% 
  knitr::kable()

Approach	MAE	MSE	MAPE	N
MLP	24.17	736.67	5.02	12
ARIMA	18.58	576.75	4.19	12
Ensemble	13.00	256.83	2.80	12

References

1. Kourentzes, N., Barrow, D. K., & Crone, S. F. (2014). Neural network ensemble operators for time series forecasting. Expert Systems with Applications, 41(9), 4235-4244.
2. Hyndman, R. J., and Y. Khandakar (2008):“Automatic time series forecasting: The forecast package for R,”. Journal of Statistical Software, 26(3).
3. Barrow, D. K., Crone, S. F., & Kourentzes, N. (2010, July). An evaluation of neural network ensembles and model selection for time series prediction. In Neural Networks (IJCNN), The 2010 International Joint Conference on (pp. 1-8). IEEE.