Demand Forecasting: SARIMA vs Prophet

Build a demand forecasting pipeline comparing SARIMA and Prophet models.

Introduction

30 Second Summary

Every retail store faces the same challenge: order too much stock and profits disappear into storage costs, order too little and customers walk away empty-handed. Accurate demand forecasting is how businesses solve this balancing act.

In this project, you will build a complete demand forecasting pipeline that compares two fundamentally different approaches to predicting future sales: SARIMA (classical statistics) and Prophet (Meta's additive decomposition model).

What You'll Build

A Jupyter Notebook pipeline that takes 5 years of daily retail sales data and produces 90-day forecasts using SARIMA and Prophet, with proper temporal splitting and quantified performance metrics.

By the end of this project, you'll have:

  • A time series decomposition showing trend, seasonality, and residuals extracted from real retail sales data.
  • Two trained forecasting models (SARIMA and Prophet) evaluated on the same holdout period using a leakage-free temporal split.
  • A visual comparison chart with quantified MAE/RMSE metrics showing which model captures the data patterns better.
  • Secret Mission: Add an XGBoost model with hand-crafted lag and rolling features to see how ML approaches compare when you must explicitly engineer what statistical models capture automatically.

Are there any prerequisites?

This project assumes you are comfortable writing Python and have used pandas and numpy before. Basic statistics knowledge (what a mean and standard deviation are) is helpful but not required.

Before We Start

You're about to build a demand forecasting pipeline that compares two fundamentally different approaches to predicting the future: SARIMA (classical statistics) and Prophet (Meta's additive decomposition model).

Before diving into the data science, you need the right tools in place. In this step, you'll install the forecasting libraries, download a 5-year retail sales dataset, and launch the Jupyter Notebook environment where you'll build the entire pipeline.

In this step, get ready to:
  • Install the forecasting and visualization libraries.
  • Download the Kaggle retail sales dataset.
  • Launch Jupyter Notebook and create your project notebook.
Install the required packages

This project uses statsmodels for classical SARIMA modeling and prophet for Meta's additive forecasting approach. You'll also need matplotlib for visualizations and jupyter for the notebook environment.

  • Open your terminal. Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal (macOS) or Command Prompt (Windows) and press Enter to open it.
  • Install all required packages by running this command:
pip install statsmodels==0.14.6 prophet==1.3.0 matplotlib jupyter

What does this command do?

This installs four packages with specific version pins to ensure reproducibility:

  • statsmodels==0.14.6 provides the SARIMAX model and time series decomposition tools.
  • prophet==1.3.0 is Meta's forecasting library with automatic seasonality detection.
  • matplotlib handles all chart rendering.
  • jupyter provides the interactive notebook interface.

The installation may take a few minutes. Prophet downloads its own copy of CmdStan (a statistical computing backend) during installation, so expect some extra output.

  • Verify the installation succeeded by running:
python -c "import statsmodels; import prophet; print('statsmodels:', statsmodels.__version__); print('prophet:', prophet.__version__)"

✔️ I see both version numbers

You should see output like:

statsmodels: 0.14.6
prophet: 1.3.0

Both libraries are ready. Move on to the next substep.

ⓧ I see an error during installation

Prophet installation can be tricky, especially on Windows. Try these fixes:

  • On Windows, Prophet requires a C++ compiler. If you see build errors, try installing via conda instead (see below).
  • Make sure your pip is up to date by running python -m pip install --upgrade pip first.
  • If you see path-related errors on Windows, make sure you're installing exactly prophet==1.3.0 and not an older version.

As a fallback on Windows, try installing Prophet through conda-forge:

conda install -c conda-forge prophet

After resolving the issue, re-run the verification command above to confirm both packages load correctly.

Download the dataset

You'll use the Kaggle Store Item Demand Forecasting Challenge dataset. It contains 913,000 rows of daily sales data across 10 stores and 50 items, spanning 2013 to 2017.

What's in this dataset?

The file has four columns: date, store, item, and sales. Each row records how many units of a specific item sold at a specific store on a given day.

You'll aggregate all stores and items into a single daily total to create one clean univariate time series.

  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder and a data subdirectory by running these commands:
mkdir demand-forecasting
cd demand-forecasting
mkdir data
  • Move the downloaded train.csv file into the data/ folder inside demand-forecasting. You can drag and drop it using Finder (macOS) or File Explorer (Windows).
  • Confirm the file is in place by running:
ls data/

You should see train.csv listed in the output.

Don't see train.csv?

  • Make sure you downloaded train.csv specifically, not the full zip archive. If you downloaded a zip, unzip it first.
  • Double-check the file landed inside demand-forecasting/data/ on your Desktop, not in your Downloads folder.
  • On Windows, use dir data\ instead of ls data/ if you're using Command Prompt.

Help me find my downloaded dataset

Launch Jupyter Notebook

Jupyter Notebook gives you an interactive environment where you can run code cell by cell and see visualizations inline. This is ideal for exploratory data science work.

  • Make sure you're inside the demand-forecasting directory, then launch Jupyter by running:
jupyter notebook

A browser window opens showing the Jupyter file browser. You should see your data/ folder listed.

  • Click New in the top-right corner of the Jupyter interface.
  • Select Python 3 from the dropdown to create a new notebook.
  • Click Untitled at the top of the notebook to rename it.
  • Type demand_forecasting and press Enter to save the name.

Why a notebook instead of a script?

Time series forecasting is exploratory. You'll want to see decomposition plots, check model diagnostics, and compare forecasts visually.

Notebooks let you run code in small chunks and keep outputs visible. This is much faster for iteration than rerunning an entire script each time.

Browser didn't open?

  • Check your terminal output for a URL starting with http://localhost:8888. Copy and paste it into your browser manually.
  • If you see a command not found error, the jupyter package may not be on your PATH. Try running python -m jupyter notebook instead.

Help me troubleshoot Jupyter

Your environment is ready and your notebook is open. Next up, you'll load the raw sales data into a pandas DataFrame and explore the time series patterns hiding inside 5 years of retail transactions.

Load, Aggregate, and Explore the Data

Your environment is ready and the dataset is waiting. Before you can forecast anything, you need to understand what the data looks like.

The raw CSV has 913,000 rows covering 10 stores and 50 items. For a univariate forecasting pipeline, you need a single time series of total daily sales. In this step, you will aggregate the data, plot it, and decompose it into its underlying components.

In this step, get ready to:
  • Import all libraries and load the dataset into a single daily sales time series.
  • Visualize the raw time series to spot trends and seasonality.
  • Decompose the series into trend, seasonal, and residual components.
Import libraries and load the dataset

You will import every library the project needs up front. This includes data manipulation, visualization, statistical modeling, and evaluation tools.

  • In your open notebook, click into the first empty cell and paste the following imports:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from prophet import Prophet
from sklearn.metrics import mean_absolute_error, mean_squared_error

What does this code do?

Each import serves a specific role in the pipeline:

  • pandas and numpy handle data loading and numerical operations.
  • matplotlib produces all the charts and visualizations.
  • seasonal_decompose and adfuller come from statsmodels for decomposition and stationarity testing.
  • SARIMAX is the classical statistical forecasting model.
  • Prophet is Meta's additive decomposition forecasting library.
  • mean_absolute_error and mean_squared_error from scikit-learn evaluate forecast accuracy.
  • Run the cell by pressing Shift+Enter.
  • The cell should complete without errors. If it runs successfully, you will see no output and the cursor moves to the next cell.

Seeing a ModuleNotFoundError?

  • Make sure you installed the packages in the same Python environment that Jupyter is using. Run !pip install statsmodels prophet inside a notebook cell to install directly into the active kernel.
  • If Prophet fails to install on Windows, try restarting your terminal and running the install command again.

Why is my import failing?

Now you will load the CSV and aggregate it into a single daily time series.

  • In the next empty cell, paste this code to load and aggregate the data:
# --- Step 2: Load and Explore ---
df = pd.read_csv("data/train.csv", parse_dates=["date"])
daily_sales = df.groupby("date")["sales"].sum()
daily_sales.index = pd.DatetimeIndex(daily_sales.index, freq="D")

What does this code do?

  • pd.read_csv("data/train.csv", parse_dates=["date"]) loads the CSV and converts the date column into proper datetime objects.
  • df.groupby("date")["sales"].sum() sums sales across all 10 stores and 50 items for each day, producing one number per date.
  • pd.DatetimeIndex(daily_sales.index, freq="D") explicitly sets the index frequency to daily. statsmodels requires a known frequency to fit time series models.
  • Run the cell with Shift+Enter.
  • Verify by adding daily_sales.shape in the next cell and running it. You should see (1826,) confirming 1,826 days of data from 2013 through 2017.

Seeing a FileNotFoundError?

  • Confirm that train.csv is inside a data/ folder in the same directory as your notebook.
  • Check that you launched jupyter notebook from the project directory, not from a parent folder.

Help me fix my file path for loading the CSV.

Visualize the raw time series

Plotting the full series helps you visually identify patterns before applying any mathematical methods. You are looking for trend (does it go up or down over time?) and seasonality (are there repeating cycles?).

  • Plot the total daily sales over all 5 years by pasting this into the next cell:
# Plot raw time series
plt.figure(figsize=(14, 5))
plt.plot(daily_sales, linewidth=0.7)
plt.title("Total Daily Sales (All Stores, All Items)")
plt.xlabel("Date")
plt.ylabel("Total Sales")
plt.tight_layout()
plt.show()

What does this code do?

  • plt.figure(figsize=(14, 5)) creates a wide chart so you can see all 5 years without cramping.
  • linewidth=0.7 thins the line so individual daily fluctuations stay visible within the overall trend.
  • plt.tight_layout() prevents axis labels from being cut off at the edges.
  • Run the cell with Shift+Enter.
  • You should see a line chart showing daily sales from 2013 to 2017. Notice the clear upward trend and the repeating wave pattern that cycles every year.

Plot not showing up?

  • If you see a blank output, try adding %matplotlib inline at the top of your notebook (in the imports cell) and re-run both cells.
  • Make sure the previous cell ran without errors. If your data did not load correctly, the plot cell will fail.

Why is my matplotlib plot not displaying in Jupyter?

Decompose the series into components

Time series decomposition separates your data into three distinct signals: the long-term trend, repeating seasonal patterns, and leftover noise (residuals). This confirms whether the patterns you spotted visually are real and measurable.

  • Decompose the series and display all four components by pasting this into the next cell:
# Decompose
decomposition = seasonal_decompose(daily_sales, model="additive", period=365)
decomposition.plot()
plt.tight_layout()
plt.show()

What does this code do?

  • seasonal_decompose(daily_sales, model="additive", period=365) breaks the series into trend + seasonal + residual. The additive model assumes these components add up to the observed value.
  • period=365 tells the function to look for a pattern that repeats every 365 days (yearly seasonality).
  • decomposition.plot() produces a four-panel chart: observed data on top, then trend, seasonal cycle, and residuals at the bottom.

Why additive instead of multiplicative?

An additive model works when the seasonal fluctuations stay roughly the same size regardless of the trend level. In this dataset, the daily ups and downs do not grow proportionally with the trend. If the seasonal swings grew larger as sales increased, a multiplicative model would be the better choice.

  • Run the cell with Shift+Enter.
  • You should see a four-panel plot. The trend panel shows a clear upward line. The seasonal panel shows a repeating yearly wave. The residuals panel shows random noise with no obvious pattern.

Seeing a ValueError about period?

  • This error means your daily_sales index does not have its frequency set. Go back to the data-loading cell and confirm the line daily_sales.index = pd.DatetimeIndex(daily_sales.index, freq="D") is present and ran without error.
  • If you see a warning about NaN values at the edges of the decomposition plot, that is normal. The trend uses a moving average that cannot compute values for the first and last ~182 days.

Help me fix my seasonal_decompose error.

✔️ Awesome, I've got everything!

Great. Double check that your notebook is saved before continuing.

ⓧ I'd like to double check the full code

Here is the complete code for all cells so far. Compare each cell in your notebook against this:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from prophet import Prophet
from sklearn.metrics import mean_absolute_error, mean_squared_error

# --- Step 2: Load and Explore ---
df = pd.read_csv("data/train.csv", parse_dates=["date"])
daily_sales = df.groupby("date")["sales"].sum()
daily_sales.index = pd.DatetimeIndex(daily_sales.index, freq="D")

# Plot raw time series
plt.figure(figsize=(14, 5))
plt.plot(daily_sales, linewidth=0.7)
plt.title("Total Daily Sales (All Stores, All Items)")
plt.xlabel("Date")
plt.ylabel("Total Sales")
plt.tight_layout()
plt.show()

# Decompose
decomposition = seasonal_decompose(daily_sales, model="additive", period=365)
decomposition.plot()
plt.tight_layout()
plt.show()

You now have a clean daily time series and visual proof of the trend and seasonality baked into it. Next up, you will test whether the series is stationary and split it into training and test sets for fair model evaluation.

Test Stationarity and Create the Temporal Split

You have a clear picture of your time series from the decomposition. Before you can model it, you need to know whether the data is stationary. A stationary series has a constant mean and variance over time. Most forecasting models (including SARIMA) assume stationarity or apply differencing to achieve it.

In this step, you will test for stationarity, create a proper time-aware train/test split, and define the evaluation metrics you will use to compare models later.

In this step, get ready to:
  • Run the Augmented Dickey-Fuller test to check stationarity.
  • Split the series into training and test sets using a temporal boundary.
  • Define MAE and RMSE helper functions for model evaluation.
Run the Augmented Dickey-Fuller test

The Augmented Dickey-Fuller (ADF) test is a statistical hypothesis test for stationarity. The null hypothesis is that the series has a unit root (is non-stationary). A low p-value (below 0.05) means you can reject that null hypothesis and conclude the series is stationary.

  • Add a new cell in your notebook and run the ADF test on your daily_sales series:
# --- Step 3: Stationarity and Split ---
adf_result = adfuller(daily_sales.dropna())
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
if adf_result[1] > 0.05:
    print("Series is NON-STATIONARY (p > 0.05). Differencing needed.")
else:
    print("Series is STATIONARY (p <= 0.05).")

What does this code do?

  • adfuller(daily_sales.dropna()) runs the ADF test on the series. The .dropna() ensures no missing values interfere with the calculation.
  • The function returns a tuple where index 0 is the test statistic and index 1 is the p-value.
  • The if block interprets the result. If the p-value exceeds 0.05, the series is non-stationary and SARIMA will need differencing (the d and D parameters in its order).
  • Run the cell. You should see the ADF statistic, p-value, and an interpretation message printed below the cell.

What should I see?

Depending on your data, the output will tell you whether the series is stationary or non-stationary. For retail sales data with an upward trend, you will likely see a message saying differencing is needed. This confirms why SARIMA uses the differencing parameter.

Seeing a NameError for adfuller?

Make sure you ran the imports cell from Step 2 first. The import from statsmodels.tsa.stattools import adfuller must have executed successfully before this cell.

Help me fix this NameError

Create the temporal train/test split

In time series forecasting, you cannot randomly shuffle data into train and test sets. Random splits leak future information into training, which inflates accuracy and produces useless models in production.

Instead, you use a temporal split: everything before a cutoff date is training data, and everything after is the test set. This simulates a real forecasting scenario where you only have the past to predict the future.

  • Add a new cell and create the temporal split, using the last 92 days (October through December 2017) as your test set:
# Temporal split: last 92 days as test
split_date = "2017-09-30"
train = daily_sales[:split_date]
test = daily_sales[split_date:].iloc[1:]  # 2017-10-01 onward

print(f"Training: {train.index[0].date()} to {train.index[-1].date()} ({len(train)} days)")
print(f"Test:     {test.index[0].date()} to {test.index[-1].date()} ({len(test)} days)")

What does this code do?

  • split_date = "2017-09-30" defines the cutoff. Everything up to and including this date becomes training data.
  • daily_sales[:split_date] slices the series from the beginning through September 30, 2017.
  • daily_sales[split_date:].iloc[1:] takes everything from the split date onward, then drops the first row (September 30 itself) so the test set starts cleanly on October 1.
  • The print statements confirm the exact date ranges and number of days in each set.
  • Run the cell. You should see training dates from 2013-01-01 to 2017-09-30 and test dates from 2017-10-01 to 2017-12-31 (92 days).

Why not a random split?

If you randomly sampled 92 days from anywhere in the 5-year range, your model would "see" December 2017 data during training and then be asked to predict other December 2017 days. That is data leakage. The model appears accurate but fails completely on genuinely unseen future data.

A temporal split guarantees the model never sees data from the test period during training. This is the only valid approach for time series evaluation.

  • Visualize the split by adding this code in a new cell:
# Plot split
plt.figure(figsize=(14, 5))
plt.plot(train, label="Train", linewidth=0.7)
plt.plot(test, label="Test", linewidth=0.7, color="orange")
plt.axvline(pd.Timestamp(split_date), color="red", linestyle="--", label="Split")
plt.legend()
plt.title("Temporal Train/Test Split")
plt.tight_layout()
plt.show()

What does this plot show?

The chart draws the training data in blue and the test data in orange, with a red dashed vertical line marking the split point. This makes it visually obvious that the test period is a continuous block at the end of the timeline, not randomly scattered throughout.

  • Run the cell. You should see a time series plot with blue training data on the left, a red dashed split line, and orange test data on the right covering the final 92 days.

Test set showing more or fewer than 92 days?

Double check that split_date is set to exactly "2017-09-30" and that you included .iloc[1:] when creating the test set. Without .iloc[1:], September 30 would be included in both sets.

Help me debug my temporal split

Define evaluation metrics

To compare SARIMA and Prophet fairly, you need consistent metrics calculated the same way for both models. MAE (Mean Absolute Error) gives the average magnitude of errors in the same units as sales. RMSE (Root Mean Squared Error) penalizes large errors more heavily.

  • Add a new cell and define the helper function that calculates both metrics:
def calc_metrics(actual, predicted):
    mae = mean_absolute_error(actual, predicted)
    rmse = np.sqrt(mean_squared_error(actual, predicted))
    return mae, rmse

What does this function do?

  • mean_absolute_error from scikit-learn calculates the average absolute difference between actual and predicted values.
  • mean_squared_error calculates the average squared difference. Wrapping it in np.sqrt() gives RMSE, which is in the same units as the original data.
  • Returning both as a tuple lets you unpack them later with a single call: mae, rmse = calc_metrics(test, predictions).
  • Run the cell. You will not see any output yet. This function will be called in the next steps when you generate forecasts from each model.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is all the code your notebook should contain through the end of this step:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from prophet import Prophet
from sklearn.metrics import mean_absolute_error, mean_squared_error

# --- Step 2: Load and Explore ---
df = pd.read_csv("data/train.csv", parse_dates=["date"])
daily_sales = df.groupby("date")["sales"].sum()
daily_sales.index = pd.DatetimeIndex(daily_sales.index, freq="D")

# Plot raw time series
plt.figure(figsize=(14, 5))
plt.plot(daily_sales, linewidth=0.7)
plt.title("Total Daily Sales (All Stores, All Items)")
plt.xlabel("Date")
plt.ylabel("Total Sales")
plt.tight_layout()
plt.show()

# Decompose
decomposition = seasonal_decompose(daily_sales, model="additive", period=365)
decomposition.plot()
plt.tight_layout()
plt.show()

# --- Step 3: Stationarity and Split ---
adf_result = adfuller(daily_sales.dropna())
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
if adf_result[1] > 0.05:
    print("Series is NON-STATIONARY (p > 0.05). Differencing needed.")
else:
    print("Series is STATIONARY (p <= 0.05).")

# Temporal split: last 92 days as test
split_date = "2017-09-30"
train = daily_sales[:split_date]
test = daily_sales[split_date:].iloc[1:]  # 2017-10-01 onward

print(f"Training: {train.index[0].date()} to {train.index[-1].date()} ({len(train)} days)")
print(f"Test:     {test.index[0].date()} to {test.index[-1].date()} ({len(test)} days)")

# Plot split
plt.figure(figsize=(14, 5))
plt.plot(train, label="Train", linewidth=0.7)
plt.plot(test, label="Test", linewidth=0.7, color="orange")
plt.axvline(pd.Timestamp(split_date), color="red", linestyle="--", label="Split")
plt.legend()
plt.title("Temporal Train/Test Split")
plt.tight_layout()
plt.show()


def calc_metrics(actual, predicted):
    mae = mean_absolute_error(actual, predicted)
    rmse = np.sqrt(mean_squared_error(actual, predicted))
    return mae, rmse

Your stationarity test is done, your temporal split is clean, and your evaluation metrics are ready. Next up, you will fit a SARIMA model to the training data and generate your first 92-day forecast.

Fit SARIMA and Generate Forecasts

In the previous step, you confirmed the stationarity of your data and created a leakage-free temporal split. Now it's time to fit your first forecasting model and see how well it predicts those 92 holdout days.

SARIMA (Seasonal ARIMA) extends the classical ARIMA framework by adding seasonal differencing and seasonal AR/MA terms. You'll configure it to capture the weekly sales cycle, then generate point forecasts with confidence intervals for the test period.

In this step, get ready to:
  • Fit a SARIMA model with weekly seasonality on the training data.
  • Run residual diagnostics to check model assumptions.
  • Generate 92-day forecasts with 95% confidence intervals and plot them against actuals.
Fit a SARIMA model on the training data

The SARIMAX class in statsmodels implements SARIMA using a state-space representation and the Kalman filter. You'll specify both a non-seasonal order and a seasonal order tuned for the 7-day weekly cycle in the retail data.

  • In your notebook, create a new cell and fit the SARIMA model by running this code:
# --- Step 4: SARIMA ---
model_sarima = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 7),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
results_sarima = model_sarima.fit(disp=False)
print(results_sarima.summary())

What does this code do?

The order=(1, 1, 1) sets the non-seasonal ARIMA parameters: 1 autoregressive term, 1 difference, and 1 moving average term. This handles the overall trend in the data.

The seasonal_order=(1, 1, 1, 7) adds the same structure at a weekly lag. The 7 tells the model that the seasonal cycle repeats every 7 days. Seasonal differencing removes the weekly pattern so the model can learn from the deseasonalized residuals.

Setting enforce_stationarity=False and enforce_invertibility=False prevents the optimizer from failing due to numerical boundary issues. The disp=False suppresses convergence output so only the final summary prints.

  • Run the cell. This may take 30-60 seconds on ~1,700 daily observations.
  • You should see a model summary table with coefficient estimates, standard errors, and information criteria (AIC, BIC).

Cell running for more than 2 minutes?

SARIMAX fitting uses the Kalman filter, which scales with both dataset length and the seasonal period. With s=7 and ~1,700 rows, fitting should complete within a minute or two.

If it hangs, check that your train variable has a proper DatetimeIndex with freq="D". A missing frequency attribute forces the model to infer structure, which slows things down.

Help me debug my SARIMA model fitting

Run model diagnostics

A good model should have residuals that look like white noise. The diagnostics plot checks this assumption across four different views.

  • In a new cell, generate the diagnostics plot by running:
# Diagnostics
results_sarima.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

What do the four diagnostic panels show?

  • Standardized residuals: should fluctuate randomly around zero with no visible patterns or trends.
  • Histogram plus KDE: the residual distribution should roughly follow the orange normal curve.
  • Q-Q plot: points should follow the diagonal red line. Deviations at the tails indicate heavy-tailed residuals.
  • Correlogram (ACF): bars should stay within the blue shaded confidence band. Significant spikes indicate the model missed some autocorrelation structure.
  • Run the cell and examine the four-panel output.
  • You should see residuals that are approximately normally distributed with no large autocorrelation spikes.

Seeing large spikes in the correlogram?

Some residual autocorrelation is expected with a simple (1,1,1)(1,1,1,7) specification. As long as most bars stay within the shaded band and there is no obvious repeating pattern, the model is adequate for forecasting purposes.

Help me interpret my SARIMA diagnostics

Generate forecasts for the test period

With the model fitted, you can now produce predictions for the 92-day test window. The get_forecast method returns both point estimates and confidence intervals so you can see how uncertain the model is about future values.

  • In a new cell, generate the SARIMA forecast and plot it against actuals by running:
# Forecast
forecast_sarima = results_sarima.get_forecast(steps=len(test))
sarima_pred = forecast_sarima.predicted_mean
sarima_ci = forecast_sarima.summary_frame(alpha=0.05)

# Plot SARIMA forecast vs actuals
plt.figure(figsize=(14, 6))
plt.plot(test.index, test.values, label="Actual", color="black", linewidth=1.5)
plt.plot(test.index, sarima_pred.values, label="SARIMA Forecast", linestyle="--")
plt.fill_between(
    sarima_ci.index,
    sarima_ci["mean_ci_lower"],
    sarima_ci["mean_ci_upper"],
    alpha=0.1,
    label="SARIMA 95% CI",
)
plt.legend()
plt.title("SARIMA Forecast vs Actual (Test Period)")
plt.xlabel("Date")
plt.ylabel("Total Daily Sales")
plt.tight_layout()
plt.show()

What does this code do?

The get_forecast(steps=len(test)) generates predictions for 92 future time steps beyond the training data. It returns a forecast results object that holds both the predicted values and uncertainty estimates.

The predicted_mean attribute gives you the point forecast for each day. The summary_frame(alpha=0.05) produces a DataFrame with columns mean_ci_lower and mean_ci_upper representing the 95% confidence interval boundaries.

The fill_between shading makes it easy to see where the model is confident versus uncertain. Wider bands mean more uncertainty further into the future.

  • Run the cell.
  • You should see a line chart with the actual test values in black and the SARIMA forecast as a dashed line, surrounded by a light shaded confidence band.

Forecast line looks flat or unreasonable?

If your forecast is a flat line or drifts far from the actuals, check that train ends at 2017-09-30 and has the correct freq="D" attribute on its index. A misaligned or missing frequency can cause the model to produce degenerate forecasts.

Help me fix my SARIMA forecast output

✔️ Awesome, I've got everything!

Great. Double check you've saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is all the code you should have in your notebook through Step 4:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.statespace.sarimax import SARIMAX
from prophet import Prophet
from sklearn.metrics import mean_absolute_error, mean_squared_error

# --- Step 2: Load and Explore ---
df = pd.read_csv("data/train.csv", parse_dates=["date"])
daily_sales = df.groupby("date")["sales"].sum()
daily_sales.index = pd.DatetimeIndex(daily_sales.index, freq="D")

# Plot raw time series
plt.figure(figsize=(14, 5))
plt.plot(daily_sales, linewidth=0.7)
plt.title("Total Daily Sales (All Stores, All Items)")
plt.xlabel("Date")
plt.ylabel("Total Sales")
plt.tight_layout()
plt.show()

# Decompose
decomposition = seasonal_decompose(daily_sales, model="additive", period=365)
decomposition.plot()
plt.tight_layout()
plt.show()

# --- Step 3: Stationarity and Split ---
adf_result = adfuller(daily_sales.dropna())
print(f"ADF Statistic: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
if adf_result[1] > 0.05:
    print("Series is NON-STATIONARY (p > 0.05). Differencing needed.")
else:
    print("Series is STATIONARY (p <= 0.05).")

# Temporal split: last 92 days as test
split_date = "2017-09-30"
train = daily_sales[:split_date]
test = daily_sales[split_date:].iloc[1:]  # 2017-10-01 onward

print(f"Training: {train.index[0].date()} to {train.index[-1].date()} ({len(train)} days)")
print(f"Test:     {test.index[0].date()} to {test.index[-1].date()} ({len(test)} days)")

# Plot split
plt.figure(figsize=(14, 5))
plt.plot(train, label="Train", linewidth=0.7)
plt.plot(test, label="Test", linewidth=0.7, color="orange")
plt.axvline(pd.Timestamp(split_date), color="red", linestyle="--", label="Split")
plt.legend()
plt.title("Temporal Train/Test Split")
plt.tight_layout()
plt.show()


def calc_metrics(actual, predicted):
    mae = mean_absolute_error(actual, predicted)
    rmse = np.sqrt(mean_squared_error(actual, predicted))
    return mae, rmse


# --- Step 4: SARIMA ---
model_sarima = SARIMAX(
    train,
    order=(1, 1, 1),
    seasonal_order=(1, 1, 1, 7),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
results_sarima = model_sarima.fit(disp=False)
print(results_sarima.summary())

# Diagnostics
results_sarima.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()

# Forecast
forecast_sarima = results_sarima.get_forecast(steps=len(test))
sarima_pred = forecast_sarima.predicted_mean
sarima_ci = forecast_sarima.summary_frame(alpha=0.05)

# Plot SARIMA forecast vs actuals
plt.figure(figsize=(14, 6))
plt.plot(test.index, test.values, label="Actual", color="black", linewidth=1.5)
plt.plot(test.index, sarima_pred.values, label="SARIMA Forecast", linestyle="--")
plt.fill_between(
    sarima_ci.index,
    sarima_ci["mean_ci_lower"],
    sarima_ci["mean_ci_upper"],
    alpha=0.1,
    label="SARIMA 95% CI",
)
plt.legend()
plt.title("SARIMA Forecast vs Actual (Test Period)")
plt.xlabel("Date")
plt.ylabel("Total Daily Sales")
plt.tight_layout()
plt.show()

Your first forecasting model is trained and producing predictions with confidence intervals. Next up, you'll fit a Prophet model and compare both approaches on the same holdout window to see which captures the data patterns better.

Fit Prophet and Compare Both Models

Your SARIMA model is trained and producing forecasts with confidence intervals. Now you'll fit a fundamentally different model on the same data to see how a different philosophy tackles the same forecasting problem.

Prophet is Meta's additive decomposition model. Instead of the state-space approach SARIMA uses, Prophet decomposes the series into trend + seasonality + holidays using Fourier series. You'll fit it on the same training window, generate forecasts for the same 92-day holdout, and compare both models head-to-head with MAE and RMSE.

In this step, get ready to:
  • Prepare training data in Prophet's required format.
  • Fit a Prophet model and generate 92-day forecasts.
  • Compare SARIMA and Prophet with metrics and a side-by-side chart.
Prepare data for Prophet

Prophet expects a DataFrame with exactly two columns: ds (datestamp) and y (the numeric value to forecast). This is different from the pandas Series with a DatetimeIndex that statsmodels uses.

  • In a new cell in your notebook, create the Prophet-formatted training DataFrame by running this code:
# Prepare training data in Prophet's required format
train_prophet = pd.DataFrame({"ds": train.index, "y": train.values})

What does this code do?

This creates a two-column DataFrame from your existing train Series. The DatetimeIndex becomes the ds column, and the sales values become the y column.

You won't see output from this cell. The next chunk uses this DataFrame to fit the model.

Fit Prophet and generate forecasts

Prophet handles weekly and yearly seasonality automatically via Fourier terms. You'll enable both and disable daily seasonality since the data is already at a daily granularity.

  • In the next cell, fit the Prophet model and generate predictions by running:
# Initialize Prophet with weekly and yearly seasonality
m = Prophet(weekly_seasonality=True, yearly_seasonality=True, daily_seasonality=False)
m.fit(train_prophet)

# Create future dates covering the test period
future = m.make_future_dataframe(periods=len(test))
prophet_forecast = m.predict(future)

# Extract test period predictions
prophet_pred = prophet_forecast.set_index("ds").loc[test.index, "yhat"]

What does this code do?

Prophet is initialized with weekly and yearly seasonality enabled. Setting daily_seasonality=False is correct because your data points are already daily totals.

make_future_dataframe(periods=len(test)) creates a date range extending from the training data through the test period. The predict call returns a DataFrame with yhat (point forecast), yhat_lower, and yhat_upper (confidence bounds).

The last line extracts only the predictions that overlap with your test period so you can compare them directly against actuals.

  • Run the cell. Fitting should take 1-5 seconds.
  • Verify the output by checking that prophet_forecast contains columns ds, yhat, yhat_lower, and yhat_upper.

Prophet fit taking too long or erroring?

If you see a message about CmdStan compilation or a missing compiler, Prophet's backend may need to compile on first run. Wait a minute for it to complete.

If fitting fails entirely, ensure you installed Prophet with pip install prophet==1.3.0 and restart your notebook kernel before retrying.

Help me debug my Prophet installation.

Compare both models on the holdout period

Both models have now generated forecasts for the same 92-day test window. A fair comparison requires evaluating them on identical data with the same metrics.

  • Calculate MAE and RMSE for both models and print a comparison table by running this code in a new cell:
# Calculate metrics for both models
sarima_mae, sarima_rmse = calc_metrics(test.values, sarima_pred.values)
prophet_mae, prophet_rmse = calc_metrics(test.values, prophet_pred.values)

# Print comparison table
print("\n--- Model Comparison ---")
print(f"{'Model':<10} {'MAE':>10} {'RMSE':>10}")
print(f"{'SARIMA':<10} {sarima_mae:>10.2f} {sarima_rmse:>10.2f}")
print(f"{'Prophet':<10} {prophet_mae:>10.2f} {prophet_rmse:>10.2f}")

What does this code do?

calc_metrics is the helper function you defined in Step 3. It returns MAE (average absolute error) and RMSE (penalizes large errors more heavily).

Both models are evaluated against the same test.values array. The formatted print statements produce a neat side-by-side comparison table so you can see which model performed better at a glance.

  • Run the cell. You should see a printed table with MAE and RMSE values for both SARIMA and Prophet.

Now you'll create a visual comparison showing both forecasts against the actual test values.

  • Create the side-by-side forecast comparison plot by running this code in a new cell:
# Comparison plot
plt.figure(figsize=(14, 6))
plt.plot(test.index, test.values, label="Actual", color="black", linewidth=1.5)
plt.plot(test.index, sarima_pred.values, label=f"SARIMA (MAE={sarima_mae:.0f})", linestyle="--")
plt.plot(test.index, prophet_pred.values, label=f"Prophet (MAE={prophet_mae:.0f})", linestyle="--")
plt.fill_between(
    sarima_ci.index,
    sarima_ci["mean_ci_lower"],
    sarima_ci["mean_ci_upper"],
    alpha=0.1,
    label="SARIMA 95% CI",
)
plt.legend()
plt.title("Forecast Comparison: SARIMA vs Prophet (90-Day Holdout)")
plt.xlabel("Date")
plt.ylabel("Total Daily Sales")
plt.tight_layout()
plt.show()

What does this chart show?

The black line is the actual sales data from the holdout period. The two dashed lines are each model's predictions, labeled with their MAE scores in the legend.

The shaded region around the SARIMA forecast represents its 95% confidence interval. Where a dashed line stays close to the black line, that model is forecasting well. Where it diverges, the model is missing patterns in the data.

  • Run the cell. You should see a plot with the actual values in black, SARIMA as a dashed line, Prophet as a dashed line, and a light shaded confidence band around the SARIMA forecast.

Plot looks empty or lines are flat?

Make sure you ran the metrics cell first. The plot uses sarima_mae and prophet_mae in the legend labels, so those variables must exist.

If prophet_pred is empty, check that test.index and the ds column in prophet_forecast use the same date format. Both should be datetime objects.

Help me debug my forecast comparison plot.

✔️ Awesome, I've got everything!

Great. Double check you've saved your notebook with all cells executed.

ⓧ I'd like to double check the full code

Here is the complete code added in this step:

# --- Step 5: Prophet ---
train_prophet = pd.DataFrame({"ds": train.index, "y": train.values})

m = Prophet(weekly_seasonality=True, yearly_seasonality=True, daily_seasonality=False)
m.fit(train_prophet)

future = m.make_future_dataframe(periods=len(test))
prophet_forecast = m.predict(future)

# Extract test period predictions
prophet_pred = prophet_forecast.set_index("ds").loc[test.index, "yhat"]

# --- Comparison ---
sarima_mae, sarima_rmse = calc_metrics(test.values, sarima_pred.values)
prophet_mae, prophet_rmse = calc_metrics(test.values, prophet_pred.values)

print("\n--- Model Comparison ---")
print(f"{'Model':<10} {'MAE':>10} {'RMSE':>10}")
print(f"{'SARIMA':<10} {sarima_mae:>10.2f} {sarima_rmse:>10.2f}")
print(f"{'Prophet':<10} {prophet_mae:>10.2f} {prophet_rmse:>10.2f}")

# Comparison plot
plt.figure(figsize=(14, 6))
plt.plot(test.index, test.values, label="Actual", color="black", linewidth=1.5)
plt.plot(test.index, sarima_pred.values, label=f"SARIMA (MAE={sarima_mae:.0f})", linestyle="--")
plt.plot(test.index, prophet_pred.values, label=f"Prophet (MAE={prophet_mae:.0f})", linestyle="--")
plt.fill_between(
    sarima_ci.index,
    sarima_ci["mean_ci_lower"],
    sarima_ci["mean_ci_upper"],
    alpha=0.1,
    label="SARIMA 95% CI",
)
plt.legend()
plt.title("Forecast Comparison: SARIMA vs Prophet (90-Day Holdout)")
plt.xlabel("Date")
plt.ylabel("Total Daily Sales")
plt.tight_layout()
plt.show()

You now have a quantified, visual comparison of two fundamentally different forecasting approaches on the same holdout data. This is exactly how forecasting model selection works in practice.

Secret mission

Add XGBoost with Feature Engineering

SARIMA and Prophet capture lags, seasonality, and trends automatically. But what if you had to build all of that by hand? In this secret mission, you will add an XGBoost model that requires explicit feature engineering for every pattern the statistical models get for free.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

  • Project folder containing demand_forecasting.ipynb, demand_forecasting_xgboost.ipynb, and the data/ subdirectory.
  • Python packages: statsmodels, prophet, matplotlib, jupyter, and xgboost (if you completed the Secret Mission).
  • Kaggle account (free, used to download the dataset).

Keep everything running

No action needed. Choose this if you plan to extend your forecasting pipeline with new models, datasets, or longer forecast horizons.

  • Your notebooks and data stay in your project folder, ready to open anytime with jupyter notebook.
  • The installed Python packages remain available for future projects.
  • Your Kaggle account stays active for downloading other datasets.

Pause - I'll come back to this later

Shut down the Jupyter Notebook server to free up memory. Your files stay on disk for when you return.

  • In the terminal where Jupyter is running, press Ctrl+C and confirm with y to shut down the notebook server.
  • Your notebooks and data remain saved in your project folder.
  • When you want to resume, navigate back to your project folder and run jupyter notebook again.

Delete - I don't want to use this again

Remove all project files and uninstall the packages you added.

  • First, shut down Jupyter if it is still running. Press Ctrl+C in its terminal and confirm with y.
  • Delete your project folder (the one containing your notebooks and the data/ subdirectory). On macOS/Linux:
cd ~/Desktop
rm -rf demand-forecasting

On Windows:

cd ~\Desktop
rmdir /s /q demand-forecasting
  • Uninstall the Python packages you installed for this project:
pip uninstall statsmodels prophet matplotlib jupyter xgboost -y
  • If you only completed the main project and not the Secret Mission, omit xgboost from the command above.

Nice Work!

Nice Work!

You just built a complete demand forecasting pipeline that takes raw retail sales data and produces 90-day forecasts using two fundamentally different modeling approaches.

You've learned how to:

  • Properly decomposed a time series into trend, seasonality, and residuals, tested for stationarity using the ADF test, and implemented a leakage-free temporal train/test split.
  • Trained and evaluated a SARIMA model with weekly seasonal differencing and a Prophet model with automatic seasonality detection.
  • Compared models quantitatively using MAE and RMSE and visually on the same holdout period to determine which approach captures the data patterns better.
  • Secret Mission: Added an XGBoost baseline with hand-crafted lag and rolling features, showing the explicit feature engineering that ML approaches require compared to what statistical models handle automatically.

Ready to quiz yourself?