Forecast Machine Replacements with SARIMAX
Build a SARIMAX pipeline to forecast weekly machine replacements.
Introduction
30 Second Summary
When machines break down unexpectedly, maintenance teams scramble to find replacements. Predicting how many replacements you'll need each week turns reactive firefighting into proactive planning.
In this project, you will build a time series forecasting pipeline that predicts weekly machine replacements using external drivers like machine age, environmental conditions, and maintenance history.
What You'll Build
A complete SARIMAX forecasting notebook in Python that incorporates exogenous variables into weekly replacement predictions and produces multi-week forecasts with confidence intervals.
By the end of this project, you'll have:
- A fitted SARIMAX model that uses machine age, environmental conditions, and maintenance history as predictive drivers.
- Forecast accuracy metrics (RMSE, MAE, MAPE) evaluated on a held-out test set.
- Multi-week forecasts with confidence intervals that can inform maintenance planning decisions.
- Secret Mission: Compare SARIMAX against a SARIMA baseline to quantify exactly how much predictive value your external drivers add.
Are there any prerequisites?
You should be familiar with Python and pandas. A basic understanding of time series concepts like trend and seasonality will also help.
Set Up Your Environment
Before you can build your SARIMAX model, you need the right tools in place. This step installs the statistical modeling libraries and sets up your notebook with all the imports and data you'll need for the rest of the project.
You'll work with two key libraries: statsmodels provides the SARIMAX implementation for fitting time series models, and pmdarima automates the tricky process of selecting model orders.
In this step, get ready to:
- Install statsmodels and pmdarima.
- Create your forecasting notebook and import all required libraries.
- Load your weekly replacement data and verify it's ready for modeling.
Install statsmodels and pmdarima
These two libraries form the core of your forecasting pipeline. statsmodels contains the SARIMAX model class, stationarity tests, and decomposition tools. pmdarima wraps statsmodels with an automated order selection algorithm that saves you from manual trial-and-error.
- Install both libraries by running this command in your terminal:
pip install statsmodels==0.14.6 pmdarima==2.1.1- Verify the installation succeeded by running this quick import check:
python -c "from statsmodels.tsa.statespace.sarimax import SARIMAX; import pmdarima as pm; print('Libraries installed successfully')"- You should see Libraries installed successfully printed to your terminal.
Seeing an install error?
If you see a permissions error, try adding --user to the end of the pip command.
If you see a version conflict, make sure you're using the pip associated with your Python 3 installation. Try pip3 instead of pip.
Still stuck? Help me troubleshoot my statsmodels or pmdarima install error.
Create your notebook and set up imports
Your forecasting work lives in a single Python file that uses # %% cell markers. These markers let you run sections independently, just like cells in a Jupyter Notebook.
- Create a new file called forecasting_notebook.py in your project folder.
- Add the following imports to your file:
# %% [markdown]
# # Weekly Machine Replacement Forecasting with SARIMAX
# %% Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error, mean_absolute_error
import pmdarima as pm
plt.style.use('seaborn-v0_8-whitegrid')What does this code do?
- pandas and numpy handle your data manipulation and numerical operations.
- matplotlib creates all the visualizations: time series plots, decomposition charts, and forecast graphs.
- The statsmodels imports bring in SARIMAX for modeling, adfuller for stationarity testing, seasonal_decompose for breaking down your series, and plot_acf/plot_pacf for autocorrelation analysis.
- pmdarima (imported as pm) automates ARIMA order selection with its auto_arima function.
- plt.style.use('seaborn-v0_8-whitegrid') sets a clean grid background for all your charts.
- Save forecasting_notebook.py and run the imports cell.
- You should see no output and no errors. A clean run means all libraries loaded correctly.
Seeing a ModuleNotFoundError?
This means the library is not installed in the Python environment your notebook is using. Make sure the pip you used to install matches the Python interpreter running your notebook.
If you're using a virtual environment, confirm it's activated before running the cell.
Help me fix my ModuleNotFoundError for statsmodels or pmdarima.
Load your weekly replacement data
Your dataset is a CSV file with weekly timestamps and four columns: the target replacement count plus three exogenous variables. SARIMAX cannot handle missing values in exogenous columns, so you'll verify the data is complete before proceeding.
- Make sure your weekly_replacements.csv file is in the same folder as your notebook.
- Add a new cell below your imports with the following data loading code:
# %% Load and prepare data
# Replace this section with your actual data loading
df = pd.read_csv('weekly_replacements.csv', parse_dates=['week'], index_col='week')
df.index.freq = 'W'
# Verify columns
print("Columns:", df.columns.tolist())
print("Shape:", df.shape)
print("Date range:", df.index.min(), "to", df.index.max())
print("\nMissing values:")
print(df.isnull().sum())
# Target and exogenous variables
target_col = 'replacements'
exog_cols = ['avg_machine_age', 'env_condition', 'maintenance_events']What does this code do?
- pd.read_csv(...) loads the CSV and parses the week column as datetime, setting it as the index.
- df.index.freq = 'W' explicitly sets the frequency to weekly. SARIMAX requires a known frequency on the DatetimeIndex.
- The print statements let you confirm the data shape, date range, and that no values are missing.
- target_col and exog_cols store column references you'll reuse throughout the notebook.
- Save forecasting_notebook.py and run the data loading cell.
- You should see your column names (replacements, avg_machine_age, env_condition, maintenance_events), the shape, the date range, and 0 for all missing value counts.
Seeing a FileNotFoundError?
This means Python cannot find weekly_replacements.csv in the current working directory. Make sure the CSV file is in the same folder as forecasting_notebook.py.
If your missing values column shows any non-zero counts, you'll need to handle those before fitting. SARIMAX will throw an error if any exogenous values are NaN.
Help me debug my data loading issue.
Why do exogenous variables need to be complete?
SARIMAX treats exogenous variables as OLS regressors alongside the SARIMA error structure. If any value is missing, the model cannot compute the regression coefficients and will fail during fitting.
This is different from the target variable, where SARIMAX can sometimes handle gaps through its state-space formulation.
✔️ Awesome, I've got everything!
Great. Double check you've saved your file before moving on.
ⓧ I'd like to double check the full code
Here's the complete forecasting_notebook.py so far:
# %% [markdown]
# # Weekly Machine Replacement Forecasting with SARIMAX
# %% Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error, mean_absolute_error
import pmdarima as pm
plt.style.use('seaborn-v0_8-whitegrid')
# %% Load and prepare data
# Replace this section with your actual data loading
df = pd.read_csv('weekly_replacements.csv', parse_dates=['week'], index_col='week')
df.index.freq = 'W'
# Verify columns
print("Columns:", df.columns.tolist())
print("Shape:", df.shape)
print("Date range:", df.index.min(), "to", df.index.max())
print("\nMissing values:")
print(df.isnull().sum())
# Target and exogenous variables
target_col = 'replacements'
exog_cols = ['avg_machine_age', 'env_condition', 'maintenance_events']Your environment is ready and your data is loaded. Next up, you'll explore the time series visually and test whether it's stationary enough for SARIMAX to model directly.
Explore Data and Test for Stationarity
Your data is loaded and ready to analyze. Before fitting any model, you need to understand the structure hiding inside your time series.
Time series models like SARIMAX make assumptions about your data. The most critical assumption is stationarity. A stationary series has constant statistical properties (mean and variance) over time. If your data violates this, the model's predictions will be unreliable.
In this step, get ready to:
- Visualize the replacement series and exogenous variables to spot trends and seasonal patterns.
- Test stationarity with the Augmented Dickey-Fuller test and identify differencing needs.
- Check correlations between your exogenous variables and the target to confirm they add predictive value.
Visualize the time series and exogenous variables
Plotting your data is the fastest way to spot trends, seasonal cycles, and potential relationships between variables. You will create a multi-panel chart that shows the replacement count alongside each exogenous driver.
- Add a new code cell in your notebook and create the 4-panel time series plot by pasting this code:
# %% Visualize the time series
fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)
axes[0].plot(df.index, df[target_col], color='steelblue', linewidth=1.2)
axes[0].set_ylabel('Replacements')
axes[0].set_title('Weekly Machine Replacements')
for i, col in enumerate(exog_cols):
axes[i + 1].plot(df.index, df[col], color='darkorange', linewidth=1.0)
axes[i + 1].set_ylabel(col)
plt.tight_layout()
plt.show()What does this code do?
This creates a 4-row figure with a shared x-axis (time). The top panel shows your target variable (weekly replacements) in blue. The three panels below it each show one exogenous variable in orange.
By stacking them vertically with sharex=True, you can visually scan for periods where a spike in one exogenous variable lines up with a spike in replacements.
- Run the cell. You should see a 4-panel chart with the replacement count on top and each exogenous variable below it.
Not seeing the chart?
Make sure you ran the imports cell earlier (the one containing import matplotlib.pyplot as plt). Also confirm that target_col and exog_cols are defined from Step 1. If you restarted the kernel, re-run all previous cells first.
Help me debug why my time series plot isn't showing
Now you need to separate the trend, seasonality, and residual noise from your replacement series. This is called seasonal decomposition.
- Add a new cell and run the seasonal decomposition by pasting this code:
# %% Seasonal decomposition
decomposition = seasonal_decompose(df[target_col], model='additive', period=52)
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.tight_layout()
plt.show()What is seasonal decomposition?
Seasonal decomposition splits your series into three components: a long-term trend, a repeating seasonal pattern, and leftover residual noise.
The period=52 parameter tells the function your data repeats every 52 weeks (one year). The model='additive' assumes the components add together rather than multiply. This decomposition helps you confirm whether seasonality is strong enough to justify the seasonal terms in SARIMAX.
- Run the cell. You should see a 4-panel decomposition showing the observed series, trend, seasonal pattern, and residuals.
Seeing a ValueError about period?
If your dataset has fewer than 104 observations (less than 2 full years), period=52 may cause an error. Try reducing the period to 13 (quarterly seasonality) if you have a shorter dataset.
Help me fix a period error in seasonal_decompose
Test stationarity with the Augmented Dickey-Fuller test
Visual inspection gives you clues, but you need a formal statistical test to confirm whether your series is stationary. The Augmented Dickey-Fuller (ADF) test checks for the presence of a unit root. If a unit root exists, the series is non-stationary and needs differencing before modeling.
- Add a new cell and define the ADF test helper function, then run it on your replacement series:
# %% Augmented Dickey-Fuller test
def adf_test(series, name=''):
result = adfuller(series.dropna(), autolag='AIC')
print(f"ADF Test for {name}:")
print(f" Test Statistic: {result[0]:.4f}")
print(f" p-value: {result[1]:.4f}")
print(f" Lags Used: {result[2]}")
print(f" Critical Values: {result[4]}")
if result[1] < 0.05:
print(f" => Series is stationary (reject unit root)\n")
else:
print(f" => Series is NOT stationary (cannot reject unit root)\n")
return result[1]
p_value = adf_test(df[target_col], name=target_col)
# If non-stationary, difference and re-test
if p_value >= 0.05:
df['replacements_diff'] = df[target_col].diff()
adf_test(df['replacements_diff'].dropna(), name='replacements (differenced)')What does this code do?
The adf_test() function wraps the adfuller test from statsmodels. It prints the test statistic, p-value, and a human-readable verdict. If p < 0.05, you reject the null hypothesis of a unit root, meaning the series IS stationary.
The conditional block at the bottom applies first differencing (.diff()) if the original series fails the stationarity test. Differencing subtracts each value from its predecessor, removing trends. The re-test confirms whether one round of differencing was enough.
Why does stationarity matter for SARIMAX?
SARIMAX models the autocorrelation in a series. If the mean or variance is drifting over time, the model cannot distinguish real patterns from trend artifacts. The "d" in ARIMA(p,d,q) represents how many times the series needs to be differenced to become stationary.
The ADF test result tells you whether d should be 0 (already stationary) or 1+ (needs differencing). Later, auto_arima will select this automatically, but understanding it here builds your intuition.
- Run the cell. You should see the ADF test output with a test statistic, p-value, and verdict printed for your replacements series.
- If the series is non-stationary (p >= 0.05), you will also see a second test result on the differenced series.
Getting a NameError for adfuller?
Make sure your imports cell includes from statsmodels.tsa.stattools import adfuller and that you have run it. Restart your kernel and re-run all cells from the top if needed.
Help me fix my ADF test setup
Next, you will plot the ACF and PACF to visually identify candidate AR and MA orders for your model.
- Add a new cell and generate the ACF/PACF plots by pasting this code:
# %% ACF and PACF plots
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(df[target_col].dropna(), lags=40, ax=axes[0])
plot_pacf(df[target_col].dropna(), lags=40, ax=axes[1])
plt.tight_layout()
plt.show()What are ACF and PACF?
ACF (Autocorrelation Function) shows how correlated the series is with lagged versions of itself. PACF (Partial Autocorrelation Function) strips out the effect of intermediate lags.
The number of significant lags in PACF suggests the AR order (p). The number of significant lags in ACF suggests the MA order (q). You will use these visual clues to sanity-check the automated order selection in the next step.
- Run the cell. You should see two side-by-side bar plots showing correlation values at each lag, with a blue shaded confidence band.
Seeing a warning about the PACF method?
This warning is harmless and does not affect your results. The PACF plot will still render correctly. If you want to suppress it, you can add import warnings; warnings.filterwarnings('ignore') at the top of the cell.
Help me understand my ACF/PACF plots
Check exogenous variable correlations
Before feeding exogenous variables into SARIMAX, you need to verify two things. First, do they actually correlate with your target? Second, are they collinear with each other? Perfectly collinear exogenous variables cause rank deficiency errors during model fitting.
- Add a new cell and compute the correlation matrices by pasting this code:
# %% Correlation check
print("Correlation with target:")
print(df[exog_cols + [target_col]].corr()[target_col].drop(target_col))
print("\nExogenous variable correlation matrix:")
print(df[exog_cols].corr())What should I look for?
The first output shows how each exogenous variable correlates with your replacement count. Values closer to 1 or -1 indicate stronger relationships. Any variable with near-zero correlation may not add much to your model.
The second output shows correlations between exogenous variables themselves. If two exogenous variables have a correlation above 0.9, they carry nearly identical information. This collinearity can cause numerical instability when SARIMAX tries to estimate separate coefficients for each.
- Run the cell. You should see correlation values printed for each exogenous variable against replacements, followed by a 3x3 correlation matrix between the exogenous variables.
Seeing very high collinearity (above 0.9)?
If two of your exogenous variables are highly correlated, consider dropping one before fitting the model. You can remove a column from exog_cols and re-run subsequent steps with the reduced set.
Help me decide which exogenous variable to drop
✔️ Awesome, I've got everything!
Great. Double check you have saved your notebook with all cells run.
ⓧ I'd like to double check the full code
Here is the complete notebook code up to this point. Make sure your cells match:
# %% [markdown]
# # Weekly Machine Replacement Forecasting with SARIMAX
# %% Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error, mean_absolute_error
import pmdarima as pm
plt.style.use('seaborn-v0_8-whitegrid')
# %% Load and prepare data
# Replace this section with your actual data loading
df = pd.read_csv('weekly_replacements.csv', parse_dates=['week'], index_col='week')
df.index.freq = 'W'
# Verify columns
print("Columns:", df.columns.tolist())
print("Shape:", df.shape)
print("Date range:", df.index.min(), "to", df.index.max())
print("\nMissing values:")
print(df.isnull().sum())
# Target and exogenous variables
target_col = 'replacements'
exog_cols = ['avg_machine_age', 'env_condition', 'maintenance_events']
# %% Visualize the time series
fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)
axes[0].plot(df.index, df[target_col], color='steelblue', linewidth=1.2)
axes[0].set_ylabel('Replacements')
axes[0].set_title('Weekly Machine Replacements')
for i, col in enumerate(exog_cols):
axes[i + 1].plot(df.index, df[col], color='darkorange', linewidth=1.0)
axes[i + 1].set_ylabel(col)
plt.tight_layout()
plt.show()
# %% Seasonal decomposition
decomposition = seasonal_decompose(df[target_col], model='additive', period=52)
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.tight_layout()
plt.show()
# %% Augmented Dickey-Fuller test
def adf_test(series, name=''):
result = adfuller(series.dropna(), autolag='AIC')
print(f"ADF Test for {name}:")
print(f" Test Statistic: {result[0]:.4f}")
print(f" p-value: {result[1]:.4f}")
print(f" Lags Used: {result[2]}")
print(f" Critical Values: {result[4]}")
if result[1] < 0.05:
print(f" => Series is stationary (reject unit root)\n")
else:
print(f" => Series is NOT stationary (cannot reject unit root)\n")
return result[1]
p_value = adf_test(df[target_col], name=target_col)
# If non-stationary, difference and re-test
if p_value >= 0.05:
df['replacements_diff'] = df[target_col].diff()
adf_test(df['replacements_diff'].dropna(), name='replacements (differenced)')
# %% ACF and PACF plots
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(df[target_col].dropna(), lags=40, ax=axes[0])
plot_pacf(df[target_col].dropna(), lags=40, ax=axes[1])
plt.tight_layout()
plt.show()
# %% Correlation check
print("Correlation with target:")
print(df[exog_cols + [target_col]].corr()[target_col].drop(target_col))
print("\nExogenous variable correlation matrix:")
print(df[exog_cols].corr())You now have a clear picture of your data's structure. You know whether differencing is needed, you can see the seasonal patterns, and you have confirmed your exogenous variables relate to the target. Next up, you will use these insights to fit the actual SARIMAX model with automated order selection.
Fit SARIMAX with Automated Order Selection
You've confirmed whether your replacement series is stationary and identified which exogenous variables correlate with it. Now it's time to build the actual forecasting model.
In this step, you'll split the data into training and test sets, use pmdarima to automatically search for the best ARIMA orders, then fit the full SARIMAX model with statsmodels. By the end, you'll have a trained model with interpretable coefficients showing how each exogenous variable influences weekly replacements.
In this step, get ready to:
- Split your data into training and test sets.
- Run auto_arima to find optimal ARIMA and seasonal orders.
- Fit and interpret the SARIMAX model with statsmodels.
Split data into train and test sets
Before fitting any model, you need a held-out test set to evaluate forecast accuracy on data the model has never seen. Holding out the final 8 weeks gives you enough observations to calculate meaningful error metrics later.
- In your Jupyter Notebook, create a new cell and add the train/test split by running this code:
# %% Train/test split
test_weeks = 8
y_train = df[target_col].iloc[:-test_weeks]
y_test = df[target_col].iloc[-test_weeks:]
X_train = df[exog_cols].iloc[:-test_weeks]
X_test = df[exog_cols].iloc[-test_weeks:]
print(f"Training: {len(y_train)} weeks ({y_train.index.min()} to {y_train.index.max()})")
print(f"Testing: {len(y_test)} weeks ({y_test.index.min()} to {y_test.index.max()})")
What does this code do?
This splits your data into two aligned pieces using the last 8 weeks as the test set.
- y_train and y_test are the target variable (replacement counts) split at -8 weeks.
- X_train and X_test are the exogenous variable DataFrames aligned to the same split point.
- The print statements confirm how many weeks of training data you have and the date range of your test set.
- Run the cell. You should see two lines showing the number of training weeks and testing weeks, along with their date ranges.
Seeing a KeyError?
- Make sure target_col and exog_cols are defined in a cell above (they were set when you loaded your data).
- If you restarted your kernel, re-run all cells from the top to restore your variables.
Ask for help: Why am I getting a KeyError when splitting my training and test data?
Run auto_arima to find optimal orders
A SARIMAX model requires you to choose (p,d,q) for the non-seasonal component and (P,D,Q,s) for the seasonal component. Instead of manually trying every combination, auto_arima performs a stepwise search that tests candidates and selects the one with the lowest AIC (Akaike Information Criterion).
- Create a new cell and run the automated order selection:
# %% Auto ARIMA with exogenous variables
print("Running auto_arima (this may take a few minutes for weekly data)...")
auto_model = pm.auto_arima(
y_train,
X=X_train,
seasonal=True,
m=52,
stepwise=True,
suppress_warnings=True,
trace=True,
error_action='ignore',
max_p=3,
max_q=3,
max_P=2,
max_Q=2,
max_d=2,
max_D=1,
)
print(f"\nBest model: ARIMA{auto_model.order} x {auto_model.seasonal_order}")
print(f"AIC: {auto_model.aic():.2f}")
What does this code do?
- X=X_train passes your exogenous variables into the search so that candidate models are evaluated WITH external drivers included.
- m=52 tells the search that your data has a 52-week seasonal cycle (one year).
- stepwise=True uses an efficient algorithm that tests a subset of combinations instead of brute-forcing every possibility.
- trace=True prints each candidate model and its AIC as the search progresses so you can watch it work.
- The max_p, max_q, max_P, and max_Q limits keep the search space manageable for weekly data.
The final print statement shows the winning (p,d,q) x (P,D,Q,s) combination and its AIC score. Lower AIC means a better balance of model fit and complexity.
Why use X instead of exogenous?
In pmdarima 2.0+, the exogenous parameter was removed and replaced with X. Using the old name raises a TypeError. Always use X for exogenous data in pmdarima.
- Run the cell. This may take 2-5 minutes depending on your dataset size. You should see lines printing each model candidate with its AIC score, followed by the best order at the end.
Search is running very slowly?
- If you have fewer than 104 observations (less than 2 years of weekly data), m=52 can cause convergence issues. Try changing m=52 to m=13 (quarterly seasonality) for shorter datasets.
- If the cell seems stuck, wait it out. The stepwise search with m=52 is computationally intensive for weekly seasonal patterns.
Ask for help: My auto_arima search is taking too long or failing to converge. What should I try?
Fit and interpret the SARIMAX model
Now that auto_arima has found the best orders, you'll fit the model using statsmodels directly. This gives you a full statistical summary with coefficient significance tests for each exogenous variable.
- Create a new cell and fit the SARIMAX model by running this code:
# %% Fit SARIMAX with statsmodels
order = auto_model.order
seasonal_order = auto_model.seasonal_order
model = SARIMAX(
y_train,
exog=X_train,
order=order,
seasonal_order=seasonal_order,
enforce_stationarity=False,
enforce_invertibility=False,
)
results = model.fit(disp=False)
print(results.summary())
What does this code do?
- order and seasonal_order are extracted from the auto_arima results so you use the exact same (p,d,q)(P,D,Q,s) combination it selected.
- enforce_stationarity=False and enforce_invertibility=False prevent numerical errors when the optimizer explores edge cases during fitting.
- results.summary() prints the full model table including coefficient values, standard errors, and p-values for each exogenous variable.
In the summary output, look at the p-values for avg_machine_age, env_condition, and maintenance_events. A p-value below 0.05 means that variable has a statistically significant effect on weekly replacements. Also check that coefficient signs make domain sense (e.g., higher machine age should increase replacements).
- Run the cell. You should see a large summary table showing model parameters, coefficient estimates, and diagnostic statistics.
Getting a convergence warning?
- A convergence warning means the optimizer couldn't find the best parameters within its iteration limit. The model still fits, but results may be less reliable.
- This is common with m=52 seasonal models. The enforce_stationarity=False setting helps avoid most of these issues.
Ask for help: My SARIMAX model is showing convergence warnings. Should I be worried?
Now verify the model's residuals pass diagnostic checks. Good residuals look like white noise with no remaining patterns.
- Create a new cell and generate the diagnostics plot by running:
# %% Model diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()
What should I see in the diagnostics?
The diagnostics plot shows four panels. The top-left shows standardized residuals over time. The top-right is a histogram with a normal distribution curve overlay.
The bottom-left is a Q-Q plot. Points close to the diagonal line mean residuals are approximately normal. The bottom-right shows autocorrelation (correlogram) of residuals. Values within the shaded band mean no significant leftover patterns in your residuals.
- Run the cell. You should see a 4-panel diagnostic figure. Confirm that the residuals look roughly random (no obvious trends in the top-left panel) and the correlogram stays mostly within the shaded confidence bands.
Residuals show a pattern?
- If the top-left panel shows a clear trend or cycles, the model may be missing seasonal or autoregressive structure. Consider re-running auto_arima with higher max_P or max_Q limits.
- Significant spikes in the correlogram (bottom-right) at specific lags suggest the model hasn't captured all autocorrelation. This can happen if your data has multiple seasonal cycles.
Ask for help: My SARIMAX diagnostics show patterns in the residuals. How do I improve the model?
✔️ 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
# %% [markdown]
# # Weekly Machine Replacement Forecasting with SARIMAX
# %% Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error, mean_absolute_error
import pmdarima as pm
plt.style.use('seaborn-v0_8-whitegrid')
# %% Load and prepare data
# Replace this section with your actual data loading
df = pd.read_csv('weekly_replacements.csv', parse_dates=['week'], index_col='week')
df.index.freq = 'W'
# Verify columns
print("Columns:", df.columns.tolist())
print("Shape:", df.shape)
print("Date range:", df.index.min(), "to", df.index.max())
print("\nMissing values:")
print(df.isnull().sum())
# Target and exogenous variables
target_col = 'replacements'
exog_cols = ['avg_machine_age', 'env_condition', 'maintenance_events']
# %% Visualize the time series
fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)
axes[0].plot(df.index, df[target_col], color='steelblue', linewidth=1.2)
axes[0].set_ylabel('Replacements')
axes[0].set_title('Weekly Machine Replacements')
for i, col in enumerate(exog_cols):
axes[i + 1].plot(df.index, df[col], color='darkorange', linewidth=1.0)
axes[i + 1].set_ylabel(col)
plt.tight_layout()
plt.show()
# %% Seasonal decomposition
decomposition = seasonal_decompose(df[target_col], model='additive', period=52)
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.tight_layout()
plt.show()
# %% Augmented Dickey-Fuller test
def adf_test(series, name=''):
result = adfuller(series.dropna(), autolag='AIC')
print(f"ADF Test for {name}:")
print(f" Test Statistic: {result[0]:.4f}")
print(f" p-value: {result[1]:.4f}")
print(f" Lags Used: {result[2]}")
print(f" Critical Values: {result[4]}")
if result[1] < 0.05:
print(f" => Series is stationary (reject unit root)\n")
else:
print(f" => Series is NOT stationary (cannot reject unit root)\n")
return result[1]
p_value = adf_test(df[target_col], name=target_col)
# If non-stationary, difference and re-test
if p_value >= 0.05:
df['replacements_diff'] = df[target_col].diff()
adf_test(df['replacements_diff'].dropna(), name='replacements (differenced)')
# %% ACF and PACF plots
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(df[target_col].dropna(), lags=40, ax=axes[0])
plot_pacf(df[target_col].dropna(), lags=40, ax=axes[1])
plt.tight_layout()
plt.show()
# %% Correlation check
print("Correlation with target:")
print(df[exog_cols + [target_col]].corr()[target_col].drop(target_col))
print("\nExogenous variable correlation matrix:")
print(df[exog_cols].corr())
# %% Train/test split
test_weeks = 8
y_train = df[target_col].iloc[:-test_weeks]
y_test = df[target_col].iloc[-test_weeks:]
X_train = df[exog_cols].iloc[:-test_weeks]
X_test = df[exog_cols].iloc[-test_weeks:]
print(f"Training: {len(y_train)} weeks ({y_train.index.min()} to {y_train.index.max()})")
print(f"Testing: {len(y_test)} weeks ({y_test.index.min()} to {y_test.index.max()})")
# %% Auto ARIMA with exogenous variables
print("Running auto_arima (this may take a few minutes for weekly data)...")
auto_model = pm.auto_arima(
y_train,
X=X_train,
seasonal=True,
m=52,
stepwise=True,
suppress_warnings=True,
trace=True,
error_action='ignore',
max_p=3,
max_q=3,
max_P=2,
max_Q=2,
max_d=2,
max_D=1,
)
print(f"\nBest model: ARIMA{auto_model.order} x {auto_model.seasonal_order}")
print(f"AIC: {auto_model.aic():.2f}")
# %% Fit SARIMAX with statsmodels
order = auto_model.order
seasonal_order = auto_model.seasonal_order
model = SARIMAX(
y_train,
exog=X_train,
order=order,
seasonal_order=seasonal_order,
enforce_stationarity=False,
enforce_invertibility=False,
)
results = model.fit(disp=False)
print(results.summary())
# %% Model diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()
Your SARIMAX model is trained and validated. Next up, you'll use it to generate multi-week forecasts with confidence intervals and measure exactly how accurate those predictions are against your held-out test data.
Forecast and Evaluate Accuracy
Your SARIMAX model is fitted and diagnosed. Now it is time to answer the question this entire pipeline was built for: how many machine replacements should you expect in the coming weeks?
In this step, you will generate multi-week forecasts with confidence intervals, visualize how well your predictions track reality, and compute accuracy metrics that quantify your model's performance on unseen data.
In this step, get ready to:
- Generate forecasts with 95% confidence intervals on the held-out test set.
- Visualize forecast predictions against actual replacement counts.
- Compute RMSE, MAE, and MAPE to evaluate forecast accuracy.
Generate forecasts with confidence intervals
The get_forecast method produces point predictions plus a probability distribution around each point. You must supply the future exogenous values because SARIMAX cannot generate predictions without knowing what the external drivers will be during the forecast horizon.
- Add a new cell in your notebook and generate the forecast by running this code:
# %% Forecast on test set
forecast = results.get_forecast(steps=len(y_test), exog=X_test)
predictions = forecast.predicted_mean
conf_int = forecast.conf_int(alpha=0.05)What does this code do?
results.get_forecast(steps=len(y_test), exog=X_test) asks the model to predict the next 8 weeks, feeding in the matching exogenous values the model needs for each future time step.
forecast.predicted_mean extracts the single best-guess value for each week. forecast.conf_int(alpha=0.05) extracts the 95% confidence bounds, giving you a range that should contain the true value 95% of the time.
- Run the cell.
- You should see the cell execute without errors. The predictions variable now holds 8 weekly forecast values.
Seeing a ValueError about exog shape?
The number of rows in X_test must equal the steps argument. Confirm that len(X_test) == len(y_test) by printing both lengths. If they differ, revisit your train/test split from the previous step.
Help me debug a SARIMAX exog shape mismatch error
Visualize forecast vs actuals
A chart that overlays your predictions on the real data is the fastest way to judge whether the model captures the overall pattern. The confidence band shows where uncertainty grows as you forecast further into the future.
- Add a new cell and create the forecast visualization by running this code:
# %% Visualize forecast vs actuals
fig, ax = plt.subplots(figsize=(14, 6))
# Plot training data (last 30 weeks for visibility)
ax.plot(y_train.index[-30:], y_train.iloc[-30:], label='Training Data', color='steelblue')
# Plot test actuals
ax.plot(y_test.index, y_test, label='Actual', color='black', linewidth=2)
# Plot predictions
ax.plot(predictions.index, predictions, label='SARIMAX Forecast', color='red', linestyle='--', linewidth=2)
# Confidence interval
ax.fill_between(
conf_int.index,
conf_int.iloc[:, 0],
conf_int.iloc[:, 1],
alpha=0.2,
color='red',
label='95% Confidence Interval'
)
ax.axvline(x=y_test.index[0], color='gray', linestyle=':', label='Train/Test Split')
ax.set_xlabel('Week')
ax.set_ylabel('Machine Replacements')
ax.set_title('SARIMAX Forecast vs Actuals')
ax.legend()
plt.tight_layout()
plt.show()What does this visualization show?
The chart has four layers. The blue line is the last 30 weeks of training data for context. The black line is the actual test data your model has never seen. The red dashed line is the model's forecast.
The shaded red band is the 95% confidence interval. Wider bands mean more uncertainty. A vertical gray dotted line marks where the train/test split occurs so you can clearly see where predictions begin.
- Run the cell.
- You should see a line chart with the training data in blue, actual test values in black, the red dashed forecast line, and a light red shaded confidence band.
Chart appears blank or axes are missing?
Make sure you ran the forecast cell above first. The predictions and conf_int variables must exist before this cell executes. If you restarted the kernel, re-run all cells from the top.
Help me figure out why my forecast chart is blank
Compute accuracy metrics
Visual inspection tells you the shape of the forecast. Quantitative metrics tell you exactly how far off the predictions are. Three standard metrics give you different perspectives on forecast error.
- Add a new cell and compute the evaluation metrics by running this code:
# %% Evaluation metrics
rmse = np.sqrt(mean_squared_error(y_test, predictions))
mae = mean_absolute_error(y_test, predictions)
mape = np.mean(np.abs((y_test - predictions) / y_test)) * 100
print(f"SARIMAX Forecast Evaluation:")
print(f" RMSE: {rmse:.2f}")
print(f" MAE: {mae:.2f}")
print(f" MAPE: {mape:.2f}%")
print(f" AIC: {results.aic:.2f}")What do these metrics tell you?
- RMSE (Root Mean Squared Error) penalizes large errors heavily. It is in the same units as your target (number of replacements).
- MAE (Mean Absolute Error) gives the average magnitude of errors without squaring. It is easier to interpret as "on average, the forecast is off by X replacements."
- MAPE (Mean Absolute Percentage Error) expresses error as a percentage of actual values. A MAPE under 10% is generally considered a strong forecast.
- AIC (Akaike Information Criterion) measures model fit while penalizing complexity. Lower AIC means a better balance of accuracy and simplicity.
- Run the cell.
- You should see four printed metric values: RMSE, MAE, MAPE (as a percentage), and AIC.
Seeing inf or NaN in your MAPE?
MAPE divides by y_test, so any week where actual replacements equal zero causes a division-by-zero. Check y_test.min() to see if you have zero values. If so, filter them out before computing MAPE or use MAE as your primary metric instead.
Help me fix inf values in my MAPE calculation
✔️ Awesome, I've got everything!
Great. Double check you have saved your notebook with all the cells from this step.
ⓧ I'd like to double check the full code
Here is the complete notebook code through the end of this step:
# %% [markdown]
# # Weekly Machine Replacement Forecasting with SARIMAX
# %% Imports
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_squared_error, mean_absolute_error
import pmdarima as pm
plt.style.use('seaborn-v0_8-whitegrid')
# %% Load and prepare data
# Replace this section with your actual data loading
df = pd.read_csv('weekly_replacements.csv', parse_dates=['week'], index_col='week')
df.index.freq = 'W'
# Verify columns
print("Columns:", df.columns.tolist())
print("Shape:", df.shape)
print("Date range:", df.index.min(), "to", df.index.max())
print("\nMissing values:")
print(df.isnull().sum())
# Target and exogenous variables
target_col = 'replacements'
exog_cols = ['avg_machine_age', 'env_condition', 'maintenance_events']
# %% Visualize the time series
fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)
axes[0].plot(df.index, df[target_col], color='steelblue', linewidth=1.2)
axes[0].set_ylabel('Replacements')
axes[0].set_title('Weekly Machine Replacements')
for i, col in enumerate(exog_cols):
axes[i + 1].plot(df.index, df[col], color='darkorange', linewidth=1.0)
axes[i + 1].set_ylabel(col)
plt.tight_layout()
plt.show()
# %% Seasonal decomposition
decomposition = seasonal_decompose(df[target_col], model='additive', period=52)
fig = decomposition.plot()
fig.set_size_inches(14, 8)
plt.tight_layout()
plt.show()
# %% Augmented Dickey-Fuller test
def adf_test(series, name=''):
result = adfuller(series.dropna(), autolag='AIC')
print(f"ADF Test for {name}:")
print(f" Test Statistic: {result[0]:.4f}")
print(f" p-value: {result[1]:.4f}")
print(f" Lags Used: {result[2]}")
print(f" Critical Values: {result[4]}")
if result[1] < 0.05:
print(f" => Series is stationary (reject unit root)\n")
else:
print(f" => Series is NOT stationary (cannot reject unit root)\n")
return result[1]
p_value = adf_test(df[target_col], name=target_col)
# If non-stationary, difference and re-test
if p_value >= 0.05:
df['replacements_diff'] = df[target_col].diff()
adf_test(df['replacements_diff'].dropna(), name='replacements (differenced)')
# %% ACF and PACF plots
fig, axes = plt.subplots(1, 2, figsize=(14, 4))
plot_acf(df[target_col].dropna(), lags=40, ax=axes[0])
plot_pacf(df[target_col].dropna(), lags=40, ax=axes[1])
plt.tight_layout()
plt.show()
# %% Correlation check
print("Correlation with target:")
print(df[exog_cols + [target_col]].corr()[target_col].drop(target_col))
print("\nExogenous variable correlation matrix:")
print(df[exog_cols].corr())
# %% Train/test split
test_weeks = 8
y_train = df[target_col].iloc[:-test_weeks]
y_test = df[target_col].iloc[-test_weeks:]
X_train = df[exog_cols].iloc[:-test_weeks]
X_test = df[exog_cols].iloc[-test_weeks:]
print(f"Training: {len(y_train)} weeks ({y_train.index.min()} to {y_train.index.max()})")
print(f"Testing: {len(y_test)} weeks ({y_test.index.min()} to {y_test.index.max()})")
# %% Auto ARIMA with exogenous variables
print("Running auto_arima (this may take a few minutes for weekly data)...")
auto_model = pm.auto_arima(
y_train,
X=X_train,
seasonal=True,
m=52,
stepwise=True,
suppress_warnings=True,
trace=True,
error_action='ignore',
max_p=3,
max_q=3,
max_P=2,
max_Q=2,
max_d=2,
max_D=1,
)
print(f"\nBest model: ARIMA{auto_model.order} x {auto_model.seasonal_order}")
print(f"AIC: {auto_model.aic():.2f}")
# %% Fit SARIMAX with statsmodels
order = auto_model.order
seasonal_order = auto_model.seasonal_order
model = SARIMAX(
y_train,
exog=X_train,
order=order,
seasonal_order=seasonal_order,
enforce_stationarity=False,
enforce_invertibility=False,
)
results = model.fit(disp=False)
print(results.summary())
# %% Model diagnostics
results.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()
# %% Forecast on test set
forecast = results.get_forecast(steps=len(y_test), exog=X_test)
predictions = forecast.predicted_mean
conf_int = forecast.conf_int(alpha=0.05)
# %% Visualize forecast vs actuals
fig, ax = plt.subplots(figsize=(14, 6))
# Plot training data (last 30 weeks for visibility)
ax.plot(y_train.index[-30:], y_train.iloc[-30:], label='Training Data', color='steelblue')
# Plot test actuals
ax.plot(y_test.index, y_test, label='Actual', color='black', linewidth=2)
# Plot predictions
ax.plot(predictions.index, predictions, label='SARIMAX Forecast', color='red', linestyle='--', linewidth=2)
# Confidence interval
ax.fill_between(
conf_int.index,
conf_int.iloc[:, 0],
conf_int.iloc[:, 1],
alpha=0.2,
color='red',
label='95% Confidence Interval'
)
ax.axvline(x=y_test.index[0], color='gray', linestyle=':', label='Train/Test Split')
ax.set_xlabel('Week')
ax.set_ylabel('Machine Replacements')
ax.set_title('SARIMAX Forecast vs Actuals')
ax.legend()
plt.tight_layout()
plt.show()
# %% Evaluation metrics
rmse = np.sqrt(mean_squared_error(y_test, predictions))
mae = mean_absolute_error(y_test, predictions)
mape = np.mean(np.abs((y_test - predictions) / y_test)) * 100
print(f"SARIMAX Forecast Evaluation:")
print(f" RMSE: {rmse:.2f}")
print(f" MAE: {mae:.2f}")
print(f" MAPE: {mape:.2f}%")
print(f" AIC: {results.aic:.2f}")Secret mission
Quantify Exogenous Variable Value
You built a SARIMAX model with exogenous variables, but how do you know those extra variables actually help? In this secret mission, you'll fit a SARIMA baseline without exogenous inputs and compare it head-to-head against your full model to quantify exactly how much predictive value machine age, environmental conditions, and maintenance history add.
Clean Up Your Resources
Clean Up Your Resources
This project runs entirely locally, so there are no ongoing costs or cloud resources to worry about.
Decide whether to keep your files for future re-training, pause by closing your notebook, or delete everything entirely.
Resources you used:
- Jupyter Notebook files (forecasting_notebook.py and forecasting_notebook_secret_mission.py).
- Data file (weekly_replacements.csv).
- Installed Python packages (statsmodels and pmdarima).
Keep everything running
No action needed. Choose this if you plan to re-train your model as new weekly data arrives.
- Your notebook and data files stay on your machine, ready for future use.
- As new weekly observations become available, refit your SARIMAX model periodically to keep forecasts current.
- The installed packages (statsmodels, pmdarima) remain available for other time series projects.
Pause - I'll come back to this later
Shut down running processes to free up memory. Your files stay saved on disk.
- Save your notebook if you haven't already.
- Shut down the Jupyter Notebook kernel by selecting Kernel then Shut Down in the menu bar.
- Close the Jupyter Notebook browser tab.
- In your terminal, stop the Jupyter server by pressing Ctrl+C.
Delete - I don't want to use this again
Remove all project files and uninstall packages to start completely fresh.
- Shut down the Jupyter Notebook kernel and stop the Jupyter server.
- Delete your project files by running:
rm forecasting_notebook.py forecasting_notebook_secret_mission.py weekly_replacements.csv- On Windows, use this command instead:
del forecasting_notebook.py, forecasting_notebook_secret_mission.py, weekly_replacements.csv- Uninstall the Python packages by running:
pip uninstall statsmodels pmdarima -y- This removes statsmodels and pmdarima from your Python environment. Other packages like pandas and numpy remain untouched.
Nice Work!
Nice Work!
You just built a complete time series forecasting pipeline that predicts weekly machine replacements using external drivers like machine age, environmental conditions, and maintenance history.
You've learned how to:
- Fitted a SARIMAX model using statsmodels and pmdarima, incorporating exogenous variables as predictive drivers and interpreting their statistical significance.
- Generated multi-week forecasts with 95% confidence intervals that can inform real maintenance planning decisions.
- Evaluated forecast accuracy using RMSE, MAE, and MAPE on held-out test data, verifying your model produces reliable predictions.
- Secret Mission: Quantified the added predictive value of exogenous variables by benchmarking SARIMAX against a SARIMA baseline.
Ready to quiz yourself?