ARIMA

The model#

Let $L$ denote the lag operator, $L y_t = y_{t-1}$, and $\Delta = 1 - L$. An ARIMA$(p,d,q)$ process satisfies

$$ \phi(L)\,(1-L)^d\, y_t = c + \theta(L)\,\varepsilon_t $$

with

$$ \phi(L) = 1 - \phi_1 L - \dots - \phi_p L^p, \qquad \theta(L) = 1 + \theta_1 L + \dots + \theta_q L^q $$

• $y_t$ — the observation at time $t$.
• $\varepsilon_t$ — white noise, $\mathbb{E}[\varepsilon_t]=0$, $\mathrm{Var}(\varepsilon_t)=\sigma^2$, uncorrelated across $t$.
• $\phi_i$ — autoregressive coefficients; $p$ of them.
• $\theta_j$ — moving-average coefficients; $q$ of them.
• $d$ — the order of ordinary differencing.
• $c$ — a constant; when $d \ge 1$ it induces a deterministic polynomial trend of order $d$ in the levels.

The multiplicative seasonal extension SARIMA$(p,d,q)(P,D,Q)_m$ is

$$ \phi(L)\,\Phi(L^m)\,(1-L)^d (1-L^m)^D\, y_t = c + \theta(L)\,\Theta(L^m)\,\varepsilon_t $$

with $\Phi$ and $\Theta$ polynomials of degree $P$ and $Q$ in $L^m$, and $m$ the seasonal period.

Stationarity, invertibility and why they are enforced#

Stationarity of the AR part requires all roots of $\phi(z)=0$ to lie outside the unit circle in the complex plane. Equivalently, the companion matrix has spectral radius below one. For AR(1) this reduces to $|\phi_1|<1$; for AR(2) the admissible region is the triangle $\phi_1+\phi_2<1$, $\phi_2-\phi_1<1$, $|\phi_2|<1$.

Invertibility of the MA part requires all roots of $\theta(z)=0$ to lie outside the unit circle. Without it the process has no convergent AR($\infty$) representation, and — more damaging in practice — the parameters are not identified: for MA(1), $\theta$ and $1/\theta$ generate the same autocovariance function with rescaled $\sigma^2$. Enforcing invertibility picks one of the pair.

enforce_stationarity=False in the SARIMAX call above relaxes the first constraint during optimisation. It is useful when a seasonal difference has already made the series stationary and the optimiser wants to explore a boundary; it is not a licence to ignore a genuinely explosive fit.

Identification by ACF and PACF#

The autocorrelation function is $\rho(k) = \gamma(k)/\gamma(0)$. The partial autocorrelation $\alpha(k)$ is the last coefficient in a projection of $y_t$ on $y_{t-1},\dots,y_{t-k}$, computed by the Durbin-Levinson recursion.

The identification rules follow from the theory rather than from convention:

Under the null of white noise, $\hat{\rho}(k)$ is approximately $\mathrm{N}(0, 1/T)$ (Bartlett), giving the familiar $\pm 1.96/\sqrt{T}$ bands. These are pointwise, so with 20 displayed lags you expect one exceedance by chance.

Estimation#

The state-space form with the Kalman filter gives the exact Gaussian likelihood, including for series with missing observations. Writing $\mathbf{x}_t$ for the state,

$$ \mathbf{x}t = T\mathbf{x}{t-1} + R\eta_t, \qquad y_t = Z\mathbf{x}_t + \varepsilon_t $$

the prediction error decomposition yields

$$ \log \mathcal{L} = -\frac{1}{2}\sum_{t=1}^{T}\left[\log(2\pi F_t) + \frac{v_t^2}{F_t}\right] $$

with $v_t$ the one-step prediction error and $F_t$ its variance from the filter. Harvey (1989) and Durbin and Koopman (2012) give the construction. Non-stationary components are handled by exact diffuse initialisation (Koopman, 1997), which is what statsmodels uses when $d>0$.

Cost is $O(T r^3)$ for state dimension $r = \max(p+d+mP+mD,\; q+mQ+1)$. For monthly seasonal models $r$ is small; for high-frequency data with $m=336$ the state blows up, which is why SARIMA is impractical for sub-daily multiple seasonality and why Fourier-regressor approaches take over there.

Model selection#

$$ \mathrm{AIC} = -2\log\mathcal{L} + 2k, \qquad \mathrm{AICc} = \mathrm{AIC} + \frac{2k(k+1)}{T-k-1}, \qquad \mathrm{BIC} = -2\log\mathcal{L} + k\log T $$

with $k$ the number of estimated parameters including $\sigma^2$. AIC is asymptotically efficient for prediction; BIC is consistent for selecting a true finite-dimensional model. AICc corrects AIC's small-sample bias and should be the default when $T/k < 40$.

The critical restriction: the likelihoods must be computed on the same observations. Different $d$ or $D$ changes the effective sample after differencing, so those values are not comparable. Fix the differencing orders by test, then select $p, q, P, Q$ by criterion.

Hyndman and Khandakar (2008), Automatic time series forecasting: the forecast package for R, specifies the stepwise search that auto.arima and pmdarima implement: choose $D$ by a seasonal-strength or OCSB test, $d$ by KPSS, then a hill-climbing search over the remaining orders under AICc.

Residual diagnostics#

A correctly specified model leaves white-noise residuals. The Ljung-Box statistic

$$ Q = T(T+2)\sum_{k=1}^{h}\frac{\hat{\rho}k^2}{T-k} \;\sim\; \chi^2{h-k^*} $$

tests that jointly, where $k^*$ is the number of fitted ARMA parameters (the degrees of freedom must be reduced, a correction routinely omitted). Significant $Q$ means unmodelled structure remains and the prediction intervals are too narrow.

Heteroskedastic residuals with no autocorrelation in levels but strong autocorrelation in squares indicate conditional variance dynamics; ARCH and GARCH (Engle, 1982; Bollerslev, 1986) model the second moment and leave the ARIMA mean equation intact.

Where it sits today#

ARIMA remains the reference point rather than the state of the art. In M4, statistical methods including ARIMA formed the backbone of the top hybrid entries but no pure ARIMA entry placed near the top. In M5, gradient-boosted trees on engineered features dominated the hierarchical retail setting, where cross-series information and calendar covariates matter more than a single series' autocorrelation.

ARIMA's enduring advantages are honest ones: closed-form prediction intervals with known coverage properties under its assumptions, few parameters, no tuning infrastructure, and interpretable coefficients. Those matter when you must defend a forecast to a regulator.

Reading#

• Box, Jenkins, Reinsel and Ljung, Time Series Analysis: Forecasting and Control, 5th ed., Wiley 2015 — the source of the methodology and of the airline model.
• Hamilton, Time Series Analysis, Princeton 1994, chapters 3 to 5.
• Durbin and Koopman, Time Series Analysis by State Space Methods, 2nd ed., OUP 2012.
• Hyndman and Khandakar (2008), Automatic time series forecasting, Journal of Statistical Software 27(3).
• Ljung and Box (1978), On a measure of lack of fit in time series models, Biometrika 65(2).

What to learn next#

• Prophet — a different design that handles holidays and missing data without hand-tuned orders.
• Evaluating a forecast — scoring ARIMA against the baselines it must beat.
• Multivariate forecasting — bringing outside information into the model.