Chapter 6 Project 4: Demand Forecasting — Prophet + XGBoost Ensemble
Accurate demand forecasting drives inventory decisions, supply chain optimization, and revenue planning. In this chapter, we build an ensemble forecasting system that combines Facebook Prophet’s ability to model seasonality and trends with XGBoost’s power to capture complex feature interactions.
6.1 Business Context
We need to forecast daily unit sales for each SKU-warehouse combination 30 days ahead. The system must:
- Handle multiple seasonalities (weekly, yearly, holiday effects)
- Incorporate promotional calendars
- Account for stockouts (missing data that isn’t truly zero demand)
- Provide safety stock recommendations at 95% service level
6.2 Architecture Overview
Sales Data (3 years)
│
├──► Prophet ──► Trend + Seasonality ──┐
│ ├──► Weighted Ensemble ──► Forecast + Safety Stock
└──► XGBoost ──► Lag/Rolling Features ──┘6.3 Data Pipeline
library(prophet)
library(xgboost)
library(timetk)
library(modeltime)
library(tidyverse)
library(lubridate)
get_sales_data <- function(con, sku_id = NULL, warehouse_id = NULL) {
query <- tbl(con, "fact_daily_sales") %>%
filter(date >= Sys.Date() - years(3))
if (!is.null(sku_id)) query <- query %>% filter(sku_id == !!sku_id)
if (!is.null(warehouse_id)) query <- query %>% filter(warehouse_id == !!warehouse_id)
query %>%
select(date, sku_id, warehouse_id, units_sold, revenue,
promotion_flag, stockout_flag) %>%
collect() %>%
mutate(date = as.Date(date))
}6.4 Feature Engineering for Time Series
XGBoost requires carefully constructed time-based features. We create lag features, rolling statistics, and cyclical encodings.
build_ts_features <- function(df) {
df %>%
group_by(sku_id, warehouse_id) %>%
arrange(date) %>%
mutate(
# Lag features: capture recent sales momentum
sales_lag_1 = lag(units_sold, 1),
sales_lag_7 = lag(units_sold, 7), # Same day last week
sales_lag_14 = lag(units_sold, 14), # Two weeks ago
sales_lag_28 = lag(units_sold, 28), # Four weeks ago
# Rolling statistics: smooth out noise
sales_roll_mean_7 = zoo::rollmean(units_sold, 7, fill = NA, align = "right"),
sales_roll_mean_14 = zoo::rollmean(units_sold, 14, fill = NA, align = "right"),
sales_roll_mean_30 = zoo::rollmean(units_sold, 30, fill = NA, align = "right"),
sales_roll_sd_7 = zoo::rollapply(units_sold, 7, sd, fill = NA, align = "right"),
# Seasonality features
day_of_week = wday(date, label = TRUE),
month = month(date, label = TRUE),
quarter = quarter(date),
is_month_start = as.integer(day(date) <= 5),
is_month_end = as.integer(day(date) >= 25),
# Cyclical encoding: maps linear features to circles
# This prevents the model from thinking January (1) and December (12) are far apart
sin_day = sin(2 * pi * wday(date) / 7),
cos_day = cos(2 * pi * wday(date) / 7),
sin_month = sin(2 * pi * month(date) / 12),
cos_month = cos(2 * pi * month(date) / 12),
# Promo interaction: days since last promotion
days_since_promo = as.numeric(date - lag(date[promotion_flag == 1], 1)),
# Target
units_sold = units_sold
) %>%
ungroup() %>%
drop_na()
}Tip: Cyclical encoding is crucial for time features. Without it, a model might learn that December (12) is “greater than” January (1), missing that they’re adjacent months.
6.5 Prophet Model
Prophet decomposes time series into trend, seasonality, and holiday effects.
fit_prophet <- function(df) {
# Define holidays with impact windows
holidays <- tibble(
holiday = c("Black Friday", "Cyber Monday", "Christmas", "New Year"),
ds = as.Date(c("2023-11-24", "2023-11-27", "2023-12-25", "2024-01-01")),
lower_window = c(-1, -1, -2, -1),
upper_window = c(1, 1, 2, 1)
)
prophet_df <- df %>%
select(ds = date, y = units_sold, promotion_flag) %>%
mutate(promotion_flag = as.numeric(promotion_flag))
m <- prophet(
yearly.seasonality = TRUE,
weekly.seasonality = TRUE,
daily.seasonality = FALSE,
changepoint.prior.scale = 0.05, # Flexibility of trend changes
seasonality.prior.scale = 10, # Strength of seasonality
holidays = holidays
)
m <- add_regressor(m, "promotion_flag")
m <- fit.prophet(m, prophet_df)
# Generate future dataframe
future <- make_future_dataframe(m, periods = 30)
future$promotion_flag <- get_promotion_calendar(future$ds)
forecast <- predict(m, future)
list(model = m, forecast = forecast)
}6.6 XGBoost with Time Series Cross-Validation
Standard k-fold CV leaks future information into training. Time series CV respects temporal ordering.
fit_xgboost_ts <- function(df) {
feature_cols <- c(
"sales_lag_1", "sales_lag_7", "sales_lag_14", "sales_lag_28",
"sales_roll_mean_7", "sales_roll_mean_14", "sales_roll_mean_30",
"sales_roll_sd_7", "sin_day", "cos_day", "sin_month", "cos_month",
"is_month_start", "is_month_end", "promotion_flag"
)
# Time-series cross-validation
tscv <- time_series_cv(
data = df,
date_var = date,
initial = "2 years",
assess = "30 days",
skip = "30 days",
cumulative = FALSE
)
dtrain <- xgb.DMatrix(
data = as.matrix(df %>% select(all_of(feature_cols))),
label = df$units_sold
)
params <- list(
objective = "reg:squarederror",
eta = 0.05,
max_depth = 8,
subsample = 0.8,
colsample_bytree = 0.8,
eval_metric = "mae"
)
xgb_cv <- xgb.cv(
params = params,
data = dtrain,
nrounds = 500,
nfold = 5,
early_stopping_rounds = 20,
print_every_n = 50
)
final_model <- xgb.train(
params = params,
data = dtrain,
nrounds = xgb_cv$best_iteration,
watchlist = list(train = dtrain)
)
list(model = final_model, features = feature_cols)
}Warning: Never use random k-fold CV for time series. It creates data leakage by training on future periods and testing on past periods, producing optimistically biased metrics.
6.7 Ensemble Forecast
We combine Prophet and XGBoost predictions, weighting by historical performance.
ensemble_forecast <- function(prophet_fc, xgb_fc, weights = c(0.4, 0.6)) {
ensemble <- weights[1] * prophet_fc + weights[2] * xgb_fc
# Safety stock at 95% service level
# Assumes forecast errors are normally distributed
forecast_std <- sd(c(prophet_fc, xgb_fc))
safety_stock <- qnorm(0.95) * forecast_std
tibble(
date = seq.Date(Sys.Date() + 1, by = "day", length.out = 30),
point_forecast = round(pmax(ensemble, 0)), # No negative sales
safety_stock = round(safety_stock),
reorder_point = round(pmax(ensemble + safety_stock - current_inventory, 0))
)
}6.8 Full Pipeline Function
forecast_sku <- function(con, sku_id, warehouse_id) {
df <- get_sales_data(con, sku_id, warehouse_id) %>%
build_ts_features()
# Prophet forecast
prophet_result <- fit_prophet(df)
prophet_fc <- tail(prophet_result$forecast$yhat, 30)
# XGBoost forecast
xgb_result <- fit_xgboost_ts(df)
xgb_fc <- predict(
xgb_result$model,
xgb.DMatrix(as.matrix(df %>% tail(30) %>% select(all_of(xgb_result$features))))
)
# Ensemble
ensemble_forecast(prophet_fc, xgb_fc)
}6.10 Exercises
- Implement dynamic weights for the ensemble based on recent model performance (e.g., weight by inverse MAE of last 30 days).
- Add external regressors to Prophet for weather data or competitor pricing.
- Build a
modeltimeworkflow that compares Prophet, XGBoost, ARIMA, and ETS automatically. - Implement a stockout correction mechanism that imputes likely demand during stockout periods.
Next, we’ll segment customers using RFM analysis and K-Means clustering to enable targeted marketing strategies.