Module 5 — ARIMA and SARIMA

Modules 3 and 4 handed us everything we need to fit a proper statistical model: a stationarized version of the pharmacy series, a reading of its ACF and PACF, and a baseline MAE of 9.8 units per day to beat. Now we open the ARIMA family — the classical work-horse of statistical forecasting — and use it to fill the second row of our results table.

Three letters, three ideas

The name ARIMA(p, d, q) unpacks into three components.

AR($p$) — autoregressive. The current value is a linear combination of the past $p$ values:

$$y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p} + \varepsilon_t$$

The $\phi_i$ are learned coefficients; the process is stationary only if their sum is below 1. AR(1) with $\phi_1$ close to 1 behaves like a random walk with slow reversion.

I($d$) — integrated. We difference the series $d$ times before modeling. In practice $d \in \{0, 1, 2\}$; $d = 1$ removes a linear trend, and we saw in module 3 that our pharmacy series needed $d = 1$.

MA($q$) — moving average. The current value depends on the past $q$ innovations (the residuals from previous predictions):

$$y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q}$$

The MA part captures short-lived shocks that decay quickly; AR captures long memory.

An ARIMA(p, d, q) model puts the three together on the differenced series.

Seasonal extension: SARIMA

Real series often have periodic structure that ARIMA alone cannot represent. SARIMA(p, d, q)(P, D, Q)$_s$ adds a seasonal AR of order $P$, a seasonal difference of order $D$ at period $s$, and a seasonal MA of order $Q$. On our pharmacy data with weekly seasonality, $s = 7$; on monthly retail data it would be 12; on hourly electricity data, 24 or 168.

The full six-tuple notation intimidates on first read; it is worth learning because it exactly matches the diagnostic we did in module 3: read ACF and PACF, pick $p, q, P, Q$, and check with AIC.

Fitting SARIMA with statsmodels

statsmodels gives us SARIMAX — SARIMA with optional exogenous regressors (which we will use in module 7).

import pandas as pd
from statsmodels.tsa.statespace.sarimax import SARIMAX

df = pd.read_parquet("pharmacy_daily.parquet")
series = df.set_index("date")["units"].asfreq("D")
train, test = series.iloc[:-28], series.iloc[-28:]

model = SARIMAX(
    train,
    order=(2, 1, 0),
    seasonal_order=(0, 1, 1, 7),
    enforce_stationarity=False,
    enforce_invertibility=False,
)
fit = model.fit(disp=False)
print(fit.summary().tables[1])

order=(2, 1, 0) says AR(2), one differencing, no MA. seasonal_order=(0, 1, 1, 7) says no seasonal AR, one seasonal differencing at lag 7, a seasonal MA of order 1. This is exactly what the ACF/PACF reading of module 3 suggested.

The two enforce_*=False flags let the fit converge on real data where the theoretical constraints are borderline. In production we would leave them on and inspect why the fit failed if it did.

AIC-based order selection

Reading ACF and PACF is a starting point, not a proof. The Akaike Information Criterion balances model fit against complexity:

$$\text{AIC} = 2k - 2\log \hat{L}$$

where $k$ is the number of parameters and $\hat{L}$ the maximum likelihood. Lower is better. We compare candidate orders by AIC on the same training window.

pmdarima.auto_arima scans a grid of orders and returns the AIC-optimal fit:

import pmdarima as pm

auto = pm.auto_arima(
    train,
    seasonal=True, m=7,
    d=1, D=1,
    start_p=0, max_p=3, start_q=0, max_q=3,
    start_P=0, max_P=2, start_Q=0, max_Q=2,
    stepwise=True, information_criterion="aic",
    trace=True, error_action="ignore", suppress_warnings=True,
)
print(auto.summary())

On our data, auto_arima picks SARIMA(1, 1, 1)(0, 1, 1)$_7$ with an AIC of 8 402. Manual reading suggested SARIMA(2, 1, 0)(0, 1, 1)$_7$ with AIC 8 411. Close, and their forecasts differ by about 0.2 units of MAE. That is a common outcome: the reading and the AIC search converge on a small neighborhood, and the differences among them are dwarfed by other decisions (are we differencing correctly? is there a change point?).

Diagnostics: is the fit trustworthy?

A SARIMA with .summary() looking pretty can still be wrong. Three checks decide.

Residual whiteness. The Ljung-Box test (module 3) is applied to the residuals: p-value above 0.05 means we did not leave autocorrelation on the floor.

from statsmodels.stats.diagnostic import acorr_ljungbox
resid = fit.resid[10:]  # skip burn-in
print(acorr_ljungbox(resid, lags=[10, 20], return_df=True))

Residual plot. fit.plot_diagnostics(figsize=(10, 6)) produces four panels: standardized residual over time (should look like white noise), histogram with a normal overlay, Q-Q plot, and residual ACF. Any spike outside the confidence bands in the ACF signals a missed structure.

Normality of residuals. ARIMA intervals assume normally distributed innovations. If the Q-Q plot bends heavily at the tails, the point forecast is still fine but the interval is optimistic — you will see the true value fall outside your "95 %" band more than 5 % of the time. Module 10 will replace those parametric intervals with quantile intervals for that reason.

Producing a forecast with an interval

fc = fit.get_forecast(steps=28)
mean = fc.predicted_mean
ci = fc.conf_int(alpha=0.05)

import matplotlib.pyplot as plt
plt.plot(train.index[-90:], train.iloc[-90:], label="history")
plt.plot(test.index, test.values, label="actual", marker="o")
plt.plot(mean.index, mean.values, label="forecast", marker="x")
plt.fill_between(ci.index, ci.iloc[:, 0], ci.iloc[:, 1], alpha=0.2)
plt.legend()

The forecast tracks the weekly pattern, the interval widens with the horizon (as any statistically correct interval should), and we compare its MAE to the baselines.

The result on the pharmacy series

Model	MAE (units/day)
naive	24.6
seasonal_naive_7	9.8
SARIMA(1,1,1)(0,1,1)$_7$	8.4

SARIMA beats seasonal-naïve by 1.4 units, roughly 14 %. That is a decent gain for a five-line model that already produces intervals for free. Whether it is worth deploying against a two-line seasonal-naïve depends on the operational cost of the improvement, and module 10 will make that call.

When ARIMA breaks

Three failure modes recur.

Long horizons. ARIMA is trained to minimize one-step-ahead error. Its 28-day forecast is a recursion, and errors compound. On very long horizons (say 90 days) SARIMA typically underperforms models that were trained to output a whole window at once — Prophet, gradient boosting with cyclical features, or deep networks with a direct multi-step output.

Multiple seasonalities. SARIMA supports one seasonal period. For daily data with both a weekly and a yearly cycle, the yearly one must be handled another way — Fourier terms as exogenous regressors (SARIMAX accepts them), or a move to Prophet or gradient boosting in modules 6 and 7.

Change points. SARIMA has no notion of an abrupt level shift. A promotional campaign that lifts sales by 20 % for a week will bleed into the model's estimated coefficients and bias its forecast for weeks after. Module 6 (Prophet) handles this with an explicit changepoint mechanism; module 7 handles it with a dummy variable.

Do not skip the diagnostics

It is tempting to look at the AIC, note that it went down, and ship. A SARIMA whose Ljung-Box p-value is below 0.05 is one whose residuals still contain autocorrelation, and every metric you compute — including that AIC — is compromised. Ten seconds of plot_diagnostics prevent two weeks of confusion.

Summary

• ARIMA(p, d, q) combines autoregression, integration and moving-average terms; SARIMA(p, d, q)(P, D, Q)$_s$ adds seasonal counterparts at period $s$.
• The reading from module 3 (ACF and PACF) plus AIC on a small grid pick reasonable orders; pmdarima.auto_arima automates the search.
• Diagnostics are not optional: Ljung-Box for residual whiteness, plot_diagnostics for tail behavior, before trusting the intervals.
• SARIMA beats our seasonal-naïve baseline by 14 % here; it degrades on long horizons, multiple seasonalities, and change points — the reasons we move to Prophet and gradient boosting next.

Next module: exponential smoothing and Prophet — two families that trade off some of SARIMA's rigor for flexibility on holidays, changepoints and business calendars.