Time Series Case Study

Assignment: Demand Forecasting and Exploratory Data Analysis

Introduction

Welcome to the demand forecasting and exploratory data analysis assignment! In this exercise, you will have the opportunity to demonstrate your skills in analysing and forecasting demand for a retail company’s products. The data set provided contains information about sales transactions, including the country of sale, department (brand), stock keeping unit (SKU), date of sale, and the number of items sold.

Objective

Your task is to perform exploratory data analysis (EDA) and build a demand forecasting model. The ultimate goal is to help the company make informed decisions regarding inventory management, stock replenishment, and sales strategy. You are free to choose your approach and models for demand forecasting, but you must explain your reasons for selecting them.

Dataset Description:

The dataset you will be working with contains the following columns: ● COUNTRY: The country in which the items were sold. ● DEPARTMENT: The brand or department for which the items are being sold. ● SKU: The unique ID of a stock-keeping unit. ● DATE: The date of the sales transaction. ● NUM_ITEMS_SOLD: The number of items sold on a daily level.

Exploratory Data Analysis

#Import needed Libraries
library(tidyverse)
library(data.table)
library(lubridate)
library(forecast)
library(tseries)
library(zoo)
library(imputeTS)
library(modeltime)
library(tidymodels)
library(dunn.test)
library(ggstatsplot)
library(dlookr)
library(effectsize)
library(tidyquant)
library(Metrics)
library(timetk)

#Read in data
raw_data = fread("data.csv")

glimpse(raw_data)

Rows: 1,915
Columns: 5
$ COUNTRY        <chr> "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "C1", "…
$ DEPARTMENT     <chr> "Dep1", "Dep1", "Dep1", "Dep1", "Dep1", "Dep1", "Dep1",…
$ SKU            <chr> "SKU_1", "SKU_1", "SKU_1", "SKU_1", "SKU_1", "SKU_1", "…
$ DATE           <IDate> 2021-11-29, 2021-11-30, 2021-12-01, 2021-12-02, 2021-…
$ NUM_ITEMS_SOLD <int> 1, 1, 3, 2, 1, 3, 1, 8, 1, 3, 0, 0, 0, 0, 1, 3, 3, 8, 9…

Using glimpse() function, we can see an overview of our data. We have 1915 rows and 5 columns. We have 3 character type columns COUNTRY, DEPARTMENT and SKU. DATE column is a date type object and NUM_ITEMS_SOLD is an integer type of object.

Let’s perform basic exploration on each of our columns.

raw_data %>% 
  ggplot(aes(x = COUNTRY)) +
  geom_bar()

C1 seems to have the largest count on our data set. Let’s continue exploring our data.

raw_data %>% 
  ggplot(aes(x = DEPARTMENT)) +
  geom_bar()

raw_data %>% 
  ggplot(aes(x = COUNTRY, fill = DEPARTMENT)) +
  geom_bar(position = "dodge")

Interesting. Looks like Dep2 is only present at Country 3 or C3.

raw_data %>% 
  group_by(SKU) %>% 
  reframe(
    SKU_count = n()
  ) %>% 
  arrange(desc(SKU_count)) %>% 
  head(n = 10)

# A tibble: 10 × 2
   SKU     SKU_count
   <chr>       <int>
 1 SKU_15        126
 2 SKU_66         70
 3 SKU_69         69
 4 SKU_6          65
 5 SKU_94         57
 6 SKU_46         50
 7 SKU_55         49
 8 SKU_89         48
 9 SKU_110        46
10 SKU_79         45

SKU column contains a lot of unique items, so it won’t make sense for us to visualize it as bar graph. It will be too congested for us to interpret it.

raw_data %>% 
  group_by(DATE) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  ggplot(aes(x = DATE, y = NUM_ITEMS_SOLD)) +
  geom_col()

The above plot shows the distribution of Number of Items sold over time. We can see that the sales started of below 500 but notice how we experienced a spike of sales on dates at around 2022-03-2022-04. After that spike, we can see how our sales started to decrease again and interestingly, we can see that there is at least 1 quarter of high sales before it started to drop again. We can further investigate this distribution buy breaking it down using our categorical variable.

raw_data %>% 
  group_by(DATE, COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  filter(COUNTRY == "C1") %>% 
  ggplot(aes(x = DATE, y = NUM_ITEMS_SOLD)) +
  geom_col()

Sales for C1 seems to be more consistent, even though we can still see time period where there are no sales or have low count of sales.

raw_data %>% 
  group_by(DATE, COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  filter(COUNTRY == "C2") %>% 
  ggplot(aes(x = DATE, y = NUM_ITEMS_SOLD)) +
  geom_col()

Sales for C2 are more sparse. This is expected because of the low count we saw from our bar graph comparing the count of the 3 country.

raw_data %>% 
  group_by(DATE, COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  filter(COUNTRY == "C3") %>% 
  ggplot(aes(x = DATE, y = NUM_ITEMS_SOLD)) +
  geom_col()

Even if Country 3 has more count than Country 2 it is still more sparse in terms of sales over time period. We can see that there is a bit of sales concentration on the end of 2022. It is also worth noticing that the spike we saw from the distribution graph without the country break down came from country 3.

Now that we have an idea of how the sales are distributed over time, let’s try to check what it will look like on a monthly view.

We need to create a month year column first out of the DATE column we already have.

raw_data = raw_data %>% 
  mutate(
    MONTH_YEAR = format(DATE,"%Y-%m"),
    MONTH_YEAR = paste(MONTH_YEAR,"-","01",sep = ""),
    MONTH_YEAR = as.Date(MONTH_YEAR)
  )

raw_data %>% 
  group_by(MONTH_YEAR) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  ggplot(aes(x = MONTH_YEAR, y = NUM_ITEMS_SOLD)) +
  geom_col()+
  ggtitle(
    "Monthly Sales Distribution"
  )

raw_data %>% 
  group_by(MONTH_YEAR,COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  filter(COUNTRY == "C1") %>% 
  ggplot(aes(x = MONTH_YEAR, y = NUM_ITEMS_SOLD)) +
  geom_col() +
  ggtitle(
    "Country 1 Monthly Sales Distribution"
  )

raw_data %>% 
  group_by(MONTH_YEAR,COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  filter(COUNTRY == "C2") %>% 
  ggplot(aes(x = MONTH_YEAR, y = NUM_ITEMS_SOLD)) +
  geom_col() +
  ggtitle(
    "Country 2 Monthly Sales Distribution"
  )

raw_data %>% 
  group_by(MONTH_YEAR,COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  ) %>% 
  filter(COUNTRY == "C3") %>% 
  ggplot(aes(x = MONTH_YEAR, y = NUM_ITEMS_SOLD)) +
  geom_col() +
  ggtitle(
    "Country 3 Monthly Sales Distribution"
  )

Overall the monthly trend for the 3 country seems to reflect the trend for daily. We can still see that C2 and C3 are sparse in terms of distribution.

Department and Country Sales Impact

by_country_data = raw_data %>% 
  group_by(DATE,COUNTRY) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  )

by_country_data %>% 
  group_by(COUNTRY) %>% 
  normality()

# A tibble: 3 × 5
  variable       COUNTRY statistic  p_value sample
  <chr>          <chr>       <dbl>    <dbl>  <dbl>
1 NUM_ITEMS_SOLD C1          0.561 2.36e-31    428
2 NUM_ITEMS_SOLD C2          0.566 1.48e-15     98
3 NUM_ITEMS_SOLD C3          0.561 1.63e-18    137

We’ve verified that our data for each group are not normal, we can proceed applying Kruskal-Wallis test.

by_country_data %>% 
  ggplot(aes(x = COUNTRY, y = NUM_ITEMS_SOLD)) +
  geom_boxplot()

by_country_data %>% 
  ggplot(aes(x = NUM_ITEMS_SOLD, fill = COUNTRY)) +
  geom_histogram( alpha = 0.5)

impact_analysis = kruskal.test(NUM_ITEMS_SOLD ~ COUNTRY, data = by_country_data)

impact_analysis

    Kruskal-Wallis rank sum test

data:  NUM_ITEMS_SOLD by COUNTRY
Kruskal-Wallis chi-squared = 93.447, df = 2, p-value < 2.2e-16

ggbetweenstats(
  data = by_country_data,
  x = COUNTRY,
  y = NUM_ITEMS_SOLD,
  type = "nonparametric"
)

We used Kruskal-Wallis test of independence to determine if number of items sold on different countries has a significant difference on their center measures. On the plot above, we can see that our p.value is 5.11 * 10^21 which indicates a really small number. This shows that at least one of our country displays a significant difference in central measures to other country. The other part of the plot shows a p.value calculated from a post-hoc analysis to identify what specific groups on our data has significant difference on each other. We can see the C3 is both significantly different in terms of center measure to C1 and C2.

by_department = raw_data %>% 
  group_by(DATE, DEPARTMENT) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD)
  )

We need to check first if our data is not normally distributed before conducting Mann-Whitney U test.

by_department %>% 
  group_by(DEPARTMENT) %>% 
  normality(NUM_ITEMS_SOLD)

# A tibble: 2 × 5
  variable       DEPARTMENT statistic  p_value sample
  <chr>          <chr>          <dbl>    <dbl>  <dbl>
1 NUM_ITEMS_SOLD Dep1           0.608 1.56e-30    446
2 NUM_ITEMS_SOLD Dep2           0.615 1.09e-15    112

by_department %>% 
  ggplot(aes(x = DEPARTMENT, y = NUM_ITEMS_SOLD)) +
  geom_boxplot()

by_department %>% 
  ggplot(aes(x = NUM_ITEMS_SOLD, fill = DEPARTMENT))+
  geom_histogram(alpha = 0.7)

ggbetweenstats(
  data = by_department,
  x = DEPARTMENT,
  y = NUM_ITEMS_SOLD,
  type = "nonparametric"
)

interpret_rank_biserial(-0.56)

[1] "very large"
(Rules: funder2019)

The analysis above shows that there is statistically significant large difference between Department 1 and 2. We were able to prove this difference by using Mann-Whitney U test for difference in Rank Sum. To determine how large the difference is, we use the rank biserial coefficient correlation for interpretation of the difference.

Model Building

Now that we have a sense of how our data looks like. We can start exploring models for forecasting

#Prepare data for the required granularity
country_demand_data = raw_data %>% 
  group_by(DATE, COUNTRY, SKU) %>% 
  reframe(
    NUM_ITEMS_SOLD = sum(NUM_ITEMS_SOLD, na.rm = TRUE)
  )

Splitting Data Set

Let’s create a function that will do the train and test split for each Country and SKU combination.

split_train_test = function(raw_data){
  
  test = tail(raw_data, n = 7)
  train = head(raw_data, n = nrow(raw_data) - 7)
  train_test = c()
  train_test$train = train
  train_test$test = test
  
  return(train_test)
  
}

Time Series Model

We will explore a time series model called STLF or Seasonal and Trend decomposition using Loess Forecasting. We use this model for our time series data because our data is an intermittent type of time series data. Intermittent type of time series data are characterized by irregular intervals of data points. We could resort to different missing value approximation for our data set but that can lead to bias and wrongful assumptions of our prediction. The STLF can handle both outliers and irregular placed data points for forecasting.

We will create a function that will accept our raw data and will process the train and test split, modeling, and metric evaluation

country = "C3"
sku = "SKU_110"
raw_data = country_demand_data
stlf_modelling = function(raw_data, country, sku){
  filtered_data = raw_data %>% 
    filter(COUNTRY == country) %>% 
    filter(SKU == sku)
  

  data_split = split_train_test(filtered_data)
  train = data_split$train
  test = data_split$test
   train_nrows = nrow(train)
   if(train_nrows <= 20) {
    return(message("ERROR:Rows are too few to create time series model"))
  }
  
  
  ts_data <- ts(train$NUM_ITEMS_SOLD, frequency = 7)
  stlf_model <- stlf(ts_data)
  stlf_model$mean <- pmax(stlf_model$mean,0) 
  predicted = forecast::forecast(stlf_model,h = nrow(test))
  
  train = train %>% 
  mutate(
    Key = "Trained_Data"
  )
  test = test %>% 
    mutate(
      Key = "Actual"
    )
  
  predicted = tibble(
    DATE = test$DATE,
    COUNTRY = country,
    SKU = sku,
    NUM_ITEMS_SOLD =predicted$mean %>% as.numeric() %>% round(),
    Key = "Predicted"
  )
  
  forecast_data = bind_rows(train,test)
  forecast_data = bind_rows(forecast_data,predicted)

   # truncate at zero

  return_list = c()
  return_list$forecast_data = forecast_data
  return_list$ts_model = stlf_model
  return(return_list)
}

time_series_model_result = stlf_modelling(country_demand_data,"C1","SKU_66")
time_series_model_result$forecast_data %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  )+
    scale_color_tq() +
  scale_fill_tq() 

Get Evaluation Metrics

time_series_model_result$ts_model$model %>% summary()

ETS(A,N,N) 

Call:
 ets(y = na.interp(x), model = etsmodel, allow.multiplicative.trend = allow.multiplicative.trend) 

  Smoothing parameters:
    alpha = 1e-04 

  Initial states:
    l = 1.493 

  sigma:  1.7864

     AIC     AICc      BIC 
140.7793 141.7023 144.9828 

Training set error measures:
                       ME     RMSE      MAE      MPE     MAPE      MASE
Training set 0.0005600939 1.725854 1.366538 94.11672 588.2326 0.6482282
                   ACF1
Training set 0.07158276

RMSE

actual = time_series_model_result$forecast_data %>% 
  filter(Key == "Actual") %>% 
  pull(NUM_ITEMS_SOLD)

predicted = time_series_model_result$forecast_data %>% 
  filter(Key == "Predicted") %>% 
  pull(NUM_ITEMS_SOLD)


sqrt(mean((actual - predicted)^2))

[1] 1

Metrics::rmse(actual,predicted)

[1] 1

Now let’s apply it on all combination of Country and SKU

#Get all Country values
countries = country_demand_data$COUNTRY %>% unique()

consolidated_stlf_result = c()
for (country in countries){
  #Filter country
  print(paste("Country: ",country," Starting ----*",sep = ""))
  filtered_country = country_demand_data %>% 
    filter(COUNTRY == country)
  
  sku_list = filtered_country$SKU %>% unique()
  sku_temp_data = c()
  for(sku in sku_list){
    ts_model_result = stlf_modelling(filtered_country, country = country, sku = sku)
    if(is.null(ts_model_result)){
      print("Output Null, skipping ------->")
      next
    }
    sku_temp_data[[sku]] = ts_model_result$forecast_data
    print(paste("SKU: ",sku," done!",sep = ""))
  }
  sku_temp_data = bind_rows(sku_temp_data)
  consolidated_stlf_result[[country]] = sku_temp_data
}

[1] "Country: C1 Starting ----*"
[1] "SKU: SKU_15 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_76 done!"
[1] "SKU: SKU_89 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_46 done!"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_38 done!"
[1] "SKU: SKU_55 done!"
[1] "SKU: SKU_13 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_33 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_69 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_67 done!"
[1] "SKU: SKU_75 done!"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_66 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_94 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_79 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_95 done!"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_72 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_6 done!"
[1] "SKU: SKU_77 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Country: C3 Starting ----*"
[1] "SKU: SKU_110 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_118 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_115 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Country: C2 Starting ----*"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "SKU: SKU_66 done!"
[1] "SKU: SKU_63 done!"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"
[1] "Output Null, skipping ------->"

Let’s check some of SKU model from each Country

Country 1

consolidated_stlf_result$C1 %>% 
  filter(SKU == "SKU_66") %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

consolidated_stlf_result$C1 %>% 
  filter(SKU == "SKU_89") %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

consolidated_stlf_result$C1 %>% 
  filter(SKU == "SKU_33") %>% 
  filter(DATE > as.Date("2021-11-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

Country 2

consolidated_stlf_result$C2 %>% 
  filter(SKU == "SKU_66") %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

consolidated_stlf_result$C2 %>% 
  filter(SKU == "SKU_63") %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

Country 3

consolidated_stlf_result$C3 %>% 
  filter(SKU == "SKU_110") %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

consolidated_stlf_result$C3 %>% 
  filter(SKU == "SKU_118") %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

consolidated_stlf_result$C3 %>% 
  filter(SKU == "SKU_115") %>% 
  filter(DATE > as.Date("2022-10-01")) %>% 
  ggplot(aes(x = DATE,y = NUM_ITEMS_SOLD, color = Key)) +
           geom_line()+
  ggtitle(
    "Seasonal and Trend decomposition using Loess Forecasting (STLF) Result"
  ) +
    scale_color_tq() +
  scale_fill_tq() 

Model Evaluation

We are going to use Root mean squared Error to evaluate our model.

stlf_country_model_metrics = c()
for(country in countries){
  country_data = consolidated_stlf_result[[country]]
  
  sku_model_list = unique(country_data$SKU)
  sku_metric_val = c()
  for(sku in sku_model_list){
    
    filtered_model_data = country_data %>% 
      filter(SKU == sku)
    
    actual = filtered_model_data %>% 
      filter(Key == "Actual") %>% 
      pull(NUM_ITEMS_SOLD)
    predicted = filtered_model_data %>% 
      filter(Key == "Predicted") %>% 
      pull(NUM_ITEMS_SOLD)
    
    rmse_val = Metrics::rmse(actual,predicted)
    result_tibble = tibble(
      SKU = sku,
      RMSE = rmse_val
    )
    sku_metric_val[[sku]] = result_tibble
  }
  
  sku_metric_val = bind_rows(sku_metric_val)
  stlf_country_model_metrics[[country]] =sku_metric_val
}

stlf_country_model_metrics = bind_rows(stlf_country_model_metrics,.id = "COUNTRY")

stlf_country_model_metrics %>% 
  filter(COUNTRY == "C1") %>% 
  filter(SKU %in% c("SKU_66","SKU_89","SKU_33"))

# A tibble: 3 × 3
  COUNTRY SKU     RMSE
  <chr>   <chr>  <dbl>
1 C1      SKU_89 1.41 
2 C1      SKU_33 0.655
3 C1      SKU_66 1    

We can try to check the average RMSE of each Country for us to have an overview of our model performance.

stlf_country_model_metrics %>% 
  group_by(COUNTRY) %>% 
  reframe(
    Average_RMSE = mean(RMSE, na.rm = TRUE)
  )

# A tibble: 3 × 2
  COUNTRY Average_RMSE
  <chr>          <dbl>
1 C1             4.24 
2 C2             0.926
3 C3            75.1  

Regression Model

We will use another type of model to predict our daily sales data. A more common approach called Linear Regression. Unlike the common linear regression type, we are going to use a much more time-series data aligned type of regression modelling. This model is more commonly known as ARIMA or Auto Regressive Integrated Moving Average model.

c1_demand = country_demand_data %>% 
  filter(COUNTRY == "C1") %>% 
  filter(SKU == "SKU_66") %>% 
  mutate(
    DATE = as.Date(DATE)
  )

sales_split = time_series_split(
  c1_demand,
  assess = "14 days",
  cumulative = TRUE
)



sales_split %>% 
  tk_time_series_cv_plan() %>% 
  plot_time_series_cv_plan(DATE,NUM_ITEMS_SOLD)

model_arima_sales = arima_reg() %>% 
  set_engine("auto_arima") %>% 
  fit(NUM_ITEMS_SOLD ~ DATE, training(sales_split))

model_arima_sales

parsnip model object

Series: outcome 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
        mean
      1.4516
s.e.  0.3439

sigma^2 = 3.789:  log likelihood = -64.13
AIC=132.25   AICc=132.68   BIC=135.12

model_tbl_sales = modeltime_table(
  model_arima_sales
)

calib_tbl_sales = model_tbl_sales %>% 
  modeltime_calibrate(testing(sales_split))


calib_tbl_sales %>% modeltime_accuracy() %>% 
  select(`.model_id`,`.model_desc`,rmse)

# A tibble: 1 × 3
  .model_id .model_desc                      rmse
      <int> <chr>                           <dbl>
1         1 ARIMA(0,0,0) WITH NON-ZERO MEAN 0.950

calib_tbl_sales %>% 
  modeltime_forecast(
    new_data = testing(sales_split),
    actual_data = c1_demand
  ) %>% 
  plot_modeltime_forecast()

We can then wrap the steps we just did into a function and run it in a loop for each combination of Country and SKU

arima_modelling = function(raw_data, country, sku){
  cntry_demand = raw_data %>% 
  filter(COUNTRY == country) %>% 
  filter(SKU == sku) %>% 
  mutate(
    DATE = as.Date(DATE)
  )
  
  if(nrow(cntry_demand) <=25 ){
    return(message("Error:There should be at least 6 nrows in `data`"))
  }

sales_split =   tryCatch(
  {time_series_split(
  cntry_demand,
  assess = "14 days",
  cumulative = TRUE
)},
error = function(e){
 NULL
}
  )

if(is.null(sales_split)){
  return(message("Data Invalid for assessment split"))
}



model_arima_sales = arima_reg() %>% 
  set_engine("auto_arima") %>% 
  fit(NUM_ITEMS_SOLD ~ DATE, training(sales_split))


model_tbl_sales = modeltime_table(
  model_arima_sales
)

calib_tbl_sales = model_tbl_sales %>% 
  modeltime_calibrate(testing(sales_split))

evals = calib_tbl_sales %>% modeltime_accuracy()

if("rmse" %in% colnames(evals)){
  evals = evals %>% 
  select(`.model_id`,`.model_desc`,rmse)
}else{
  return(message("Data not valid for Evaluation!"))
}

forecast_plot = calib_tbl_sales %>% 
  modeltime_forecast(
    new_data = testing(sales_split),
    actual_data = cntry_demand
  ) %>% 
  plot_modeltime_forecast()


return_list = c()
return_list$plot = sales_split %>% 
  tk_time_series_cv_plan() %>% 
  plot_time_series_cv_plan(DATE,NUM_ITEMS_SOLD)
return_list$model = model_arima_sales
return_list$eval_metrics = evals
return_list$forecast_plot = forecast_plot



return(return_list)

}

sampel_result = arima_modelling(country_demand_data, "C1","SKU_66")

sampel_result$plot

sampel_result$model

parsnip model object

Series: outcome 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
        mean
      1.4516
s.e.  0.3439

sigma^2 = 3.789:  log likelihood = -64.13
AIC=132.25   AICc=132.68   BIC=135.12

sampel_result$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc                      rmse
      <int> <chr>                           <dbl>
1         1 ARIMA(0,0,0) WITH NON-ZERO MEAN 0.950

sampel_result$forecast_plot

Awesome! We were able to wrap the steps we previously did inside a function. Now we just need to loop it on all SKU under each Country

consolidated_arima_result = c()
for(country in countries){
  print(paste("Country: ",country," Starting ----*",sep = ""))
  filtered_data = country_demand_data %>% 
    filter(COUNTRY == country)
  sku_list = filtered_data$SKU %>% unique()
  consolidated_sku_models = c()
  for(sku in sku_list){
    arima_model_result = arima_modelling(country_demand_data,country = country,sku = sku)
    if(is.null(arima_model_result)) {
      next
    }
    consolidated_sku_models[[sku]] = arima_model_result
    print(paste("SKU: ",sku," done!",sep = ""))
  }
  
  consolidated_arima_result[[country]] = consolidated_sku_models
}

[1] "Country: C1 Starting ----*"
[1] "SKU: SKU_15 done!"
[1] "SKU: SKU_76 done!"
[1] "SKU: SKU_89 done!"
[1] "SKU: SKU_46 done!"
[1] "SKU: SKU_38 done!"
[1] "SKU: SKU_55 done!"
[1] "SKU: SKU_13 done!"
[1] "SKU: SKU_33 done!"
[1] "SKU: SKU_69 done!"
[1] "SKU: SKU_67 done!"
[1] "SKU: SKU_66 done!"
[1] "SKU: SKU_94 done!"
[1] "SKU: SKU_79 done!"
[1] "SKU: SKU_95 done!"
[1] "SKU: SKU_72 done!"
[1] "SKU: SKU_22 done!"
[1] "SKU: SKU_6 done!"
[1] "SKU: SKU_77 done!"
[1] "Country: C3 Starting ----*"
[1] "SKU: SKU_110 done!"
[1] "SKU: SKU_115 done!"
[1] "Country: C2 Starting ----*"
[1] "SKU: SKU_66 done!"

Now that we are done running our ARIMA model for all possible SKU under each Country, let’s look at the same SKU we visualize from the STLF model and see the difference on their forecast.

Country 1

consolidated_arima_result$C1$SKU_66$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc                      rmse
      <int> <chr>                           <dbl>
1         1 ARIMA(0,0,0) WITH NON-ZERO MEAN 0.950

consolidated_arima_result$C1$SKU_66$forecast_plot

consolidated_arima_result$C1$SKU_89$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc                      rmse
      <int> <chr>                           <dbl>
1         1 ARIMA(0,0,0) WITH NON-ZERO MEAN  1.00

consolidated_arima_result$C1$SKU_89$forecast_plot

consolidated_arima_result$C1$SKU_33$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc                      rmse
      <int> <chr>                           <dbl>
1         1 ARIMA(0,0,0) WITH NON-ZERO MEAN 0.559

consolidated_arima_result$C1$SKU_33$forecast_plot

Looks like for Country 1, ARIMA model seems to perform well compared to STLF model.

The last model we are going to explore is called Prophet Model (by Facebook). The prophet model is a tool for estimating open source time series data with daily observations and strong seasonal patterns such as holidays. It is very solid and versatile as it allows its users to easily include such elements as seasonality, holidays and other special events in their forecasting assets. Prophet is ideal for handling irregularly spaced time series problems because it can automatically detect missing points, effectively treat outliers and still give accurate forecasts. The inclusion of trend, seasonality and holiday-based components in the Prophet model structure enables this tool to capture many different intermittent patterns which makes it a widely applicable choice for many forecasting scenarios where analytical complexity may be a problem.

prophet_modelling = function(raw_data, country, sku){
  
  
    cntry_demand = raw_data %>% 
  filter(COUNTRY == country) %>% 
  filter(SKU == sku) %>% 
  mutate(
    DATE = as.Date(DATE)
  )
  
  if(nrow(cntry_demand) <=25 ){
    return(message("Error:There should be at least 6 nrows in `data`"))
  }

sales_split =   tryCatch(
  {time_series_split(
  cntry_demand,
  assess = "14 days",
  cumulative = TRUE
)},
error = function(e){
 NULL
}
  )

if(is.null(sales_split)){
  return(message("Data Invalid for assessment split"))
}

model_prophet_sales = prophet_reg(seasonality_daily =  TRUE) %>% 
  set_engine("prophet") %>% 
  fit(NUM_ITEMS_SOLD ~ DATE, training(sales_split))

model_prophet_sales

model_tbl_sales = modeltime_table(
  model_prophet_sales
)

calib_tbl_sales = model_tbl_sales %>% 
  modeltime_calibrate(testing(sales_split))

evals = calib_tbl_sales %>% modeltime_accuracy()

if("rmse" %in% colnames(evals)){
  evals = evals %>% 
  select(`.model_id`,`.model_desc`,rmse)
}else{
  return(message("Data not valid for Evaluation!"))
}

forecast_plot = calib_tbl_sales %>% 
  modeltime_forecast(
    new_data = testing(sales_split),
    actual_data = cntry_demand
  ) %>% 
  plot_modeltime_forecast()


return_list = c()
return_list$plot = sales_split %>% 
  tk_time_series_cv_plan() %>% 
  plot_time_series_cv_plan(DATE,NUM_ITEMS_SOLD)
return_list$model = model_arima_sales
return_list$eval_metrics = evals
return_list$forecast_plot = forecast_plot

return(return_list)

}

prophet_modelling(country_demand_data, "C1","SKU_66")

$plot

$model
parsnip model object

Series: outcome 
ARIMA(0,0,0) with non-zero mean 

Coefficients:
        mean
      1.4516
s.e.  0.3439

sigma^2 = 3.789:  log likelihood = -64.13
AIC=132.25   AICc=132.68   BIC=135.12

$eval_metrics
# A tibble: 1 × 3
  .model_id .model_desc  rmse
      <int> <chr>       <dbl>
1         1 PROPHET      1.60

$forecast_plot

Let’s apply the custom function for prophet on loop

consolidated_prophet_result = c()
for(country in countries){
  print(paste("Country: ",country," Starting ----*",sep = ""))
  filtered_data = country_demand_data %>% 
    filter(COUNTRY == country)
  sku_list = filtered_data$SKU %>% unique()
  consolidated_sku_models = c()
  for(sku in sku_list){
    prophet_model_result = prophet_modelling(country_demand_data,country = country,sku = sku)
    if(is.null(prophet_model_result)) {
      next
    }
    consolidated_sku_models[[sku]] = prophet_model_result
    print(paste("SKU: ",sku," done!",sep = ""))
  }
  
  consolidated_prophet_result[[country]] = consolidated_sku_models
}

[1] "Country: C1 Starting ----*"
[1] "SKU: SKU_15 done!"
[1] "SKU: SKU_76 done!"
[1] "SKU: SKU_89 done!"
[1] "SKU: SKU_46 done!"
[1] "SKU: SKU_38 done!"
[1] "SKU: SKU_55 done!"
[1] "SKU: SKU_13 done!"
[1] "SKU: SKU_33 done!"
[1] "SKU: SKU_69 done!"
[1] "SKU: SKU_67 done!"
[1] "SKU: SKU_66 done!"
[1] "SKU: SKU_94 done!"
[1] "SKU: SKU_79 done!"
[1] "SKU: SKU_95 done!"
[1] "SKU: SKU_72 done!"
[1] "SKU: SKU_22 done!"
[1] "SKU: SKU_6 done!"
[1] "SKU: SKU_77 done!"
[1] "Country: C3 Starting ----*"
[1] "SKU: SKU_110 done!"
[1] "SKU: SKU_115 done!"
[1] "Country: C2 Starting ----*"
[1] "SKU: SKU_66 done!"

consolidated_prophet_result$C1$SKU_66$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc  rmse
      <int> <chr>       <dbl>
1         1 PROPHET      1.60

consolidated_prophet_result$C1$SKU_66$forecast_plot

consolidated_prophet_result$C1$SKU_89$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc  rmse
      <int> <chr>       <dbl>
1         1 PROPHET      1.40

consolidated_prophet_result$C1$SKU_89$forecast_plot

consolidated_prophet_result$C1$SKU_33$eval_metrics

# A tibble: 1 × 3
  .model_id .model_desc  rmse
      <int> <chr>       <dbl>
1         1 PROPHET     0.614

consolidated_prophet_result$C1$SKU_33$forecast_plot

The Prophet Model seems to be the worst performing model among the 3. To have a better view let’s tabulate the result from the 3 model

final_c1_metrics = tibble(
  Country = c("C1","C1","C1"),
  SKU = c("SKU_66","SKU_89","SKU_33"),
  STLF = c(1,1.41,0.65),
  ARIMA = c(0.95,1,0.559),
  PROPHET = c(1.60,1.40,0.614)
)
final_c1_metrics

# A tibble: 3 × 5
  Country SKU     STLF ARIMA PROPHET
  <chr>   <chr>  <dbl> <dbl>   <dbl>
1 C1      SKU_66  1    0.95    1.6  
2 C1      SKU_89  1.41 1       1.4  
3 C1      SKU_33  0.65 0.559   0.614

As conclusion, ARIMA seems to perform well compared to the other 2 models. It is followed by STLF then on the last spot is Prophet.

Model Deployment

For model deployment, we can create an API end point with the best model on the back end that will return forecast based on inputted number of days, country and SKU. If we will continue to use R on the model, we can use Plumber package to transform our script into an API end point and deploy it in either AWS or Azure cloud platform. We can also convert our script into a python script and export a Pickle file that can be use as backend calculator for creating a Flask Web API. We can also upload the API on AWS, Azure or GCP.

Conclusion and Recommendation

The data we used for predictions are highly sporadic and intermittent that’s why it poses a challenge to come up with a meaningful prediction. Here are my recommendations

1. We can reach out to data steward or data management and ask what are some of the rules on data imputation. By applying data imputation we can convert our data into a more cleaner data with no missing data point interval.
2. Implement rules for SKU on forecasting maturity. This could be a task for Data Engineering team, but I will highly recommend for us to set a threshold for number of rows that are valid before we can create a model for a specific SKU. A lot of SKU from C1,C2 and C3 were rejected because of low row count and low values variance.
3. Based on the 3 model, we can use ARIMA to have a confident estimation of our demand for each SKU on each country on a daily granularity.
4. We can also explore the possibility of a monthly prediction model to eliminate noise present on data with daily granularity.