Forecasting Daily Cash Demand: A Short Comparision between Automatic Machine Learning and ARIMA Aprroach (Case 1)

R for Pleasure

Nguyen Chi Dung

#==============================
#  Exploratory Data Analysis
#==============================

# Inspect data path: 

rm(list = ls())
my_path <- dir("/home/chidung/Desktop/atm_cash_data", full.names = TRUE)

# Import data: 

library(tidyverse)
library(magrittr)

my_case1 <- read_csv(my_path[1])
my_case2 <- read_csv(my_path[2])

# A function for black theme: 

my_theme <- function(...) {
  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"))
  
}

# Daily Withdrawals of cash for case 1 and case 2: 
my_case1 %>% 
  ggplot(aes(DAYID, daily_cash)) + 
  geom_line(color = "cyan") + 
  my_theme() + 
  labs(x = NULL, y = NULL, title = "Daily Cash Withdrawals (Terminal ID: 00002009)")

my_case2 %>% 
  ggplot(aes(DAYID, daily_cash)) + 
  geom_line(color = "cyan") + 
  my_theme() + 
  labs(x = NULL, y = NULL, title = "Daily Cash Withdrawals (Terminal ID: 00002019)")

# A function for data pre-processing: 

pre_proccess <- function(your_df) {
  your_df %>% 
    timetk::tk_augment_timeseries_signature() %>% 
    mutate_if(is.ordered, as.character) %>% 
    mutate_if(is.character, as.factor)%>% 
    na.omit() %>% 
    select(-hour, -minute, -second, -am.pm, -index.num, -year.iso, -year, -month, -day) %>% 
    return()
}


# Used the function: 
pre_proccess(my_case1) ->> cash_clean

#=======================================
#  Automatic Machine Learning Approach
#=======================================

# Wrire a function to convert to h2o frame: 

convert_to_h2o <- function(your_df) {
  your_df %>% 
    h2o::as.h2o() %>% 
    return()
}

# Wrire a function for conducting One-step Ahead Prediction: 

library(h2o)
h2o.init(nthreads = 32, max_mem_size = "16G")

##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         1 hours 1 minutes 
##     H2O cluster timezone:       Asia/Ho_Chi_Minh 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.20.0.2 
##     H2O cluster version age:    3 months and 19 days !!! 
##     H2O cluster name:           H2O_started_from_R_chidung_rjy105 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   15.96 GB 
##     H2O cluster total cores:    32 
##     H2O cluster allowed cores:  32 
##     H2O cluster healthy:        TRUE 
##     H2O Connection ip:          localhost 
##     H2O Connection port:        54321 
##     H2O Connection proxy:       NA 
##     H2O Internal Security:      FALSE 
##     H2O API Extensions:         XGBoost, Algos, AutoML, Core V3, Core V4 
##     R Version:                  R version 3.5.1 (2018-07-02)

h2o.no_progress()
h2o.removeAll()

## [1] 0

my_predict <- function(k) {
  
  range_validation <- 60
  n_ahead <- 1
  
  train_tbl <- cash_clean %>% 
    slice((1 + k):(400 + k))
  
  valid_tbl <- cash_clean %>% 
    slice((400 + 1 + k):(400 + k + range_validation)) 
  
  test_tbl  <- cash_clean %>% 
    slice((400 + k + range_validation + 1):(400 + k + range_validation + n_ahead))
  
  
  train_h2o <- train_tbl %>% select(-DAYID) %>% convert_to_h2o()
  valid_h2o <- valid_tbl %>% select(-DAYID) %>% convert_to_h2o()
  test_h2o <- test_tbl %>% select(-DAYID) %>% convert_to_h2o()
  
  # Set target and features: 
  
  y <- "daily_cash"
  x <- setdiff(names(train_h2o), y)
  
  # Automatic ML: 
  automl_models_h2o <- h2o.automl(x = x, y = y, 
                                  training_frame = train_h2o,
                                  validation_frame = valid_h2o, 
                                  leaderboard_frame = test_h2o, 
                                  max_runtime_secs = 60, 
                                  max_models = 30, 
                                  stopping_metric = "MAE", 
                                  seed = 29)
  
  
  # Use for forecasting: 
  pred_h2o <- h2o.predict(automl_models_h2o, newdata = test_h2o) %>% 
    as.data.frame() %>% 
    pull(predict)
  
  actual_predicted_df_test <- test_tbl %>% 
    mutate(predicted = pred_h2o, error = daily_cash - predicted, pct = error / daily_cash) %>% 
    select(DAYID, daily_cash, predicted, everything())
  
  return(actual_predicted_df_test)
  
}

# Function for conducting One-step Ahead Prediction for (n + 1) consecutive days: 

my_predict_n_days <- function(n_days) {
  results <- lapply(0:n_days, my_predict)
  results_df <- do.call("bind_rows", results)
  return(results_df)
  
}

# Use this function: 

system.time(my_predict_n_days(30) ->> results_df)

##     user   system  elapsed 
##   51.868    1.977 2019.704

# Function visualizes predictions vs actuals: 

vis_results <- function(r_df) {
  r_df %>% 
    select(DAYID, Actual = daily_cash, Predicted = predicted) %>% 
    gather(a, b, -DAYID) %>% 
    ggplot(aes(DAYID, b, color = a)) + 
    geom_line() + 
    geom_point() + 
    scale_color_manual(values = c("purple", "orange"), name = "") + 
    my_theme()
  
}

# Visualize predicted - actual results
results_df %>% 
  vis_results() + 
  labs(x = NULL, y = NULL, subtitle = "Approach Used: Automatic Machine Learning", 
       title = "One-step Ahead Prediction for Cash Withdrawals (31 Consecutive Days)")

# Function presents predicted - actual results: 

print_results <- function(df) {
  df %>% 
    select(DAYID, Actual = daily_cash, Predicted = predicted) %>% 
    mutate_if(is.numeric, function(x) {round(x, 2)}) %>% 
    knitr::kable()
}

results_df %>% 
  print_results()

# Compare total actual vs predicted demands (31 days): 

results_df %>% 
  select(DAYID, daily_cash, predicted) %>% 
  gather(a, value, -DAYID) %>% 
  group_by(a) %>% 
  summarise(cash = sum(value))

## # A tibble: 2 x 2
##   a           cash
##   <chr>      <dbl>
## 1 daily_cash 2314.
## 2 predicted  2166.

#========================
#     ARIMA Approach
#========================

library(fpp2)

predict_from_arima <- function(k) {
  
  range_validation <- 60
  n_ahead <- 1
  
  train_tbl <- cash_clean %>% slice((1 + k):(400 + k))
  valid_tbl <- cash_clean %>% slice((400 + 1 + k):(400 + k + range_validation)) 
  test_tbl  <- cash_clean %>% slice((400 + k + range_validation + 1):(400 + k + range_validation + n_ahead))
  
  train_arima <- bind_rows(train_tbl, valid_tbl) %>% select(1:2)
  test_arima <- test_tbl %>% select(1:2)
  
  # Automatic ARIMA model as proposed by Hyndman and Khandakar (2008): 
  my_arima <- auto.arima(train_arima[, 2] %>% ts(start = 1))
  
  # Use the model for forecasting: 
  predicted_arima <- forecast(my_arima, h = 1)$mean %>% as.vector()
  
  actual_predicted_df_test <- test_arima %>% 
    mutate(predicted = predicted_arima) 
  
  return(actual_predicted_df_test)
  
}

# Use this function: 

system.time(lapply(0:30, predict_from_arima) ->> arima_results)

##    user  system elapsed 
##  13.460  21.515   5.335

arima_results <- do.call("bind_rows", arima_results)

# Compare cash demands predicted by ARIMA approach and actuals: 

arima_results %>% 
  vis_results() + 
  labs(x = NULL, y = NULL, subtitle = "Approach Used: ARIMA", 
       title = "One-step Ahead Prediction for Cash Withdrawals (31 Consecutive Days)")

#  Presents predicted - actual results: 
arima_results %>% 
  print_results()

#====================================================================================
#  Comparision between the two approaches besed on some Forecast Accuracy Measures
#====================================================================================

# Function calculate some accuracy measures: 

get_accuracy_measures <- function(your_result_df) {
  
  act <- your_result_df %>% pull(daily_cash)
  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, trim = 5)
  
  # Return results: 
  return(data.frame(MAE = mae, MSE = mse, MAPE = mape))
}


# Use above function for comparing: 

bind_rows(results_df %>% get_accuracy_measures(), arima_results %>% get_accuracy_measures) %>% 
  mutate(Approach = c("Automatic ML", "ARIMA")) %>% 
  select(Approach, everything()) %>% 
  mutate_if(is.numeric, function(x) {round(x, 2)}) %>% 
  knitr::kable()