Kaggle Username: chasefarha
Public Score: 0.11910
Private Score: 0.16725

1. Introduction

This project addresses the Forecasting Sticker Sales Kaggle competition, where the objective is to forecast daily sticker sales across a multi-level retail structure involving countries, stores, and products. The dataset contains seven years (2010-2016) of daily sales observations for each unique product–store combination, along with corresponding location information.

The goals of this assignment were to (1) clean and structure the data into a hierarchical time-series (HTS) framework, (2) implement and compare at least four forecasting models, (3) perform time-based cross-validation, (4) reconcile forecasts across the hierarchy, and (5) submit final coherent forecasts to Kaggle. My personal goal with this assignment is to fit a model that accounts for the seasonality, holiday variation, and overall volatility of this data.

A 3-year cross-validation (CV) period was chosen intentionally because the test set also spans 3 years, and using the same horizon for model selection yields more realistic long-horizon accuracy estimates. My approach ultimately involved statistical, regression-based, and machine learning models, along with several reconciliation techniques.

2. Methods

2.1 Data Processing and Feature Engineering

Daily sales were first aggregated using the three-level hierarchy:

  • Country
    • Store
      • Product

The “Holographic Goose” sticker was not sold at the “Discount Stickers” store in Canada or Kenya in the dataset. These missing values were filled with zeros so they could be aggregated into the hierarchical structure properly.

Next, I tested the top-level aggregated data for stationarity. The ACF indicates the raw aggregated data are not stationary, but the ndiffs test yields a result of 1 (not shown). This indicates that adding trend components to models (like the dynamic model and ETS) and using models that first-difference the data (like ARIMA) will be sufficient to achieve stationarity.
ACF Plot of Top-Level Aggregated Sales

ACF Plot of Top-Level Aggregated Sales

Additional temporal features were added, including:

  • Date-based components (day of week, month, year, day of year, weekend, Christmas, New Years, New Years Eve, month of December, Christmas/New Years window(12/22 - 1/3))
  • Quarter variables
  • Fourier terms for complex seasonality

2.2 Hierarchy Structure and Reconciliation

Before modeling, the dataset was converted into an hts-compatible format using the fable/fpp framework.

A three-level summing matrix was constructed using:

  • Countries (6)
  • Stores (3)
  • Products (5)

To ensure coherent forecasts across levels, I use bottom-up, middle-out, top-down, and minT reconciliation methods. Importantly, I attempt the minT method because it improves forecast accuracy and produces coherent, minimum-variance predictions (Hyndman et. al, 2011).

Below are shown the sticker sales aggregated by country in the first figure, and then by product in each panel with different lines for each store in the second figure.
Sticker Sales by Country and Product

Sticker Sales by Country and Product

Sticker Sales by Product and Store

Sticker Sales by Product and Store

2.3 Forecasting Models

I implemented the following models:

1. ETS

Automatic exponential smoothing using state-space ETS models. The error, trend, and seasonality parameters were auto-selected for each product-store-country combination and each aggregation level, creating custom and bespoke models and forecasts for each level. ETS and ARIMA were included because Rambing et. al (2025) show that these models remain strong benchmarks for hierarchical and grouped time series, and typically form the base models used in optimal-combination frameworks such as MinT.

2. ARIMA / SARIMA

I fit two versions of ARIMA models. The first is a stepwise ARIMA selected via AICc minimization, including seasonal components where appropriate. Again, this is done at every aggregation level and for each combination, creating bespoke model parameters for each. The fpp3 package auto-selected weekly (7-day) seasonal models for the vast majority of combinations.

The second model is a SARIMA model. This specification searches across autoregressive orders 0–2, differencing 0–1, and moving‐average orders 0–2 in the nonseasonal component, while allowing seasonal AR, differencing, and MA terms of order 0–1, ensuring there will be some level of seasonality in the parameters for this model.

I then compare the mean and median AICc for each combination for the auto-fitted ARIMA models and the SARIMA model. The auto-fitted model has a lower mean AICc, so this is the ARIMA model I will use going forward.

3. Dynamic Regression (ARIMAX) with Fourier Terms

Next, I fit several dynamic models, incorporating date-based regressors, seasonal and holiday dummies, and Fourier series. I include a trend component into all the linear models here to account for the non-stationarity mentioned above. First, I fit a basic linear model with dummies for day of the week, day of the year, and month.

Then, I fit a second model that adds dummies for Christmas, New Years, New Years Eve, and a full window for the holiday period ranging from December 22 to January 3. The third model keeps all of the previous regressors and adds Fourier series. Specifically, it adds a weekly (7-day) seasonal series with parameter K = 3 and a yearly cyclical series (365.25 days) with parameter K = 5. Adding holiday indicators, date-based regressors, and high-frequency seasonality via Fourier series has been empirically proven to be useful in forecasting retail sales (Taylor & Letham, 2018).

I then fit two additional ARIMAX models, one ARIMA model with all of the dummy variables mentioned above and another with all the dummies plus the two Fourier series. Next, I compare the mean AICc across all combinations to evaluate which model is best. The dynamic linear model with the Fourier series and the dummies has the lowest AICc, so I will go forward with this one.

4. XGBoost

A machine-learning model using lag features, rolling means, holiday windows, Fourier features, and categorical encodings. This aligns with evidence from the M4 competition showing that tree-based ML models can outperform classical approaches when supplied with engineered temporal features (Makridakis et al., 2018).

The features include 1-day, 7-day, and 30-day lagged sales; 7-day and 30-day rolling mean; 7-day rolling standard deviation; year; month; quarter; day of week; day of year; weekend; New Years; Christmas; holiday window (same as above); month of December; a Q4 dummy; and the country, store and product for each combination.

The figure below plots the feature importances, showing that the 7-day lag, 7-day rolling mean, and 30-day rolling mean are the most insightful features for the model. This makes sense given the weekly pattern observed.
XGBoost Feature Importance

XGBoost Feature Importance

2.4 Cross-Validation Strategy

As mentioned in the introduction, a 3-year cross-validation (CV) period was chosen intentionally because the test set also spans 3 years, and using the same horizon for model selection yields more realistic long-horizon accuracy estimates. I trained the models on data from 2010-2013 and cross-validated it on data from 2014-2016 before refitting the best-performing models on the full train set.

I selected the best-performing version for each of the 4 types of models by AICc to determine which versions to move forward with for cross-validation. I used the default auto-select ETS, the auto-select ARIMA, the dynamic linear regression with Fourier series and dummies (dynreg_fourier), and the XGBoost specified above. I compare these four models’ performance using MAPE in the CV, and compare the performance of each reconciliation method for the first three traditional models.

2.5 Model Selection and Refitting

After CV MAPE comparison, the top models (XGBoost and bottom-up reconciled dynamic Fourier regression) were re-fit to the entire training dataset and used to generate test-set predictions.

3. Results

3.1 Cross-Validation Performance

Within each of the three traditional models, the bottom-up reconciliation performed best by MAPE in the CV set. This makes some intuitive sense given that we are aggregating up from the bottom at the product level. Additionally, given that the Kaggle score is based on predictions at the bottom level, it makes sense to stick with the bottom-up approach.

The table below summarizes the average MAPE across all series in the 3-year CV window. XGBoost has the lowest MAPE, with the bottom_up dynamic Fourier model in second.

Model Comparison by MAPE on CV Set
Model Mean MAPE (%)
XGBoost 6.42
Dynamic Regression (Bottom-Up) 18.16
Dynamic Regression Fourier 18.16
ARIMA (Auto) 26.63
ARIMA (Bottom-Up) 27.64
ETS (Bottom-Up) 42.36
ETS (Auto) 43.84

Although XGBoost achieved the lowest CV MAPE, it produced an excessive number of 0 forecasts, especially at the product-store daily granularity. That behavior is hard to correct without major feature redesign or specialized zero-inflated modeling. Because of this instability, I did not select XGBoost as my final model.

Instead, I selected: Dynamic Regression with Fourier Terms (Bottom-Up Reconciled)

This model delivered the second-best CV MAPE, provided stable long-horizon behavior, and produced realistic daily forecasts across the hierarchy. Given the high-frequency seasonal fluctuations in this data, it makes sense that a Fourier regression performs well. This model also forecasts some zeros (which makes sense given that some stickers are not sold at certain stores), but its distribution (shown below) is far more evenly spread across multiple values.
Forecast Distribution for Selected Model

Forecast Distribution for Selected Model

3.2 Forecast Results

The figure below shows the historical data along with my forecasts from the dynamic Fourier model aggregated by country. Although we do not have the observed data here for comparison, the forecasts visually appear to be in line with the levels, trends, and seasonality of the historical data.
Forecasts by Country

Forecasts by Country

3.3 Kaggle Scores

These forecasts earned a public score of 0.11910 and a private score of 0.16725. This can be viewed under my Kaggle username, chasefarha.

The private leaderboard score was higher than the public score, indicating weaker performance in the final evaluation window. One possibility is some alignment of modeling choices with patterns specific to the public leaderboard period, although this risk is limited given the simplicity and high regularization of the final specification, which relies on a linear trend, fixed Fourier seasonality, and a small set of calendar and holiday indicators. An alternative explanation is that the private test window differs structurally from the public window, particularly in how seasonal and calendar effects align over time. Because the model imposes fixed seasonal structure rather than adapting flexibly to recent dynamics, shifts in demand intensity, holiday timing, or product mix could plausibly lead to reduced performance in the private holdout. Overall, the observed gap is consistent with a combination of these factors rather than clear-cut overfitting.
Kaggle Submission Screenshot

Kaggle Submission Screenshot

4. Discussion

This project highlighted several practical forecasting challenges. First, high-frequency retail data contains long periods of zero sales, irregular patterns, and sudden shifts. Classical ETS and ARIMA models, while strong in many contexts, struggled to capture complex multiple seasonalities at the weekly, monthly, and yearly level. Dynamic regression, especially with Fourier terms, handled these patterns more effectively.

XGBoost performed extremely well in cross-validation, consistent with findings from the M4 competition, but at the cost of producing unrealistic runs of zeros. This demonstrates the tradeoff between pure predictive accuracy and forecast plausibility. A retailer would almost certainly prefer slightly higher MAPE in exchange for operationally usable forecasts.

I ultimately used the bottom-up reconciliation method because they performed well across the traditional models, and it makes sense to use the bottom-level forecasts as a base given the nature of this competition. Namely, that the scores are graded on the accuracy of the bottom-level forecasts.

5. Conclusion

This project required building a complete forecasting pipeline—from preprocessing and feature engineering to model estimation, hierarchical reconciliation, cross-validation, and final submission. My final model, a dynamic regression with Fourier terms and bottom-up aggregation, was selected for its balance of accuracy, interpretability, and stable long-horizon behavior. While XGBoost had the best raw CV accuracy, its forecast distribution was not suitable for final deployment.

The final Kaggle scores (public 0.119, private 0.167) demonstrate strong performance for a coherent hierarchical forecasting approach.

If I were to continue improving this solution, I would explore hybrid methods combining ML base forecasts with MinT reconciliation. I would also compute an ensemble forecast that is a simple average of all the methods I used above, and try error-based ensembling across reconciled models. I think another major avenue for improvement would be to implement rolling-origin time-based cross-validation, rather than a single cutoff date. While using a static 3-year period is sufficient for this competition, evaluating the models on multiple CV datasets would certainly improve my ability to pick the best model and fine-tune parameters.

Overall, this project reinforced modern principles of hierarchical forecasting, reconciliation, and long-horizon validation, and provided a realistic end-to-end forecasting experience. The model produced realistic and applicable forecasts that would be of use to retailers, capturing secular trends, seasonality, and volatility.

6. References

Hyndman, R.J., Ahmed, R.A., Athanasopoulos, G., Shang, H.L. (2011). Optimal combination forecasts for hierarchical time series. Computational Statistics & Data Analysis, 55(9), 2579-2589. https://doi.org/10.1016/j.csda.2011.03.006

Makridakis S., Spiliotis E., & Assimakopoulos, V. (2018). Statistical and Machine Learning forecasting methods: Concerns and ways forward. PLOS ONE 13(3): e0194889. https://doi.org/10.1371/journal.pone.0194889

Rambing, D.H., Sugiharto, A., & Kusumaningrum, R. (2025). A Comparative Study of ETS, ARIMA, and Reconciliation Techniques for Hierarchical Product Forecasting in Retail SMEs. 2025 International Conference on Smart Computing, IoT and Machine Learning, SIML. DOI: 10.1109/SIML65326.2025.11081058

Taylor, S.J., & Letham, B. (2018). Forecasting at Scale. The American Statistician, 72(1), 37–45. https://doi.org/10.1080/00031305.2017.1380080