# ETL
library(tidyverse)
library(sparkline)
library(data.table)
library(DT)
library(scales)
# Forecasting
library(forecast)
# Charts
library(highcharter)
library(RColorBrewer)
# Tables
library(reactable)45 Arima Forecasts
As usual, we start by uploading the libraries we’re going to use.
46 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 263Transform
# 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 JanSelect 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 Jun47 Time Series Forecasting (Arima)
Let’s calculate the Statistical Forecasts, using the Time Serie approach Arima .
47.1 Create Historical Time Serie object
We consider the actuals over the Full Historical Window, and transform them into a time series object .
47.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 145147.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 47.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.0824The 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"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-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 56547.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 NA47.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()fromhighcharter.
hchart(fore) |>
hc_title(text = "ARIMA method Sales Forecasting") |>
hc_subtitle(text = "in units") |>
hc_add_theme(hc_theme_google())48 Train | Test | Validate
48.1 Define Train dataset
We stop 6 months earlier : i.e. Dec 2018.
48.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 145148.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 145148.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 70248.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.216848.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 95648.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.061043285Let’s calculate the MAPE over this horizon :
# calculate the mean
MAPE <- mean(df2$abs_delta_pc)
MAPE <- percent(MAPE)
MAPE[1] "10%"