library(dplyr)
library(lubridate)
library(ggplot2)
library(prophet)

Background

In this exercise we will create a time series model using a sales dataset provided in Kaggle. The dataset you will be using is stored under data_input/dive-deeper folder in a file called train.csv. This dataset consist of 4 years of store-item sales data. Along this dive deeper exercise we will perform a time series analysis and forecasting using Prophet model.

Data Preparation

In this section we will prepare our raw dataset into a proper time series format. Please read in our train.csv dataset using read.csv() function we have learned and perform a simple inspection of the data:

sales_t <- read.csv("data_input/dive-deeper/train.csv")

# glimpse, head, tail,
glimpse(sales_t)

## Observations: 730,500
## Variables: 4
## $ date  <fct> 2013-01-01, 2013-01-02, 2013-01-03, 2013-01-04, 2013-01-05, 201…
## $ store <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ item  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ sales <int> 13, 11, 14, 13, 10, 12, 10, 9, 12, 9, 9, 7, 10, 12, 5, 7, 16, 7…

head(sales_t)

##         date store item sales
## 1 2013-01-01     1    1    13
## 2 2013-01-02     1    1    11
## 3 2013-01-03     1    1    14
## 4 2013-01-04     1    1    13
## 5 2013-01-05     1    1    10
## 6 2013-01-06     1    1    12

tail(sales_t)

##              date store item sales
## 730495 2016-12-26    10   50    61
## 730496 2016-12-27    10   50    60
## 730497 2016-12-28    10   50    43
## 730498 2016-12-29    10   50    68
## 730499 2016-12-30    10   50    63
## 730500 2016-12-31    10   50    64

# Your code here
class(sales_t$date)

## [1] "factor"

There are 730,500 observations with 4 variables, and the data is about sales information include date, where the store is, and what kind of item it is in each row. the class of date needs to be change to a date type.

# Your code here

# import data

# Transform the time stamp column into date type 
sales_t <- sales_t %>% 
  mutate(date= ymd(date))

# Decide the frequency of sampling (hourly / daily / weekly / monthly / etc.) 

# Perform a data aggregation based on the sampling frequency using `group_by()` and `summarise()` 
daily_demand <- sales_t %>% 
  group_by(date) %>% 
  summarise(demand = sum(sales))
daily_demand

## # A tibble: 1,461 x 2
##    date       demand
##    <date>      <int>
##  1 2013-01-01  13696
##  2 2013-01-02  13678
##  3 2013-01-03  14488
##  4 2013-01-04  15677
##  5 2013-01-05  16237
##  6 2013-01-06  17291
##  7 2013-01-07  11769
##  8 2013-01-08  13560
##  9 2013-01-09  13878
## 10 2013-01-10  14642
## # … with 1,451 more rows

# Sort the data frame based on its time of sample

Visual Exploratory Analysis

# Your code here
daily_demand %>% 
  ggplot(aes(x=date, y=demand)) + 
  geom_line(color = "pink") + 
  labs( 
    title = "Daily Software Sales", 
    subtitle = "Store 31", 
    caption = "1C Company", 
    x = "Date", 
    y = "Total Sales" 
  ) + 
  theme_minimal()

Baseline Prophet Model

# Your code here
train_daily <- daily_demand %>% 
  rename(
    ds = "date",
    y = "demand"
    )

model_prophet <- prophet() %>%
  fit.prophet(train_daily)

future_prophet <- make_future_dataframe(model_prophet, periods = 365, freq = "day")

# Your code here
forecast_prophet <- predict(model_prophet, future_prophet)

prophet_plot_components(model_prophet, forecast_prophet)

We can see in the weekly seasonality that Monday is the lowest point of demand. While yearly seasonality, July resulting the highest point through the year. This is a bit confusing because on the information before, we know that the demand is high around Nov to Dec. This initial model coundn’t catch the seasonality. We should tune the model to see if the additional seasonality (for instance, monthly, quarterly) can affect the model to sensitively catch the seasonality better.

Model Fine Tuning

Trend (Linear)
- changepoint.prior.scale
- changepoints =
Seasonality
- Daily, Weekly, Yearly (TRUE/FALSE)
- yearly.seasonality = FALSE
- daily.seasonality = FALSE –> Per jam
- weekly.seasonality = FALSE
- Add non-regular seasonality
  - period
  - fourier.order (3-10)
Holiday

holiday <- 
  data.frame( 
    holiday = "newyeareve", 
    ds = dmy(c("31-12-2013", "31-12-2014", "31-12-2015", "31-12-2016")), 
    lower_window = -5, 
    upper_window = 0 
  ) 
holiday

##      holiday         ds lower_window upper_window
## 1 newyeareve 2013-12-31           -5            0
## 2 newyeareve 2014-12-31           -5            0
## 3 newyeareve 2015-12-31           -5            0
## 4 newyeareve 2016-12-31           -5            0

# Your code here
train_daily <- daily_demand %>% 
  rename(
    ds = "date",
    y = "demand"
    )

model_new <- prophet(yearly.seasonality = TRUE, 
                             changepoint.prior.scale = 0.5,
                     changepoints = ymd("2014-01-01"),
                     holidays = holiday) %>%
  add_seasonality(name = "monthly", period = 30.5, fourier.order = 5) %>% 
  fit.prophet(train_daily)

future_new <- make_future_dataframe(model_new, periods = 365, freq = "day")

forecast_new <- predict(model_new, future_new)

prophet_plot_components(model_new, forecast_new)

Model Evaluation

test <- read.csv("data_input/dive-deeper/test.csv")
glimpse(test)

## Observations: 182,500
## Variables: 4
## $ date  <fct> 2017-01-01, 2017-01-02, 2017-01-03, 2017-01-04, 2017-01-05, 201…
## $ store <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ item  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ sales <int> 19, 15, 10, 16, 14, 24, 14, 20, 18, 11, 14, 17, 7, 16, 29, 15, …

# Your code here
test <- test %>% 
  mutate(date= ymd(date))

test_ <- test %>% 
  group_by(date) %>% 
  summarise(demand = sum(sales))

test_data <- test_ %>% rename( 
    "ds" = date, 
    "y" = demand 
  )

evalx <- test_data %>% 
  mutate( 
    ds = as.POSIXct(ds) 
  ) %>% 
  left_join(forecast_prophet) %>% 
  select(ds, y, yhat, yhat_upper, yhat_lower) 
evalx

## # A tibble: 365 x 5
##    ds                      y   yhat yhat_upper yhat_lower
##    <dttm>              <int>  <dbl>      <dbl>      <dbl>
##  1 2017-01-01 00:00:00 23709 26755.     28089.     25334.
##  2 2017-01-02 00:00:00 15772 16667.     18023.     15236.
##  3 2017-01-03 00:00:00 18650 19998.     21405.     18716.
##  4 2017-01-04 00:00:00 18510 20017.     21442.     18627.
##  5 2017-01-05 00:00:00 19895 21598.     23003.     20139.
##  6 2017-01-06 00:00:00 20994 23282.     24511.     21790.
##  7 2017-01-07 00:00:00 22591 24913.     26331.     23511.
##  8 2017-01-08 00:00:00 23700 26610.     27901.     25092.
##  9 2017-01-09 00:00:00 15797 16508.     17898.     15023.
## 10 2017-01-10 00:00:00 18608 19834.     21142.     18423.
## # … with 355 more rows

eval <- test_data %>% 
  mutate( 
    ds = as.POSIXct(ds) 
  ) %>% 
  left_join(forecast_new) %>% 
  select(ds, y, yhat, yhat_upper, yhat_lower) 
eval

## # A tibble: 365 x 5
##    ds                      y   yhat yhat_upper yhat_lower
##    <dttm>              <int>  <dbl>      <dbl>      <dbl>
##  1 2017-01-01 00:00:00 23709 26974.     28498.     25411.
##  2 2017-01-02 00:00:00 15772 16911.     18426.     15405.
##  3 2017-01-03 00:00:00 18650 20290.     21944.     18715.
##  4 2017-01-04 00:00:00 18510 20351.     21910.     18761.
##  5 2017-01-05 00:00:00 19895 21971.     23539.     20386.
##  6 2017-01-06 00:00:00 20994 23670.     25157.     22137.
##  7 2017-01-07 00:00:00 22591 25295.     26833.     23718.
##  8 2017-01-08 00:00:00 23700 26975.     28600.     25423.
##  9 2017-01-09 00:00:00 15797 16861.     18321.     15366.
## 10 2017-01-10 00:00:00 18608 20211.     21696.     18700.
## # … with 355 more rows

# Your code here
mape <- function(y, yhat) { 
  return(mean(abs(y - yhat)/ y)) 
} 
mape(eval$y, eval$yhat)

## [1] 0.05254495

mape(evalx$y, evalx$yhat)

## [1] 0.05550715

The MAPE on forecast_new return 0.05254495 which is slightly better than the forecast_prophet on 0.05550715. It means the tuning model to add seasonality on monthly and adding Holidays argument is quite significant to make the model perform better than without tuning.

# Your code here
eval %>% 
  mutate( 
    month = month(ds, label = T) 
  ) %>% 
  group_by(month) %>% 
  summarise( 
    mape = mean(abs(y-yhat)/y) 
  ) %>% 
  ggplot(aes(x = month, y = mape)) + 
  geom_line(group=1) + 
  geom_point() + 
  theme_minimal()

eval %>% 
  mutate( 
    flag_extreme = ifelse(y > yhat_upper | y < yhat_lower, "1", "0"), 
    month = month(ds, label = T), 
  ) %>% 
ggplot(aes(x=ds, y=y)) + 
  geom_ribbon(aes(ymin=yhat_lower, ymax=yhat_upper), alpha = 0.3, fill="dodgerblue4") + 
  geom_line(aes(x=ds, y=yhat)) + 
  geom_point(aes(color=flag_extreme), alpha = 0.5, size = 2) + 
  facet_wrap(~month, scales = "free_x")

The seasonality of monthly added to the model_new shows a bigger MAPE value on Jan, slightly on July, and increasing during the last Quarter of the year. This seasonality gives us a hint that the sales could have some outliers when the demand are high. With this insight, we can tell to the company they should be well-prepared with the stock so they can be ready for increasing demand on Jan, Feb, Jul, Nov, and Dec each year.

It is interesting that a model can recognize seasonality and even better: it is flexible enough to be tuned on the seasonality and adding a changepoint, (also Holiday!), where there will always be changing behaviour over time. The Prophet package can be a good tool to show business people understand better what can seasonality do to our data, and also to our target (forecasted sales). It can also make a business run to be more effective as the forecasted data could be a valuable information to prepare the warehouse and the worker to catching on the demand.

Business Report: Store Demand Forecasting

Andina Septia

May 14, 2020