44  Arima Forecasts

As usual, we start by uploading the libraries we’re going to use.

# ETL
library(tidyverse)
library(sparkline)
library(data.table)
library(DT)
library(scales)

# Forecasting
library(forecast)

# Charts
library(highcharter)
library(RColorBrewer)

# Tables
library(reactable)

45 Set Up Raw Data

Upload Raw Data

We’re going to upload a dataset from the file “actual_sales_data.csv”. This data frame contains historical sales quantities per product and per country.

We will then calculate some statistical forecasts, using the time serie approach ARIMA (Auto Regressive Integrated Moving Average).

Several steps are pretty similar to the ones we went through in the previous chapter, so we will go quickly through them.

# Upload dataset

# Define the URL of the raw CSV file
url <- "https://raw.githubusercontent.com/nguyennico/sales_forecasting_practice/main/actual_sales_data.csv"

# Read the CSV file from the URL
df1 <- read.csv(url)

# pivot
df1 <- df1 |> gather(key = "period", 
                     value = "actuals_qty", 
                     3:length(df1))

# remove the "X" in front of the period
df1$period <- gsub("X", "",df1$period)

# Format Date
df1$period <- as.Date(df1$period, format = '%m.%d.%Y')

# keep results
initial_data <- df1

head(initial_data)

  country product_description     period actuals_qty
1   Spain            ProductA 2016-01-01        1370
2  Mexico            ProductB 2016-01-01           0
3  Brazil            ProductB 2016-01-01          12
4  Brazil            ProductC 2016-01-01         369
5 Germany            ProductD 2016-01-01           0
6  France            ProductA 2016-01-01         263

Transform

# set a working df
df1 <- initial_data

# create a field DFU
df1$DFU <- paste(df1$country, df1$product_description, sep="_")

# keep only needed columns and get results
df1 <- df1 |> select(DFU, period, actuals_qty)


#--------------
# add some Calendar Features
#--------------

# add the calendar_year
df1$calendar_year <- year(df1$period)

# add the calendar_month
df1$calendar_month <- month(df1$period)

# transform the Month number into a month abbreviation
# note that month.abb comes with [ and not (
df1$calendar_month_abb <- month.abb[df1$calendar_month]


# create a Factor for the calendar_month_abb
# setting the beginning to January
df1$calendar_month_abb <- factor(df1$calendar_month_abb,
                              levels= c("Jan","Feb","Mar","Apr","May","Jun",
                                        "Jul","Aug","Sep", "Oct", "Nov", "Dec"))

# Get Results
actuals_data <- df1

head(actuals_data)

               DFU     period actuals_qty calendar_year calendar_month
1   Spain_ProductA 2016-01-01        1370          2016              1
2  Mexico_ProductB 2016-01-01           0          2016              1
3  Brazil_ProductB 2016-01-01          12          2016              1
4  Brazil_ProductC 2016-01-01         369          2016              1
5 Germany_ProductD 2016-01-01           0          2016              1
6  France_ProductA 2016-01-01         263          2016              1
  calendar_month_abb
1                Jan
2                Jan
3                Jan
4                Jan
5                Jan
6                Jan

Select Product

To illustrate our calculation we will look at one product : “Spain_ProductA” .

# select DFU to be forecasted
df1 <- actuals_data  |> filter(DFU == "Spain_ProductA")

# keep results
sales_data <- df1

head(sales_data)

             DFU     period actuals_qty calendar_year calendar_month
1 Spain_ProductA 2016-01-01        1370          2016              1
2 Spain_ProductA 2016-02-01        1528          2016              2
3 Spain_ProductA 2016-03-01        1101          2016              3
4 Spain_ProductA 2016-04-01         738          2016              4
5 Spain_ProductA 2016-05-01        1229          2016              5
6 Spain_ProductA 2016-06-01        1451          2016              6
  calendar_month_abb
1                Jan
2                Feb
3                Mar
4                Apr
5                May
6                Jun

46 Time Series Forecasting (Arima)

Let’s calculate the Statistical Forecasts, using the Time Serie approach Arima .

46.1 Create Historical Time Serie object

We consider the actuals over the Full Historical Window, and transform them into a time series object .

46.1.1 Get actuals

# set a working df
df1 <- sales_data

# select Period (historical window)
df1 <- df1 |> filter(period >= "2016-01-01" & period <= "2019-06-01")

# keep only needed columns
df1 <- df1 |> select(period, actuals_qty)

# keep a dataframe for charts later on
selected_actuals_data <- df1

head(df1)

      period actuals_qty
1 2016-01-01        1370
2 2016-02-01        1528
3 2016-03-01        1101
4 2016-04-01         738
5 2016-05-01        1229
6 2016-06-01        1451

46.1.2 Convert into a time series object

# set a working df
df1 <- selected_actuals_data

# rename
df1 <- df1 |> rename(Sales = actuals_qty)

# remove the period
df1 <- df1 |> select(-period)

# Convert it to a time series object
sales_data_ts <- ts(df1, 
                    start = c(2016,01),
                    end = c(2019,06),
                    frequency = 12)

# display
sales_data_ts

      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
2016 1370 1528 1101  738 1229 1451  879 1505 1375 1146 1325 1156
2017 1081 1258  894  700 1289 1207  926 1476 1254 1111 1175  881
2018  916 1225  682  739 1056 1117  727 1233  839  862  862  702
2019  740  934  603  468  903  901                              

46.2 Calculate Forecasts

Now, let’s calculate the statistical forecasts :

• using the function auto.arima().

• with 2 levels of confidence of 95%.

The function auto.arima() is from the R package forecast.

# create Forecasts
fore <- forecast(auto.arima(sales_data_ts), level = 95)

# see results
fore

         Point Forecast      Lo 95     Hi 95
Jul 2019       541.4113  328.32923  754.4933
Aug 2019      1081.5089  868.42685 1294.5910
Sep 2019       783.6805  570.59840  996.7625
Oct 2019       714.3688  501.28679  927.4509
Nov 2019       764.7717  551.68964  977.8537
Dec 2019       565.4435  352.36142  778.5255
Jan 2020       593.2256  380.94240  805.5088
Feb 2020       811.3800  599.09686 1023.6632
Mar 2020       413.1535  200.87027  625.4367
Apr 2020       299.8389   87.55569  512.1221
May 2020       730.0668  517.78358  942.3500
Jun 2020       752.8781  540.59487  965.1613
Jul 2020       377.6661  138.00660  617.3256
Aug 2020       917.7637  678.10423 1157.4232
Sep 2020       619.9353  380.27577  859.5948
Oct 2020       550.6237  310.96417  790.2832
Nov 2020       601.0265  361.36701  840.6860
Dec 2020       401.6983  162.03879  641.3578
Jan 2021       429.4804  190.53090  668.4299
Feb 2021       647.6349  408.68536  886.5844
Mar 2021       249.4083   10.45877  488.3578
Apr 2021       136.0937 -102.85581  375.0432
May 2021       566.3216  327.37208  805.2711
Jun 2021       589.1329  350.18337  828.0824

The results displays a calculated forecasts as well as a range of confidence.

The new object “fore” is a “forecast object”. We will need to transform it into a data frame in order to use it.

class(fore)

[1] "forecast"

46.3 Get the calculated Statistical Forecasts

# transform Forecasts File as a dataframe
df1 <- as.data.frame(fore)
    
# get the rownames
df1 <- setDT(df1, keep.rownames = TRUE)[]

# rename columns
colnames(df1)[1] <-"time"
colnames(df1)[2] <-"forecasts_qty"
    
# extract Year and Month
df1$year <- substr(df1$time, 5, 8)
df1$month_name <- substr(df1$time, 1, 3)
    
# create a period and a formatted period
df1$period <- paste(df1$month_name,"01", df1$year, sep = "/")
df1$period <- as.Date(df1$period, format = '%b/%d/%Y')


# keep only the needed columns
df1 <- df1 |> select(period, forecasts_qty)

# eventually format forecasts_qty into integer for a better display
df1$forecasts_qty <- as.integer(df1$forecasts_qty)
    
# replace negative values (sometimes Sales Forecasts become negative) by zero
df1$forecasts_qty <- if_else(df1$forecasts_qty < 0, 0, df1$forecasts_qty)

# Get Results
calculated_forecasts_data <- df1

head(df1)

       period forecasts_qty
       <Date>         <num>
1: 2019-07-01           541
2: 2019-08-01          1081
3: 2019-09-01           783
4: 2019-10-01           714
5: 2019-11-01           764
6: 2019-12-01           565

46.4 Prepare data for Chart

No we are going to combine the actuals and the calculated Forecasts.

#----------------------------
# prepare data for charting
#----------------------------

# selected_actuals_data
selected_actuals_data$forecasts_qty <- NA

# calculated_forecasts_data
calculated_forecasts_data$actuals_qty <- NA

# organize columns for stacking
calculated_forecasts_data <- calculated_forecasts_data |> select(period,
                                                                 actuals_qty,
                                                                 forecasts_qty
                                                                 ) 

# stack data
df1 <- rbind(selected_actuals_data,
             calculated_forecasts_data)

# display results
head(df1)

      period actuals_qty forecasts_qty
1 2016-01-01        1370            NA
2 2016-02-01        1528            NA
3 2016-03-01        1101            NA
4 2016-04-01         738            NA
5 2016-05-01        1229            NA
6 2016-06-01        1451            NA

46.5 Chart

We can display a line chart with those 2 components, using the library highcharter .

highchart() |>
  
  hc_add_series(name = "Actuals", 
                color = "steelblue", 
                data = df1$actuals_qty) |>
  
  hc_add_series(name = "Forecasts", 
                color = "limegreen", 
                data = df1$forecasts_qty) |>
  
  hc_title(text = "MLR [Multi Linear Regression: Trend + Season] Sales Forecasting") |>
  hc_subtitle(text = "in units") |>
  hc_xAxis(categories = df1$period) |> 
  hc_add_theme(hc_theme_google())

Note:

• We also can chart directly the forecast object “fore”, displaying the level of confidence.

• for this we will use the function hchart() from highcharter.

hchart(fore) |> 
  hc_title(text = "ARIMA method Sales Forecasting") |>
  hc_subtitle(text = "in units") |> 
  hc_add_theme(hc_theme_google())

47 Train | Test | Validate

47.1 Define Train dataset

We stop 6 months earlier : i.e. Dec 2018.

47.1.1 Get actuals for MAPE

We select the actuals up to June 2019.

# set a working df
df1 <- sales_data

# select Period (historical window)
df1 <- df1 |> filter(period >= "2016-01-01" & period <= "2019-06-01")

# keep only needed columns
df1 <- df1 |> select(period, actuals_qty)

# keep a dataframe for charts later on
original_actuals_data <- df1

head(df1)

      period actuals_qty
1 2016-01-01        1370
2 2016-02-01        1528
3 2016-03-01        1101
4 2016-04-01         738
5 2016-05-01        1229
6 2016-06-01        1451

47.1.2 Get actuals for Statistical calculation

We select the actuals up to December 2018.

# set a working df
df1 <- sales_data

# select Period (historical window)
df1 <- df1 |> filter(period >= "2016-01-01" & period <= "2018-12-01")

# keep only needed columns
df1 <- df1 |> select(period, actuals_qty)

# keep a dataframe for charts later on
selected_actuals_data <- df1

head(df1)

      period actuals_qty
1 2016-01-01        1370
2 2016-02-01        1528
3 2016-03-01        1101
4 2016-04-01         738
5 2016-05-01        1229
6 2016-06-01        1451

47.1.3 Convert into a time series object

# set a working df
df1 <- selected_actuals_data

# rename
df1 <- df1 |> rename(Sales = actuals_qty)

# remove the period
df1 <- df1 |> select(-period)

# Convert it to a time series object
sales_data_ts <- ts(df1, 
                    start = c(2016,01),
                    end = c(2018,12), # December 2016
                    frequency = 12)

# display results
sales_data_ts

      Jan  Feb  Mar  Apr  May  Jun  Jul  Aug  Sep  Oct  Nov  Dec
2016 1370 1528 1101  738 1229 1451  879 1505 1375 1146 1325 1156
2017 1081 1258  894  700 1289 1207  926 1476 1254 1111 1175  881
2018  916 1225  682  739 1056 1117  727 1233  839  862  862  702

47.2 Calculate Forecasts

# create Forecasts
fore <- forecast(auto.arima(sales_data_ts), level = 95)

# see results
fore

         Point Forecast     Lo 95     Hi 95
Jan 2019        755.875 508.76464 1002.9854
Feb 2019       1064.875 817.76464 1311.9854
Mar 2019        521.875 274.76464  768.9854
Apr 2019        578.875 331.76464  825.9854
May 2019        895.875 648.76464 1142.9854
Jun 2019        956.875 709.76464 1203.9854
Jul 2019        566.875 319.76464  813.9854
Aug 2019       1072.875 825.76464 1319.9854
Sep 2019        678.875 431.76464  925.9854
Oct 2019        701.875 454.76464  948.9854
Nov 2019        701.875 454.76464  948.9854
Dec 2019        541.875 294.76464  788.9854
Jan 2020        595.750 246.28317  945.2168
Feb 2020        904.750 555.28317 1254.2168
Mar 2020        361.750  12.28317  711.2168
Apr 2020        418.750  69.28317  768.2168
May 2020        735.750 386.28317 1085.2168
Jun 2020        796.750 447.28317 1146.2168
Jul 2020        406.750  57.28317  756.2168
Aug 2020        912.750 563.28317 1262.2168
Sep 2020        518.750 169.28317  868.2168
Oct 2020        541.750 192.28317  891.2168
Nov 2020        541.750 192.28317  891.2168
Dec 2020        381.750  32.28317  731.2168

47.3 Get the calculated Statistical Forecasts

# transform Forecasts File as a dataframe
df1 <- as.data.frame(fore)
    
# get the rownames
df1 <- setDT(df1, keep.rownames = TRUE)[]

# rename columns
colnames(df1)[1] <-"time"
colnames(df1)[2] <-"forecasts_qty"
    
# extract Year and Month
df1$year <- substr(df1$time, 5, 8)
df1$month_name <- substr(df1$time, 1, 3)
    
# create a period and a formatted period
df1$period <- paste(df1$month_name,"01", df1$year, sep = "/")
df1$period <- as.Date(df1$period, format = '%b/%d/%Y')


# keep only the needed columns
df1 <- df1 |> select(period, forecasts_qty)

# eventually format forecasts_qty into integer for a better display
df1$forecasts_qty <- as.integer(df1$forecasts_qty)
    
# replace negative values (sometimes Sales Forecasts become negative) by zero
df1$forecasts_qty <- if_else(df1$forecasts_qty < 0, 0, df1$forecasts_qty)

# Get Results
calculated_forecasts_data <- df1

head(df1)

       period forecasts_qty
       <Date>         <num>
1: 2019-01-01           755
2: 2019-02-01          1064
3: 2019-03-01           521
4: 2019-04-01           578
5: 2019-05-01           895
6: 2019-06-01           956

47.4 Calculation of the MAPE

# merging original_actuals_data and calculated_forecasts_data
df1 <-left_join(original_actuals_data, calculated_forecasts_data)

# calculate monthly difference
df1$delta <- df1$forecasts_qty - df1$actuals_qty
df1$delta_pc <- df1$delta / df1$actuals_qty


#----------------------------
# calculate the Mean absolute percentage error (MAPE)
#----------------------------

# get the absolute percentage of error
df1$abs_delta_pc <- abs(df1$delta_pc)

# define timeframe for the measurement of the MAPE

# start date
Start.Date<- '2019-01-01'

# end date
End.Date<- '2019-06-01'

# select Period (historical window)
df2 <- df1 |> filter(period >= Start.Date & df1$period <= End.Date)

# display
df2

      period actuals_qty forecasts_qty delta     delta_pc abs_delta_pc
1 2019-01-01         740           755    15  0.020270270  0.020270270
2 2019-02-01         934          1064   130  0.139186296  0.139186296
3 2019-03-01         603           521   -82 -0.135986733  0.135986733
4 2019-04-01         468           578   110  0.235042735  0.235042735
5 2019-05-01         903           895    -8 -0.008859358  0.008859358
6 2019-06-01         901           956    55  0.061043285  0.061043285

Let’s calculate the MAPE over this horizon :

# calculate the mean
MAPE <- mean(df2$abs_delta_pc)
MAPE <- percent(MAPE)

MAPE

[1] "10%"