はじめに
https://yuhi-sa.github.io/posts/20260226_arima/1/ では statsmodels.tsa.arima.model.ARIMA を使ってARMA/ARIMAモデルをPythonで推定する方法を解説し、関連記事欄で「https://yuhi-sa.github.io/posts/20260224_kalman_filter/1/ は状態空間モデルに基づく時系列推定手法であり、ARIMAとの比較が有益」と触れるにとどめていました。実は、statsmodelsのARIMA推定は内部でカルマンフィルタを使って正確なガウス尤度を計算しています。本記事では、ARMAモデルを状態空間表現に変換し、カルマンフィルタの予測誤差分解から尤度を導出し、フルスクラッチで最尤推定を実装します。そして、statsmodelsの推定結果と数値的に比較し、両者が高精度で一致することを確認します。
ARMAモデルの状態空間表現（Harvey表現）
ARMA(p, q)モデル
\[ 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} \]は、\(m = \max(p, q+1)\) 次元の状態ベクトルを使うと、以下のコンパニオン形式の状態空間モデルに書き換えられます（Harvey, 1993）。
\[ \boldsymbol{\alpha}_{t} = T \boldsymbol{\alpha}_{t-1} + R \varepsilon_{t}, \qquad y_t = Z \boldsymbol{\alpha}_t \tag{2} \]本記事ではARMA(1,1)（\(m = \max(1, 2) = 2\) ）を例に、具体的に
\[ 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} \]とおきます。状態の第1成分 \(\alpha_t^{(1)} = x_t\) 、第2成分 \(\alpha_t^{(2)} = \theta \varepsilon_t\) です。実際に式(2)(3)を展開すると
\[ \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} \]となり、式(1)のARMA(1,1)と一致することが確認できます。https://yuhi-sa.github.io/posts/20260224_kalman_filter/1/ の記法と対応させると、状態遷移行列 \(A=T\) 、観測行列 \(H=Z\) 、プロセスノイズ共分散 \(Q = R R^\top \sigma^2\) （ランク1）、観測ノイズ \(R_{\text{obs}} = 0\) （観測は状態の線形結合そのもので、追加の測定誤差を持たない）という特殊ケースになります。
定常初期分布：離散リアプノフ方程式
ARMAモデルが定常（\(|\phi| < 1\) など）であれば、状態ベクトルの無条件平均は \(\mathbf{0}\) 、無条件共分散 \(P_0\) は次の離散リアプノフ方程式の解です。
\[ P_0 = T P_0 T^\top + R R^\top \sigma^2 \tag{5} \]scipy.linalg.solve_discrete_lyapunov でこの行列方程式を直接解くことで、カルマンフィルタの初期共分散を正しく設定できます（ゼロ行列や大きな値で初期化すると、収束するまでのバーンイン期間で尤度が歪みます）。
カルマンフィルタによる正確な対数尤度
カルマンフィルタを状態空間モデル(2)に適用すると、各時刻の予測誤差（イノベーション） \(e_t = y_t - Z\boldsymbol{\alpha}_{t|t-1}\) とその分散 \(F_t = Z P_{t|t-1} Z^\top\) が得られます。ガウス性の仮定の下、観測列の同時密度は予測誤差分解（prediction error decomposition）により
\[ \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} \]と書けます。これがARMA/ARIMAモデルの正確な（近似なしの）ガウス対数尤度であり、statsmodelsを含む多くの統計ソフトウェアがこの式をカルマンフィルタで計算しています。式(6)を \(\boldsymbol{\theta} = (\phi, \theta, \sigma^2)\) について数値最適化すれば最尤推定量が得られます。
Python実装
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) # 定常初期共分散
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 # カルマンゲイン
a = a_pred + K * e
P = P_pred - np.outer(K, Z @ P_pred)
return -ll # 最小化のため符号反転
# MLE: 初期値 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) が式(5)を解いて定常初期共分散 \(P_0\)
を求め、ループ内の K = P_pred @ Z.T / F がカルマンゲイン、ll += ... が式(6)の対数尤度を逐次加算しています。
数値実験：statsmodelsとの一致検証
真のパラメータ \(\phi=0.7\)
, \(\theta=0.4\)
, \(\sigma^2=1.0\)
でARMA(1,1)データを500サンプル生成し、上記のフルスクラッチKF-MLEと statsmodels.tsa.arima.model.ARIMA(order=(1,0,1)) を同一データに適用しました。
| パラメータ | フルスクラッチKF | statsmodels | 差分 |
|---|---|---|---|
| \(\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⁻⁶ |
| 対数尤度 | -687.733333 | -687.733333 | 1.0×10⁻⁸ |
パラメータ推定値・対数尤度とも6桁の精度で一致しました（残差はNelder-Mead法とstatsmodelsの内部最適化アルゴリズムの違いによる数値誤差レベルで、理論的には同一の目的関数を最大化しています）。両者とも真値（\(\phi=0.7\) , \(\theta=0.4\) ）から有意にずれていますが、これは500サンプルという有限サンプルでの推定誤差であり、標準誤差の範囲内です（statsmodelsの出力では \(\phi\) の標準誤差は0.036）。
さらに、推定したパラメータを使って5期先までフィルタ状態から予測した値と、fit.forecast(5) の出力を比較すると、最大絶対誤差は \(3.7 \times 10^{-6}\)
でした。尤度・パラメータ・予測のすべてが一致しており、状態空間表現とカルマンフィルタがARIMA推定の正しい基盤であることが数値的に裏付けられます。
なぜこの理解が実務で役立つか
- 欠測値への対応：カルマンフィルタは観測が欠測した時刻でも更新ステップをスキップするだけで自然に扱えます。statsmodelsのARIMAが欠測値を含む系列でもそのまま
fit()できるのはこの性質によるものです。 - 状態空間モデルへの一般化：観測ノイズ \(R_{\text{obs}} > 0\) を導入すれば「ARMA信号 + 観測ノイズ」のモデルにそのまま拡張できます。これはhttps://yuhi-sa.github.io/posts/20260224_kalman_filter/1/の一般的な状態空間モデルの特殊ケースがARIMAである、という関係を直接示しています。
- 時変パラメータへの拡張：https://yuhi-sa.github.io/posts/20260715_rls_adaptive_filter/1/で見たRLSとカルマンフィルタの等価性と同様、ARMA係数を時変にしたい場合はカルマンフィルタの枠組みで自然に拡張できます（時変ARモデル、状態空間版のTVP-VARなど）。
関連記事
- ARIMA・SARIMAの理論とPython実装 - statsmodelsによるARIMA実装の基礎。本記事はその内部で使われる推定原理の解説です。
- GARCH(1,1)モデルの理論とPython実装 - 本記事と同じ予測誤差分解の枠組みを使い、平均ではなく条件付き分散をモデル化する姉妹記事です。
- カルマンフィルタの理論とPython実装 - 本記事の状態空間モデル・カルマンフィルタの記法の基盤です。
- RLS（逐次最小二乗法）アルゴリズムの理論とPython実装 - 同じくカルマンフィルタの特殊ケースとして適応フィルタを導出した姉妹記事です。
- PACF・ARモデル次数同定 - AR次数 \(p\) の同定手法。本記事の状態空間モデルの次元 \(m\) の決定に直結します。
- 時系列データの異常検知：統計的手法からカルマンフィルタまで - 状態空間モデルによる時系列異常検知の応用例です。
- 自己相関関数の理論とPython実装 - ARMAモデルの理論的自己相関構造の理解に役立ちます。
- 機械学習による時系列予測・分類・異常検知の選び方ハブ - ARIMA/Kalmanを含む7手法の選択指針です。
おすすめ書籍
ARMA/ARIMAモデルの理論を体系的に学べる定番書です。状態空間表現の章もあり、本記事の理論的背景を補完します。
※ 上記は Amazon アソシエイトのリンクです。
参考文献
- 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.