Introduction
https://yuhi-sa.github.io/en/posts/20260226_arima/1/ covered fitting ARMA/ARIMA models in Python with statsmodels.tsa.arima.model.ARIMA, noting in its related-articles section only that “https://yuhi-sa.github.io/en/posts/20260224_kalman_filter/1/ is a state-space estimation method, and comparing it to ARIMA is worthwhile” — without elaborating. In fact, statsmodels’ ARIMA estimation internally computes the exact Gaussian likelihood using a Kalman filter. This article converts an ARMA model into state-space form, derives the likelihood from the Kalman filter’s prediction-error decomposition, and implements maximum likelihood estimation from scratch. We then compare the result numerically against statsmodels and confirm the two agree to high precision.
State-Space Representation of ARMA Models (Harvey Form)
The ARMA(p, q) model
\[ x_t = \phi_1 x_{t-1} + \cdots + \phi_p x_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}, \qquad \varepsilon_t \sim \mathcal{N}(0, \sigma^2) \tag{1} \]can be rewritten as a state-space model in companion form using an \(m = \max(p, q+1)\) -dimensional state vector (Harvey, 1993):
\[ \boldsymbol{\alpha}_{t} = T \boldsymbol{\alpha}_{t-1} + R \varepsilon_{t}, \qquad y_t = Z \boldsymbol{\alpha}_t \tag{2} \]For ARMA(1,1) (\(m = \max(1, 2) = 2\) ), this becomes concretely
\[ T = \begin{bmatrix} \phi & 1 \\ 0 & 0 \end{bmatrix}, \qquad R = \begin{bmatrix} 1 \\ \theta \end{bmatrix}, \qquad Z = \begin{bmatrix} 1 & 0 \end{bmatrix} \tag{3} \]The state’s first component is \(\alpha_t^{(1)} = x_t\) and the second is \(\alpha_t^{(2)} = \theta \varepsilon_t\) . Expanding equations (2)-(3) directly gives
\[ \alpha_t^{(1)} = \phi \alpha_{t-1}^{(1)} + \alpha_{t-1}^{(2)} + \varepsilon_t = \phi x_{t-1} + \theta \varepsilon_{t-1} + \varepsilon_t \tag{4} \]which matches the ARMA(1,1) equation (1) exactly. Mapping onto the notation from https://yuhi-sa.github.io/en/posts/20260224_kalman_filter/1/: transition matrix \(A=T\) , observation matrix \(H=Z\) , process-noise covariance \(Q = R R^\top \sigma^2\) (rank 1), and observation noise \(R_{\text{obs}} = 0\) (the observation is exactly a linear combination of the state, with no additional measurement error) — a special case of the general Kalman-filter setup.
The Stationary Initial Distribution: A Discrete Lyapunov Equation
If the ARMA model is stationary (e.g., \(|\phi| < 1\) ), the state vector’s unconditional mean is \(\mathbf{0}\) , and its unconditional covariance \(P_0\) solves the discrete Lyapunov equation
\[ P_0 = T P_0 T^\top + R R^\top \sigma^2 \tag{5} \]Solving this matrix equation directly with scipy.linalg.solve_discrete_lyapunov gives the Kalman filter its correct initial covariance (initializing with a zero matrix or a large arbitrary value distorts the likelihood during the burn-in period before convergence).
Exact Log-Likelihood via the Kalman Filter
Applying the Kalman filter to state-space model (2) produces, at each time step, the prediction error (innovation) \(e_t = y_t - Z\boldsymbol{\alpha}_{t|t-1}\) and its variance \(F_t = Z P_{t|t-1} Z^\top\) . Under the Gaussian assumption, the joint density of the observations factors via the prediction error decomposition as
\[ \log L(\boldsymbol{\theta}) = -\frac{1}{2}\sum_{t=1}^{N}\left[\log(2\pi F_t) + \frac{e_t^2}{F_t}\right] \tag{6} \]This is the exact (not approximate) Gaussian log-likelihood for the ARMA/ARIMA model, and it’s what statsmodels and most other statistical software compute via the Kalman filter under the hood. Numerically maximizing (6) over \(\boldsymbol{\theta} = (\phi, \theta, \sigma^2)\) yields the maximum likelihood estimate.
Python Implementation
import numpy as np
from scipy.optimize import minimize
from scipy.linalg import solve_discrete_lyapunov
def kf_loglik(params, y):
phi, theta, log_sigma2 = params
sigma2 = np.exp(log_sigma2)
T = np.array([[phi, 1.0], [0.0, 0.0]])
R = np.array([1.0, theta])
Z = np.array([1.0, 0.0])
Q = np.outer(R, R) * sigma2
P = solve_discrete_lyapunov(T, Q) # stationary initial covariance
a = np.zeros(2)
ll = 0.0
for t in range(len(y)):
a_pred = T @ a
P_pred = T @ P @ T.T + Q
F = Z @ P_pred @ Z.T
if F <= 0:
return 1e10
e = y[t] - Z @ a_pred
ll += -0.5 * (np.log(2 * np.pi * F) + e**2 / F)
K = P_pred @ Z.T / F # Kalman gain
a = a_pred + K * e
P = P_pred - np.outer(K, Z @ P_pred)
return -ll # negate for minimization
# MLE, starting from phi=0, theta=0, log(sigma^2)=0
x0 = np.array([0.0, 0.0, 0.0])
res = minimize(kf_loglik, x0, args=(x,), method="Nelder-Mead",
options={"xatol": 1e-8, "fatol": 1e-8, "maxiter": 5000})
phi_hat, theta_hat, sigma2_hat = res.x[0], res.x[1], np.exp(res.x[2])
solve_discrete_lyapunov(T, Q) solves equation (5) for the stationary initial covariance \(P_0\)
; K = P_pred @ Z.T / F inside the loop is the Kalman gain, and ll += ... accumulates the log-likelihood from equation (6).
Numerical Experiment: Agreement with statsmodels
We simulated 500 samples of an ARMA(1,1) process with true parameters \(\phi=0.7\)
, \(\theta=0.4\)
, \(\sigma^2=1.0\)
, then fit both the from-scratch KF-MLE above and statsmodels.tsa.arima.model.ARIMA(order=(1,0,1)) to the same data.
| Parameter | From-scratch KF | statsmodels | Difference |
|---|---|---|---|
| \(\phi\) | 0.683662 | 0.683659 | 2.6×10⁻⁶ |
| \(\theta\) | 0.488372 | 0.488371 | 1.2×10⁻⁶ |
| \(\sigma^2\) | 0.914022 | 0.914016 | 6.1×10⁻⁶ |
| Log-likelihood | -687.733333 | -687.733333 | 1.0×10⁻⁸ |
Parameter estimates and log-likelihood agree to six decimal places. The residual differences are numerical-optimization noise from Nelder-Mead vs. statsmodels’ internal optimizer — both are maximizing the same objective function. Both estimates deviate noticeably from the true values (\(\phi=0.7\) , \(\theta=0.4\) ), which is ordinary finite-sample estimation error at \(N=500\) (statsmodels reports a standard error of 0.036 for \(\phi\) , so the estimate is well within one standard error of the truth).
We also compared 5-step-ahead forecasts computed from the fitted filtered state against fit.forecast(5): the maximum absolute difference was \(3.7 \times 10^{-6}\)
. Likelihood, parameters, and forecasts all agree, numerically confirming that the state-space representation and Kalman filter are the correct computational foundation for ARIMA estimation.
Why This Understanding Matters in Practice
- Handling missing observations: The Kalman filter naturally skips the update step at any time index with a missing observation. This is exactly why statsmodels’ ARIMA can fit a series with missing values without any special preprocessing.
- Generalizing to state-space models: Introducing observation noise \(R_{\text{obs}} > 0\) immediately extends the model to “ARMA signal plus measurement noise.” This directly demonstrates that ARIMA is a special case of the general state-space model covered in https://yuhi-sa.github.io/en/posts/20260224_kalman_filter/1/.
- Extending to time-varying parameters: Just as https://yuhi-sa.github.io/en/posts/20260715_rls_adaptive_filter/1/ showed RLS is a special case of the Kalman filter, allowing the ARMA coefficients to vary over time (time-varying AR models, state-space TVP-VAR) fits naturally into the same Kalman-filter framework.
Related Articles
- ARIMA and SARIMA: Theory and Python Implementation - The foundation for fitting ARIMA models in statsmodels. This article explains the estimation machinery working underneath it.
- GARCH(1,1) in Python - A companion article using the same prediction-error-decomposition framework to model conditional variance instead of the mean.
- Kalman Filter: Theory and Python Implementation - The state-space model and Kalman filter notation used throughout this article.
- The RLS Algorithm in Python: Recursive Least Squares and Its Equivalence to the Kalman Filter - A companion article deriving adaptive filtering as another special case of the Kalman filter.
- PACF and AR Order Identification - Methods for choosing the AR order \(p\) , directly relevant to sizing the state dimension \(m\) here.
- Time Series Anomaly Detection: From Statistical Methods to Kalman Filters - An application of state-space models to anomaly detection.
- Autocorrelation Function: Theory and Python Implementation - Useful background on the theoretical autocorrelation structure of ARMA models.
- Choosing Among Time-Series Methods: Forecasting, Classification, and Anomaly Detection - A selection guide across seven methods including ARIMA and Kalman filters.
References
- Harvey, A. C. (1993). Time Series Models (2nd ed.). MIT Press. Chapter 3 (State space form of ARMA models).
- Hamilton, J. D. (1994). Time Series Analysis. Princeton University Press. Chapter 13 (The Kalman Filter).
- Durbin, J., & Koopman, S. J. (2012). Time Series Analysis by State Space Methods (2nd ed.). Oxford University Press.
- Seabold, S., & Perktold, J. (2010). statsmodels: Econometric and statistical modeling with python. Proceedings of the 9th Python in Science Conference.