# ETL
library(tidyverse)
library(sparkline)
library(data.table)
library(DT)
library(scales)
# Forecasting
library(forecast)
# Charts
library(highcharter)
library(RColorBrewer)
# Tables
library(reactable)44 Holts Winter Forecasts
As usual, we start by uploading the libraries we’re going to use.
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 Holts Winter (also called exponential smoothing).
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 Jun46 Time Series Forecasting (Holts Winter)
Let’s calculate the Statistical Forecasts, using the Time Serie approach Holts Winter.
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 145146.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 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 46.2 Calculate Forecasts
Now, let’s calculate the statistical forecasts :
using the function
ets().based on the trend and seasonality.
over the next 24 months horizon.
with 2 levels of confidence of +/- 10% and +/- 5%.
The function ets() is from the R package forecast. ets() is used to fit exponential time series models.
# create model MLR
fit <- ets(sales_data_ts)
# create Forecasts
fore <- forecast(fit, h = 24,
level = c(10, 5))
# see results
fore Point Forecast Lo 5 Hi 5 Lo 10 Hi 10
Jul 2019 598.8637 595.2686 602.4587 591.6594 606.0679
Aug 2019 967.7825 961.9601 973.6049 956.1146 979.4503
Sep 2019 793.6741 788.8863 798.4620 784.0795 803.2688
Oct 2019 706.8857 702.6073 711.1640 698.3121 715.4592
Nov 2019 753.5705 748.9915 758.1494 744.3945 762.7465
Dec 2019 607.2250 603.5179 610.9322 599.7962 614.6539
Jan 2020 624.0588 620.2276 627.8899 616.3814 631.7362
Feb 2020 746.7285 742.1143 751.3426 737.4820 755.9750
Mar 2020 489.6070 486.5587 492.6554 483.4983 495.7158
Apr 2020 399.9195 397.4077 402.4313 394.8859 404.9531
May 2020 665.0583 660.8390 669.2776 656.6031 673.5136
Jun 2020 677.0438 672.6989 681.3887 668.3369 685.7507
Jul 2020 455.1809 452.2214 458.1403 449.2503 461.1114
Aug 2020 730.8494 726.0272 735.6715 721.1860 740.5127
Sep 2020 595.3196 591.3262 599.3129 587.3170 603.3221
Oct 2020 526.4636 522.8662 530.0610 519.2546 533.6726
Nov 2020 557.0530 553.1674 560.9385 549.2666 564.8394
Dec 2020 445.3540 442.1760 448.5321 438.9854 451.7227
Jan 2021 453.9207 450.5991 457.2424 447.2643 460.5772
Feb 2021 538.4141 534.3639 542.4643 530.2976 546.5305
Mar 2021 349.7708 347.0590 352.4826 344.3365 355.2051
Apr 2021 282.9140 280.6473 285.1808 278.3716 287.4565
May 2021 465.6179 461.7520 469.4839 457.8708 473.3651
Jun 2021 468.8052 464.7600 472.8505 460.6988 476.9117The results display 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 598
2: 2019-08-01 967
3: 2019-09-01 793
4: 2019-10-01 706
5: 2019-11-01 753
6: 2019-12-01 60746.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
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 NA46.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 = "ETS method Sales Forecasting") |>
hc_subtitle(text = "in units") |>
hc_add_theme(hc_theme_google())47 Train | Test | Validate
47.1 Define Training 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 145147.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 145147.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
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 70247.2 Calculate Forecasts
# create model MLR
fit <- ets(sales_data_ts)
# create Forecasts
fore <- forecast(fit, h = 24,
level = c(10, 5))
# see results
fore Point Forecast Lo 5 Hi 5 Lo 10 Hi 10
Jan 2019 829.5889 824.3744 834.8034 819.1393 840.0385
Feb 2019 977.7680 971.6220 983.9140 965.4517 990.0843
Mar 2019 652.3825 648.2816 656.4834 644.1646 660.6004
Apr 2019 548.7463 545.2966 552.1959 541.8334 555.6592
May 2019 863.2602 857.8328 868.6876 852.3840 874.1364
Jun 2019 900.1029 894.4430 905.7628 888.7607 911.4451
Jul 2019 605.2166 601.4101 609.0232 597.5886 612.8447
Aug 2019 1006.9898 1000.6544 1013.3252 994.2940 1019.6857
Sep 2019 818.3623 813.2115 823.5130 808.0404 828.6841
Oct 2019 724.2141 719.6536 728.7746 715.0751 733.3531
Nov 2019 783.6456 778.7076 788.5835 773.7502 793.5410
Dec 2019 633.5601 629.5647 637.5555 625.5535 641.5667
Jan 2020 674.8268 670.5670 679.0866 666.2903 683.3633
Feb 2020 792.4822 787.4738 797.4906 782.4456 802.5189
Mar 2020 526.7733 523.4394 530.1071 520.0924 533.4542
Apr 2020 441.3682 438.5702 444.1662 435.7611 446.9753
May 2020 691.5382 687.1456 695.9308 682.7357 700.3407
Jun 2020 718.0339 713.4625 722.6053 708.8731 727.1947
Jul 2020 480.6971 477.6285 483.7656 474.5478 486.8463
Aug 2020 796.1937 791.0954 801.2921 785.9769 806.4105
Sep 2020 644.0106 639.8720 648.1492 635.7170 652.3042
Oct 2020 567.1318 563.4723 570.7914 559.7983 574.4654
Nov 2020 610.5438 606.5856 614.5020 602.6117 618.4758
Dec 2020 490.9867 487.7865 494.1869 484.5737 497.399747.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 829
2: 2019-02-01 977
3: 2019-03-01 652
4: 2019-04-01 548
5: 2019-05-01 863
6: 2019-06-01 90047.4 Calculation of the MAPE
Finally we can calculate our MAPE.
We find a MAPE of 8%, which is very low. We can conclude that our model is statistically very relevant, using the Holts Winter model.
When we look at the values of the variable delta_pc, we 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 829 89 0.120270270 0.120270270
2 2019-02-01 934 977 43 0.046038544 0.046038544
3 2019-03-01 603 652 49 0.081260365 0.081260365
4 2019-04-01 468 548 80 0.170940171 0.170940171
5 2019-05-01 903 863 -40 -0.044296788 0.044296788
6 2019-06-01 901 900 -1 -0.001109878 0.001109878Let’s calculate the MAPE over this horizon :
# calculate the mean
MAPE <- mean(df2$abs_delta_pc)
MAPE <- percent(MAPE)
MAPE[1] "8%"