Stationarity and classical models

The big idea

A forecasting model learns patterns from the past and assumes they continue. That only works if the patterns are stable. A series whose average keeps climbing is a moving target: last year's average says nothing about next year's. The classical toolkit solves this in two steps. First transform the series until its behaviour no longer depends on when you look at it; then model the stable leftovers.

Think of predicting your commute time. The raw number keeps creeping up as the city grows, which is hard to model directly. But "how much longer than the same day last week" hovers around zero with a predictable spread. Model that difference, then add it back to last week's time. That is differencing, and it is half of ARIMA.

How it works

Step 1: Stationarity

A stationaryA time series whose statistical properties, such as mean, variance and autocorrelation, do not depend on the time at which it is observed. time series is, in the words of Forecasting: Principles and Practice, "one whose statistical properties do not depend on the time at which the series is observed". A series with a trend or seasonality is not stationary: its expected value depends on the date. Surprisingly, a series with irregular cycles can be stationary, because you cannot tell in advance where the peaks will fall.

The ACF from Time series fundamentals is a quick test: for non-stationary data it "decreases slowly", and $r_1$ is often large and positive.

Step 2: Differencing

DifferencingReplacing each value of a time series by its change from the previous value (or from the same season one period earlier), which removes trends or seasonal patterns. replaces each value by its change: $y_t^{\prime }=y_t-y_{t-1}$. A series 100, 103, 105, 109 becomes 3, 2, 4. A steady trend becomes a constant, which is stationary.

• Seasonal differencing compares each value with the same season one period earlier: $y_t^{\prime }=y_t-y_{t-m}$, such as this March minus last March.
• Second-order differencing differences twice, for trends that themselves change. The textbook notes "it is almost never necessary to go beyond second-order differences".
• If the differenced series is white noise, you have a random walkA series where each value is the previous value plus an unpredictable shock; its best forecast is the last observed value., $y_t=y_{t-1}+{\epsilon }_t$, and the best forecast is simply the last value.

Step 3: Testing for a unit root

How do you know if you have differenced enough? The ADF testAugmented Dickey-Fuller test: a hypothesis test whose null hypothesis is that a series has a unit root (is non-stationary in the random-walk sense). (augmented Dickey-Fuller) checks for a unit rootA property of a series, like a random walk, where shocks never die out, so the level wanders without returning to a fixed mean; differencing removes it., the random-walk kind of non-stationarity. The statsmodels documentation states: "The null hypothesis of the Augmented Dickey-Fuller is that there is a unit root, with the alternative that there is no unit root." So a small p-value ( Hypothesis testing and p-values ) is evidence of stationarity, and a large one means you cannot rule out a unit root. On the monthly CO2 series from the last lesson, the raw data give a p-value of 0.999; after one difference it drops below 0.001.

Step 4: AR, MA and ARIMA

Once the series is stationary, two simple ideas model what is left:

• An autoregressiveA model that predicts a value as a linear combination of the series' own previous p values plus noise. model AR(p) predicts from the series' own past values, like a regression ( Linear regression intuition ) on lags: $y_t=c+{\varphi }_1{y}_{t-1}+\cdots +{\varphi }_p{y}_{t-p}+{\epsilon }_t$
• A moving averageA time series model that predicts a value from the current and previous q forecast errors (shocks); unrelated to a moving-average smoother. model MA(q) predicts from the past forecast errors (shocks): $y_t=c+{\epsilon }_t+{\theta }_1{\epsilon }_{t-1}+\cdots +{\theta }_q{\epsilon }_{t-q}$. This MA has nothing to do with the moving-average smoother of decomposition.

ARIMAAutoRegressive Integrated Moving Average: an ARMA model fitted to a series differenced d times, written ARIMA(p, d, q). stands for AutoRegressive Integrated Moving Average, where "integration" is the reverse of differencing. ARIMA(p, d, q) differences the series d times, then fits AR(p) and MA(q) terms together. Many familiar models are special cases:

• ARIMA(0,0,0) with no constant: white noise.
• ARIMA(0,1,0): a random walk (with a constant: a random walk with drift).
• ARIMA(p,0,0): pure autoregression; ARIMA(0,0,q): pure moving average.

For seasonal data, SARIMASeasonal ARIMA, written ARIMA(p,d,q)(P,D,Q)m: adds seasonal differencing and AR and MA terms at multiples of the season length m. adds a second set of terms at multiples of the season length: ARIMA(p,d,q)(P,D,Q)m, where D is the number of seasonal differences and P and Q are AR and MA terms at lags m, 2m and so on.

Choosing orders. The partial autocorrelationThe correlation between a series and its value k steps earlier after removing the effects of the lags in between. (PACF) measures the link between $y_t$ and $y_{t-k}$ after removing the effects of lags 1 to k - 1. On the differenced data, a PACF with a significant spike at lag p and none beyond, with a decaying ACF, suggests ARIMA(p,d,0); an ACF that cuts off after lag q, with a decaying PACF, suggests ARIMA(0,d,q). When both p and q are positive the plots do not settle it, so in practice you compare candidates with an information criterion such as AIC, which automatic tools do for you.

Step 5: Exponential smoothing and Holt-Winters

A different classical family needs no differencing. Simple exponential smoothingForecasting with a weighted average of past observations whose weights decrease exponentially with age, controlled by a smoothing parameter alpha. forecasts with a weighted average of past values, where the weights decay exponentially:

With $\alpha =0.2$ the latest value gets weight 0.2, the one before 0.16 and the value five periods back still 0.0655. With $\alpha =0.8$ they are 0.8, 0.16 and 0.0003: the forecast follows the latest value closely. At $\alpha =1$ it is the naive forecast.

Holt added a smoothed trend, and Holt-WintersExponential smoothing with three smoothed components, level, trend and season, each with its own smoothing parameter. adds a smoothed season. The additive version keeps three running estimates, each updated by blending the new evidence with the old estimate:

The forecast h steps ahead is level plus h times trend plus the matching seasonal value. Use the additive form when seasonal swings are roughly constant, the multiplicative form when they grow with the level. A damped trend flattens out over long horizons; the textbook notes that a damped trend with multiplicative seasonality often gives accurate and robust forecasts for seasonal data.

Going deeper

Code: CO2, tested and forecast

import warnings; warnings.filterwarnings("ignore")
import numpy as np
from statsmodels.datasets import co2
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.holtwinters import ExponentialSmoothing

y = co2.load_pandas().data["co2"].resample("MS").mean().interpolate()   # monthly CO2, ppm

# 1. Stationarity: ADF test (null hypothesis: a unit root, i.e. not stationary)
for name, s in [("raw", y), ("diff", y.diff()), ("diff + seasonal diff", y.diff().diff(12))]:
    stat, p = adfuller(s.dropna(), autolag="AIC")[:2]
    print(f"ADF {name:22s} statistic {stat:7.2f}  p-value {p:.3f}")

# 2. Forecast the last 24 months from everything before
train, test = y[:"1999-12"], y["2000-01":]
h = len(test)
mae = lambda f: np.mean(np.abs(np.asarray(f) - test.values))

naive = np.tile(train.values[-12:], 2)                         # seasonal naive: same month last year
sarima = ARIMA(train, order=(1, 1, 1), seasonal_order=(0, 1, 1, 12)).fit()
hw = ExponentialSmoothing(train, trend="add", seasonal="add", seasonal_periods=12).fit()
last = np.full(h, train.values[-1])                            # naive: repeat the last value

print(f"\nMAE over {h} months (ppm):")
for name, f in [("naive (last value)", last), ("seasonal naive", naive),
                ("Holt-Winters additive", hw.forecast(h)), ("SARIMA(1,1,1)(0,1,1)12", sarima.forecast(h))]:
    print(f"  {name:24s} {mae(f):.2f}")
print(f"\nHolt-Winters smoothing: alpha {hw.params['smoothing_level']:.2f}, "
      f"beta {hw.params['smoothing_trend']:.2f}, gamma {hw.params['smoothing_seasonal']:.2f}")
ar, ma, sma = sarima.params[["ar.L1", "ma.L1", "ma.S.L12"]]
print(f"SARIMA coefficients: ar.L1 {ar:.2f}, ma.L1 {ma:.2f}, ma.S.L12 {sma:.2f}")
ci = sarima.get_forecast(h).conf_int(alpha=0.05).iloc[-1]
print(f"95% interval for 2001-12: {ci.iloc[0]:.1f} to {ci.iloc[1]:.1f} ppm (actual {test.iloc[-1]:.1f})")

ADF raw                    statistic    2.23  p-value 0.999
ADF diff                   statistic   -4.75  p-value 0.000
ADF diff + seasonal diff   statistic   -8.40  p-value 0.000

MAE over 24 months (ppm):
  naive (last value)       2.42
  seasonal naive           1.88
  Holt-Winters additive    0.24
  SARIMA(1,1,1)(0,1,1)12   0.35

Holt-Winters smoothing: alpha 0.57, beta 0.01, gamma 0.14
SARIMA coefficients: ar.L1 0.28, ma.L1 -0.60, ma.S.L12 -0.85
95% interval for 2001-12: 369.3 to 373.0 ppm (actual 371.0)

The ADF results confirm the picture: raw CO2 has a unit-root p-value of 0.999, while the differenced series and the doubly differenced series reject the unit root (p below 0.001).

The forecasting test trains on monthly data up to the end of 1999 and forecasts the 24 months of 2000 and 2001, scoring with mean absolute error in ppm:

• Naive (repeat the last value): 2.42. It ignores both the season and the trend.
• Seasonal naive (same month last year): 1.88. It captures the season but is a whole year of trend behind, about 1.5 to 2 ppm.
• Holt-Winters additive: 0.24.
• SARIMA(1,1,1)(0,1,1)12: 0.35.

The two classical models are five to eight times more accurate than the seasonal naive baseline on this series. The fitted Holt-Winters trend smoothing is tiny (${\beta }^{\ast }=0.01$), meaning the trend is very stable, which is what makes CO2 so forecastable. The SARIMA forecast also comes with a 95 percent prediction interval: for December 2001 it is 369.3 to 373.0 ppm, and the actual value, 371.0, falls inside. Prediction intervals matter as much as point forecasts ( Backtesting and forecast evaluation ).

Assumptions and limits

• ARIMA residuals should look like white noise. Check their ACF; leftover spikes mean the model missed structure.
• These are univariate models: they see only the series' own past. Promotions, prices and weather need regression terms (ARIMAX, dynamic regression) or the ML methods in ML and deep learning for forecasting .
• One model per series. With 50,000 products you fit 50,000 models, and new products with no history get nothing.
• Exponential smoothing and ARIMA overlap: some models in each family produce identical forecasts. Choose by out-of-sample error, not by family loyalty.

Trade-offs

Common mistakes

• Fitting ARIMA to a trending series without differencing (d = 0), then trusting the forecasts.
• Reading the ADF p-value backwards: a small p-value rejects the unit root, it does not prove non-stationarity.
• Over-differencing: differencing a series that is already stationary adds artificial negative autocorrelation.
• Confusing the MA model (past errors) with a moving-average smoother (past values).
• Judging a model by in-sample fit instead of forecasts on a later, held-out period.

In the wild

ARIMA and exponential smoothing remain the default baselines in demand planning, inventory control, capacity planning and economic statistics, because they are fast, need little data and have well-understood prediction intervals. Python's statsmodels provides ARIMA, SARIMAX and ExponentialSmoothing, R's fable package fits ETS and ARIMA models automatically across many series, and libraries such as statsforecast offer fast automatic versions. A team that cannot beat a tuned Holt-Winters or SARIMA model with a complex ML system usually has a data problem, not a model problem.

Check yourself

Recap