Time series forecasting project to predict daily sales for a large Ecuadorian supermarket chain (Corporación Favorita), using the Kaggle Store Sales dataset.
- Problem: Forecast daily unit sales for 1,782 store × product-family series across 54 stores and 33 product families
- Result: RMSLE 0.381 (XGBoost), a ~28% improvement over a strong naive baseline
- Value: Reliable short-term forecasts for inventory planning, waste reduction, and promotion-scenario simulation
- Problem Statement
- Business Value
- Dataset
- Data Challenges and Transformations
- Exploratory Data Analysis
- Time Series Analysis
- Panel Data Analysis
- Methodology
- Feature Engineering
- Models and Hyperparameter Optimisation
- Results and Evaluation
- Explainability — SHAP Analysis
- Problems Encountered
- Conclusions
- Possible Improvements
- Requirements
Corporación Favorita operates dozens of supermarkets across Ecuador, each selling hundreds of products grouped into families. Demand varies by store, product family, day of the week, season, promotions, and external factors such as oil prices (Ecuador is an oil-dependent economy) and public holidays. Poor demand forecasts lead to stockouts (lost sales) or overstock (waste, especially for perishables).
The question this project addresses is:
Given the sales history and contextual information for each store–product-family combination, can we accurately forecast its daily unit sales?
This is framed as a time series regression problem over panel data (many series observed over time).
The target metric is RMSLE (Root Mean Squared Logarithmic Error), the official competition metric. The target is transformed as log(1 + sales), so RMSLE reduces to RMSE on the log-transformed target. RMSLE penalises relative error, giving equal weight to a 10% miss regardless of a store's sales volume — essential when entities operate at very different scales.
Accurate daily sales forecasts have direct operational value for a retail chain:
- Inventory optimisation. Forecasting demand per store and product family enables better stock allocation, reducing both stockouts and waste — critical for perishables (groceries, produce, beverages).
- Promotion planning. The model quantifies the sales lift of promotions, allowing the business to simulate scenarios ("what if we promote more items in this store-family?") before committing budget.
- Resource and logistics planning. Anticipating demand peaks (weekends, paydays, holidays) supports staffing and distribution decisions.
- Decision support, not replacement. Forecasts are probabilistic estimates that inform, not replace, human planning.
All predictive features (recent sales, promotions, calendar, oil price) are available ahead of time, making daily operational forecasting feasible.
- Source: Kaggle — Store Sales: Time Series Forecasting (Corporación Favorita)
- Period: 2013-01-01 to 2017-08-15 (~4.5 years, 1,684 days)
- Granularity: One row per store × product-family × day
- Size: 3,000,888 rows (54 stores × 33 families × 1,684 days)
- Structure: Five source tables merged into a single panel:
| Table | Key columns |
|---|---|
train |
sales (target), onpromotion (items on promotion) |
stores |
city, state, type, cluster |
transactions |
transactions (daily transactions per store) |
oil |
dcoilwtico (daily WTI oil price) |
holidays_events |
holiday type, locale, transferred flag |
The five source files were merged on store_nbr, date, and family to build a single panel dataset indexed by store-family and date.
holidays_events has a non-trivial structure: a single date can carry multiple holiday records (national + local, transfers, bridges, additional days). It was aggregated per date into boolean flags (is_national, is_regional, is_local, is_holiday, is_transfer, is_event, is_additional, is_bridge) using an any aggregation. Only 312 unique dates carry a holiday flag; the rest are regular days.
dcoilwtico— missing on weekends and public holidays (markets closed). Forward-filled to carry the last known price.transactions— missing for store-days with no recorded transactions; filled with 0.- Holiday flags — missing for non-holiday dates; filled with
False.
31% of all store–family–day combinations have zero sales. These are not missing values — they are genuine zeros (a product family not sold that day in that store). This zero-inflation drove key modelling decisions: an additive decomposition (multiplicative breaks with near-zero values) and the log(1 + sales) transformation (defined at zero).
Sales are highly right-skewed (median 11, max 124,717). The target was transformed with np.log1p (log(1 + sales)) to compress the long tail and align with the RMSLE metric.
Nine business questions guided the EDA:
- Sales concentration by store — a small group of stores (44, 45, 47, 3) concentrates a disproportionate share of total sales; store type and city correlate with volume.
- Zero-sales stores — store 52 stands out with almost entirely zero sales (a store that opened late in the period), confirming the zeros are structural, not random.
- Top product families — GROCERY I and BEVERAGES dominate aggregate sales across nearly all stores.
- Promotion effect — promoting a product measurably increases its sales, with the lift varying by family (from ~10% to ~95%).
- Seasonality — clear weekly pattern (Saturday peak, midweek trough); annual pattern with December peaks.
- Payday effect — sales rise in the second half of the month, consistent with Ecuador's biweekly public-sector paydays.
- Oil price — strong negative correlation on levels (-0.79), but ~0 on first differences → the level correlation is largely spurious (both series trend), though oil still proxies macro conditions.
- Holidays — all holiday types show higher average sales than non-holidays; national holidays +15% (driven by Christmas and Carnival).
- 2016 earthquake — the earthquake week was a sales trough followed by a spike the week after; HOME APPLIANCES dropped then recovered strongly (reconstruction demand).
Beyond the EDA, the aggregate series was analysed with classical time series tools:
- Decomposition (
seasonal_decompose, period=7) — an additive model was chosen over multiplicative, because the 31% structural zeros break the multiplicative form. It cleanly separates a rising trend (200 → 500), stable weekly seasonality, and residuals. - ACF / PACF — the ACF's slow decay signals non-stationarity in levels; peaks at lags 7, 14, 21, 28 confirm weekly seasonality. The PACF shows a dominant lag-1 (~0.75) and direct weekly effects at lags 6–7 → an AR structure.
- ADF stationarity test — with
regression='ct'(constant + trend), p ≈ 1e-6 → the series is trend-stationary, not difference-stationary. Implication: no differencing needed; explicit trend/seasonal features suffice (no ARIMA with d=1).
Before the ML models, a PanelOLS regression with entity fixed effects (one intercept per store-family) was fit as an interpretable econometric benchmark, on log(1 + sales).
Key findings:
- F-test for Poolability (F = 6,389, p < 0.001) → entity fixed effects are justified; each store-family has a distinct baseline that a pooled OLS would confound with the regressors.
- Within R² = 0.16 → the contextual variables alone explain only 16% of within-entity variation; the missing 84% is temporal (lags, seasonality) — motivating the feature engineering that follows.
- Coefficient effects (as
(exp(β) − 1) × 100):onpromotion+1.6% per additional promoted item;dcoilwtico−1.7% per unit (t = −640, highly significant); transferred holidays +50%, additional +36%, events +33%; Saturday +28% vs Friday. - Entity fixed effects rank store-families by baseline demand: GROCERY I and BEVERAGES in large stores at the top; BOOKS, LADIESWEAR, BABY CARE at the floor (non-essential, structurally sparse).
The project follows a structured time series workflow:
- Data loading and merging — five tables into a single panel.
- Missing value treatment — forward-fill oil, fill transactions/holiday flags.
- EDA — nine business questions, visualised.
- Time series analysis — decomposition, ACF/PACF, ADF.
- Panel data modelling — PanelOLS with entity fixed effects.
- Feature engineering — lags, rolling means, calendar features.
- Temporal train/test split — never random for time series.
- Naive baseline — establish a benchmark to beat.
- Model training — LightGBM and XGBoost, tuned with Optuna over time-series cross-validation.
- Evaluation — RMSLE, R², MAE, and bootstrap confidence intervals.
- Explainability — native feature importance and SHAP.
The PanelOLS analysis showed that contextual variables alone are insufficient; temporal features carry the signal.
| Feature group | Features | Rationale |
|---|---|---|
| Lags | sales_log_lag1, lag7, lag14 |
Capture the AR structure (lag-1 ≈ 0.75) and weekly cycle from the ACF/PACF |
| Rolling means | rolling_mean_7, rolling_mean_28 (shifted) |
Smoothed recent level and monthly trend per entity |
| Calendar | dayofweek, month, n_week |
Weekly and annual seasonality |
| Context | onpromotion, dcoilwtico, holiday flags |
Promotions, macro conditions, calendar events |
| Entity | store_nbr, family, city, state, type, cluster |
Identity (label-encoded for tree models) |
Leakage prevention: lags use groupby(entity).shift(k) so each entity only sees its own past; rolling means use .shift(1).rolling(w).mean() so the current day is never included in its own feature. Categorical encoders are fit on train only.
Split: temporal — train < 2017-06-01, test ≥ 2017-06-01 (the most recent ~2.5 months held out). A random split would leak the future into training.
Predicts sales[t] = sales[t-7] (same weekday last week). This is a strong benchmark because of the weekly autocorrelation, and any model must beat it to justify its complexity.
Two global gradient-boosting models (one model learning all 1,782 series jointly, using entity identity and engineered features). The global approach lets the models exploit patterns shared across series and scales far better than fitting 1,782 local models.
- Hyperparameter optimisation: Optuna, 50 trials each, optimising
n_estimators,max_depth(+num_leavesfor LightGBM),learning_rate,subsample,colsample_bytree,reg_alpha,reg_lambda. - Cross-validation: time-series CV with custom date-based folds (expanding window). Because the panel is sorted by entity, a plain
TimeSeriesSplitwould split by entity rather than time — so the folds were built manually from date cutoffs to guarantee "train on past, validate on future" across all entities.
ARIMA was deliberately excluded: it is inherently local (one model per series) and cannot transfer knowledge across the 1,782 series.
Reading the metrics: RMSLE is the primary metric (lower is better) — error on the log scale, i.e. relative error. R² measures variance explained. MAE is the mean absolute error on the log target. Bootstrap CIs (10,000 resamples) show how precisely each metric is estimated.
| Model | CV RMSLE | Test RMSLE | Test R² | Test MAE | Bootstrap 95% CI |
|---|---|---|---|---|---|
| Naive baseline (lag-7) | — | 0.530 | 0.956 | 0.343 | — |
| LightGBM (Optuna) | 0.444 | 0.390 | 0.976 | 0.266 | [0.388, 0.393] |
| XGBoost (Optuna) | 0.441 | 0.381 | 0.977 | 0.258 | [0.378, 0.383] |
- Both models beat the baseline by ~26–28% (RMSLE 0.530 → 0.38), confirming the value of the engineered features over a simple "same day last week" rule.
- XGBoost is the top scorer (0.381 vs 0.390). The non-overlapping bootstrap CIs show the difference is real, but tiny (~0.01 RMSLE; both R² ≈ 0.977).
- No overfitting: train, CV, and test errors are consistent. The CV RMSLE (~0.44) exceeds the test RMSLE (~0.38) because the expanding-window folds validate on 2016 with less training history, not because of a generalisation problem.
- Two independent gradient-boosting implementations converging to almost the same score signals that the engineered features, not the algorithm, set the performance ceiling. Optuna confirmed this empirically — tuning yielded only marginal gains over sensible defaults.
LightGBM is the recommended model for production. It is within ~2.5% of XGBoost's accuracy while training noticeably faster and using less memory — a worthwhile trade for a system that retrains frequently as new sales data arrives. XGBoost is retained as the best-performing model on paper.
Native feature importance (gain) and SHAP (TreeExplainer) were applied to the LightGBM model.
sales_log_rolling_mean_7dominates — the recent 7-day average sales level is by far the strongest predictor. A high recent level pushes predictions up, a low one down (clean, monotonic relationship).sales_log_lag1— second: yesterday's sales strongly predict today's, confirming the AR(1) structure.onpromotion— the key actionable driver. It ranks 3rd in SHAP despite being near-last in native importance: nativegainundervalues it because it is correlated with the dominant rolling-mean feature, while SHAP correctly credits it. More promoted items → higher predicted sales.dayofweek— weekend values push predictions up (weekly seasonality).
Native gain concentrated ~85% of the weight on the rolling mean, masking correlated features. SHAP gives an honest, interaction-aware ranking — the promotion insight only surfaces with SHAP. This illustrates why native importance should not be used alone for feature attribution.
Least relevant: oil price, holiday flags, and store/city identity have minimal impact once the recent sales history is known — the autoregressive features absorb most of the macro and structural signal.
Sales are driven primarily by their own recent momentum (7-day average + previous day). The most important lever the business controls is promotions — the only high-impact feature that is not autoregressive.
The first rolling-mean implementation (rolling(7).mean()) included the current day in its own window — leaking the target. Fixed by shifting first: .shift(1).rolling(7).mean(), so each feature only uses strictly past values. This is invisible in the metrics (it just inflates them), making it a dangerous silent bug.
The panel is sorted by entity (all dates of one store-family, then the next). A plain TimeSeriesSplit splits by row position, so its folds were groups of entities, not time periods — completely defeating temporal validation. Fixed with custom folds built from date cutoffs.
An early version of the holiday EDA accidentally used a single store's data (93% zeros) instead of the full dataset, producing misleading conclusions. Caught and corrected — the fixed analysis reversed the finding (all holiday types are above baseline).
An early feature step computed np.log1p(1 + sales) = log(2 + sales) instead of log(1 + sales). Corrected to np.log1p(sales).
This project demonstrates that daily sales for a large, heterogeneous retail panel (1,782 store–family series) can be forecast accurately using engineered temporal features and gradient boosting.
The best model (XGBoost) reaches RMSLE 0.381 / R² 0.977, a ~28% improvement over a strong naive baseline. LightGBM matches it closely (0.390) and is recommended for production for its speed. The convergence of two independent boosting implementations shows that the engineered features — recent sales momentum above all — set the performance ceiling; Optuna hyperparameter tuning confirmed this, yielding only marginal gains over defaults.
The escalation of models tells a complete story: a naive baseline establishes the bar, PanelOLS provides an interpretable econometric view of which variables move sales and by how much, and gradient boosting delivers the accurate forecasts. Explainability (SHAP) closes the loop, showing that predictions are driven by recent sales momentum, with promotions as the main controllable lever.
Even with a large, zero-inflated, multi-scale real-world panel, a well-engineered feature set plus gradient boosting produces a production-ready forecasting system — accurate, fast to retrain, and interpretable.
- Native categorical handling. Use LightGBM/XGBoost native categorical support (pandas
categorydtype) instead of label encoding, avoiding the false ordinal structure imposed onfamily,city, etc. - Cyclical encoding for
monthandn_week(sin/cos) so December and January are treated as adjacent. - More features — promotion lags, longer rolling windows, days-to/from holiday, explicit payday indicators.
- Sequence models done right — an LSTM or Temporal Fusion Transformer fed with raw daily sequences (not pre-digested lags) could capture temporal patterns the trees miss.
- Segmented models — separate models per product family or store cluster for the most heterogeneous groups.
- Ensemble — blend LightGBM and XGBoost for a small robustness gain.
- Wider Optuna search — more trials and a larger search space, now that the pipeline is validated.
pandas
numpy
matplotlib
seaborn
plotly
statsmodels
linearmodels
scikit-learn
lightgbm
xgboost
optuna
shap
joblib
Install all dependencies:
pip install pandas numpy matplotlib seaborn plotly statsmodels linearmodels scikit-learn lightgbm xgboost optuna shap joblibData source: https://www.kaggle.com/competitions/store-sales-time-series-forecasting