43  Holts Winter 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)

44 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         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

45 Time Series Forecasting (Holts Winter)

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

45.1 Create Historical Time Serie object

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

45.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

45.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                              

45.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.9117

The 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"

45.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           607

45.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            NA

45.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 = "ETS method Sales Forecasting") |>
  hc_subtitle(text = "in units") |> 
  hc_add_theme(hc_theme_google())

46 Train | Test | Validate

46.1 Define Training dataset

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

46.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

46.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

46.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  702

46.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.3997

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-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           900

46.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.001109878

Let’s calculate the MAPE over this horizon :

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

MAPE

[1] "8%"