# ETL
library(tidyverse)
library(sparkline)
library(data.table)
library(DT)
library(scales)
# Forecasting
library(forecast)
# Charts
library(highcharter)
library(RColorBrewer)
# Tables
library(reactable)41 MLR Forecasts
As usual, we start by uploading the libraries we’re going to use.
42 Set Up Raw Data
Upload Raw Data
We’re going to upload a data frame from the file “actual_sales_data.csv”. This dataset contains historical sales quantities per product and per country.
We will then calculate some statistical forecasts, using the time serie approach MLR (Multi Linear Regression).
# 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
We create a new variable DFU (stands for Demand Forecast Unit) and add some dates features (calendar year, calendar month, etc).
# 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 Jan43 Visualize Historical Data
Let’s visualize the data, taking the example of the DFU “Spain_ProductA”.
We’re going to create 3 classic views :
(full) historical sales.
MTM Line Chart (seasonality).
Bar Chart (YTD trend).
Select Product
# 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 Jun43.1 Historical Sales
Always interesting to see the full historical sales.
# set working dataframe
df1 <- sales_data
#--------------------------------------
# chart
highchart() |>
hc_add_series(name = "Actuals",
color = "steelblue",
data = df1$actuals_qty) |>
hc_title(text = "Actual Sales") |>
hc_subtitle(text = "in units") |>
hc_xAxis(categories = df1$period) |>
#hc_yAxis(title = list(text = "Sales (units)")) |>
hc_add_theme(hc_theme_google())43.2 MTM Line Chart (seasonality)
Let’s check whether we have a recurrent seasonality.
For this :
we can display the sales through a MTM chart, and identify some month with a regular pattern (i.e. some increases and de decreases) which occur each year.
we need to have at least 2 full years of historical data to be able to identify a regular seasonality, if any.
# set working dataframe
df1 <- sales_data
# keep only the needed columns
df1 <- df1 |> select(calendar_month_abb,
calendar_year,
actuals_qty)
# spread
df1 <- df1 |> spread(calendar_year, actuals_qty)
#------------------------
# chart
highchart() |>
hc_add_series(name = "2016",
color = "lightgray",
data = df1$`2016`) |>
hc_add_series(name = "2017",
color = "lightblue",
data = df1$`2017`) |>
hc_add_series(name = "2018",
color = "lightgreen",
data = df1$`2018`) |>
hc_add_series(name = "2019",
color = "gold",
data = df1$`2019`) |>
hc_title(text = "MTM Actual Sales") |>
hc_subtitle(text = "in units") |>
hc_xAxis(categories = df1$calendar_month_abb) |>
#hc_yAxis(title = list(text = "Sales (units)")) %>%
hc_add_theme(hc_theme_google())We can see that, though the level of sales of each year varies, there is a regular pattern over the different years :
sales peaks in February, June and August, followed by a decrease.
low sales in April and July.
43.3 Bar Chart (YTD trend)
Now, let’s see whether we have a trend. We can look at the sales over the first 6months, during different years.
# set working dataframe
df1 <- sales_data
# select only the 1st 6 months
df1 <- df1 |> filter(calendar_month <= 6)
# aggregate
df1 <- df1 |> group_by(calendar_year) |>
summarise(actuals_qty = sum(actuals_qty)
)
#----------------
# chart
highchart() |>
hc_title(text = "YTD Sales Volumes") |>
hc_subtitle(text = "in units") |>
hc_add_theme(hc_theme_google()) |>
hc_xAxis(categories = df1$calendar_year) |>
hc_add_series(name = "Sales",
color = "mediumseagreen",
dataLabels = list(align = "center", enabled = TRUE),
data = df1$actuals_qty) |>
hc_chart(type = "column") We notice that year over year, over a similar horizon of time, the sales are decreasing. There is then a trend, that we will consider when we calculate the Statistical Forecasts.
44 Time Series Forecasting (MLR)
Let’s calculate the Statistical Forecasts, using a Time Serie approach called MLR (Multi Linear Regression).
44.1 Create Historical Time Serie object
We consider the actuals over the Full Historical Window, and transform them into a time series object .
44.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 145144.1.2 Convert into a time series object
We use the function ts() from the R package forecast .
# 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 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
2019 740 934 603 468 903 901 44.2 Calculate Forecasts
Now, let’s calculate the statistical forecasts :
using the function
tslm()of the R packageforecast.based on the trend and seasonality.
over the next 24 months horizon.
with a level of confidence of +/- 10%.
The function tslm() is used to fit linear models to time series including trend and seasonality components.
# create model MLR
train.lm.trend.season <- tslm(sales_data_ts ~ trend + season)
# create Forecasts MLR
train.lm.trend.season.pred <- forecast(train.lm.trend.season,
h = 24,
level = 10)
# see results
train.lm.trend.season.pred Point Forecast Lo 10 Hi 10
Jul 2019 519.7143 506.36565 533.0629
Aug 2019 1080.3810 1067.03231 1093.7296
Sep 2019 831.7143 818.36565 845.0629
Oct 2019 715.3810 702.03231 728.7296
Nov 2019 796.3810 783.03231 809.7296
Dec 2019 588.7143 575.36565 602.0629
Jan 2020 621.3929 608.18400 634.6017
Feb 2020 830.8929 817.68400 844.1017
Mar 2020 414.6429 401.43400 427.8517
Apr 2020 255.8929 242.68400 269.1017
May 2020 713.8929 700.68400 727.1017
Jun 2020 763.6429 750.43400 776.8517
Jul 2020 357.5714 343.67772 371.4651
Aug 2020 918.2381 904.34439 932.1318
Sep 2020 669.5714 655.67772 683.4651
Oct 2020 553.2381 539.34439 567.1318
Nov 2020 634.2381 620.34439 648.1318
Dec 2020 426.5714 412.67772 440.4651
Jan 2021 459.2500 445.38304 473.1170
Feb 2021 668.7500 654.88304 682.6170
Mar 2021 252.5000 238.63304 266.3670
Apr 2021 93.7500 79.88304 107.6170
May 2021 551.7500 537.88304 565.6170
Jun 2021 601.5000 587.63304 615.3670The results displays a calculated forecasts as well as a range of confidence.
The new object “train.lm.trend.season.pred” is a “forecast object”. We will need to transform it into a data frame in order to use it.
class(train.lm.trend.season.pred)[1] "forecast"44.3 Get the calculated Statistical Forecasts
Let’s transform the object train.lm.trend.season.pred into a data frame and :
get a proper period variable.
adjust the calculated forecasts, replacing the negative values by zero, if any.
# transform Forecasts File as a dataframe
df1 <- as.data.frame(train.lm.trend.season.pred)
# 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 519
2: 2019-08-01 1080
3: 2019-09-01 831
4: 2019-10-01 715
5: 2019-11-01 796
6: 2019-12-01 58844.4 Prepare data for Chart
When we look at the sales data, it’s good to have the actuals and the calculated forecasts altogether.
So now we are going to combine those 2 variables into a common data frame.
#----------------------------
# 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
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 NA44.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())We will need to go further into our analysis : how do we know how good are our calculated forecasts?
To answer this question we will need to perform a Train | Test | Validate calculation and check (validate) how accurate is our model, comparing the calculated forecasts versus some actual values.
45 Train | Test | Validate
45.1 Define Training dataset
The methodology is pretty simple :
we stop the historical values 6 months earlier : i.e. Dec 2018.
- this will become our training dataset.
we calculate the forecasts over the next 6 months : from Jan 2019 until Jun 2019.
we compare those 6 months of forecasts versus the historical values.
- we test how relevant are the forecasts using a MAPE (Moving Average Percentage of Error) indicator and validate if the forecasts are relevant.
45.1.1 Get actuals for MAPE
We select the actuals up to June 2019. We will later on use those data to test our model.
# 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 145145.1.2 Get actuals for Statistical calculation
We select the actuals up to December 2018. We will use those data to train our model.
# 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 145145.1.3 Convert into a time series object
We will use this time series object called sales_data_ts to calculate the statistical forecasts.
# 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
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 70245.2 Calculate Forecasts
We apply the function tslm() and calculate the statistical forecasts for the next 24 months, displaying a confidence level of +/- 10% .
We obtain an object called train.lm.trend.season.pred .
# create model MLR
train.lm.trend.season <- tslm(sales_data_ts ~ trend + season)
# create Forecasts MLR
train.lm.trend.season.pred <- forecast(train.lm.trend.season,
h = 24,
level = 10)
# see results
train.lm.trend.season.pred Point Forecast Lo 10 Hi 10
Jan 2019 802.0833 787.1510 817.0157
Feb 2019 1016.7500 1001.8177 1031.6823
Mar 2019 572.0833 557.1510 587.0157
Apr 2019 405.4167 390.4843 420.3490
May 2019 871.0833 856.1510 886.0157
Jun 2019 938.0833 923.1510 953.0157
Jul 2019 523.7500 508.8177 538.6823
Aug 2019 1084.4167 1069.4843 1099.3490
Sep 2019 835.7500 820.8177 850.6823
Oct 2019 719.4167 704.4843 734.3490
Nov 2019 800.4167 785.4843 815.3490
Dec 2019 592.7500 577.8177 607.6823
Jan 2020 641.9583 626.0227 657.8939
Feb 2020 856.6250 840.6894 872.5606
Mar 2020 411.9583 396.0227 427.8939
Apr 2020 245.2917 229.3561 261.2273
May 2020 710.9583 695.0227 726.8939
Jun 2020 777.9583 762.0227 793.8939
Jul 2020 363.6250 347.6894 379.5606
Aug 2020 924.2917 908.3561 940.2273
Sep 2020 675.6250 659.6894 691.5606
Oct 2020 559.2917 543.3561 575.2273
Nov 2020 640.2917 624.3561 656.2273
Dec 2020 432.6250 416.6894 448.560645.3 Get the calculated Statistical Forecasts
We transform the object train.lm.trend.season.pred into a data frame and :
get a proper period variable.
adjust the calculated forecasts, replacing the negative values by zero, if any.
# transform Forecasts File as a dataframe
df1 <- as.data.frame(train.lm.trend.season.pred)
# 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 802
2: 2019-02-01 1016
3: 2019-03-01 572
4: 2019-04-01 405
5: 2019-05-01 871
6: 2019-06-01 93845.4 Calculation of the MAPE
And finally we can calculate our MAPE.
We find a MAPE of 7%, which is very low. We can conclude that our model is statistically very relevant, using a Multi Linear Regression model.
When we look at the values of the variable delta_pc, we indeed can notice that for most of the months the calculated forecasts are pretty close to the actuals.
# 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 802 62 0.08378378 0.08378378
2 2019-02-01 934 1016 82 0.08779443 0.08779443
3 2019-03-01 603 572 -31 -0.05140962 0.05140962
4 2019-04-01 468 405 -63 -0.13461538 0.13461538
5 2019-05-01 903 871 -32 -0.03543743 0.03543743
6 2019-06-01 901 938 37 0.04106548 0.04106548Let’s calculate the MAPE over this horizon :
# calculate the mean
MAPE <- mean(df2$abs_delta_pc)
MAPE <- percent(MAPE)
MAPE[1] "7%"