はじめに
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 でこの行列方程式を直接解くことで、カルマンフィルタの初期共分散を正しく設定できます（ゼロ行列や大きな値で初期化すると、収束するまでのバーンイン期間で尤度が歪みます。実測値は後述の数値実験の節で示します）。
注意点：非定常・単位根近傍でのリアプノフ方程式の破綻
solve_discrete_lyapunov が返す \(P_0\)
が意味を持つのは \(|\phi| < 1\)
という定常性が成り立つときだけです。この境界がどれほど脆いか、\(\theta=0.4\)
, \(\sigma^2=1\)
を固定して \(\phi\)
を単位円に近づけながら実際に呼び出すと、次のようになります。
| \(\phi\) | solve_discrete_lyapunov の結果 \(P_0[0,0]\) | 備考 |
|---|---|---|
| 0.99 | 98.09 | 定常だが単位根に近く分散が急増 |
| 1.00 | （例外） | LinAlgError: A singular matrix detected |
| 1.01 | -97.91 | 非定常なのに例外を出さず負の「分散」を返す |
\(\phi=1.01\)
（爆発的過程）では solve_discrete_lyapunov は例外を送出せず、対角成分が負という共分散行列としてありえない行列を黙って返します。幸い、この値をそのままカルマンフィルタに渡すと、初回の予測分散も同じ符号で負になります（実測値 \(F_1 = -97.91\)
、\(\phi=0.99\)
側の巨大な正の値とは対照的に符号が反転しているだけで、絶対値はほぼ同じです）。本記事のPython実装（後述）にある if F <= 0: return 1e10 というガードは、単なる数値安全策ではなく、Nelder-Mead法などの数値最適化が非定常領域に迷い込んだ場合に目的関数へ大きなペナルティを与えて弾き返すという具体的な役割を果たしています。\(\phi\)
を \((-1, 1)\)
に制約する再パラメータ化（例えば \(\phi = \tanh(u)\)
）を使えば、そもそも探索空間からこの領域を除外できます。
なお、ARIMA(\(p, d, q\)
) のように \(d > 0\)
の差分を含むモデルでは、差分後のARMA部分は定常であっても、差分前の元の系列を状態に含める定式化では単位根成分そのものが状態に現れるため、式(5)の定常共分散が存在しません。実務ではこの非定常成分に「拡散事前分布（diffuse prior）」、すなわち理論的には無限大の分散を与えて初期化し、収束後の尤度計算からその影響を除外する手法が使われます（Rのstats::arimaやKFASパッケージの厳密拡散初期化などが該当）。本記事はこの発展的な話題には立ち入らず、\(d=0\)
の定常ARMA部分に絞って解説します。
カルマンフィルタによる正確な対数尤度
カルマンフィルタを状態空間モデル(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)\) について数値最適化すれば最尤推定量が得られます。
注意点：観測ノイズが厳密にゼロの場合の縮退
式(3)のARMA状態空間表現には測定誤差が存在しません（\(R_{\text{obs}}=0\) ）。このため、更新後の状態共分散 \(P_{t\vert t}\) の \((1,1)\) 成分（\(y_t\) に対応する状態成分の分散）は、理論上も数値上も厳密に \(0\) に潰れます。実際に推定パラメータ（後述の \(\hat\phi=0.676518\) , \(\hat\theta=0.420873\) , \(\hat\sigma^2=0.959929\) ）でカルマンフィルタを実行すると
| \(t\) | \(P_{t\vert t}[0,0]\) | \(\hat\alpha_{t\vert t}^{(1)}\) （更新後の状態） | \(y_t\) （観測値） |
|---|---|---|---|
| 0 | 0.0 | 0.4967141530 | 0.4967141530 |
| 1 | 0.0 | 0.4081212671 | 0.4081212671 |
| 2 | 0.0 | 0.8780677046 | 0.8780677046 |
のように、更新後のフィルタ状態は観測値 \(y_t\) と（表示桁数の範囲で）完全に一致します。観測が状態の第1成分を誤差なく明かすので、これは当然の帰結です。一見「フィルタが何もしていない」ように見えますが、実際に情報を運んでいるのは更新前の予測 \(a_{t\vert t-1}\) とその分散 \(F_{t\vert t-1}\) の方であり、これが式(6)の尤度と後述の数値実験の予測区間を構成します。カルマンフィルタの「更新後」の値だけを見ると、このモデルでは自明な結果しか得られない点に注意してください。
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\)
で statsmodels.tsa.arima_process.ArmaProcess から（再現性のため np.random.seed(42) を固定して）ARMA(1,1)データを500サンプル生成し、上記のフルスクラッチKF-MLEと statsmodels.tsa.arima.model.ARIMA(order=(1,0,1)) を同一データに適用しました。
注意点：ARIMAのデフォルト定数項に要注意
statsmodelsの ARIMA(y, order=(1,0,1)) はデフォルトで定数項（trend="c"）を推定します。本記事の状態空間表現(2)(3)はゼロ平均のARMAを前提としているため、定数項を含めたまま比較すると、ar.L1 係数とフルスクラッチKFの \(\hat\phi\)
が対応していても、統計的に有意でない小さな定数項（実測で約 \(0.035\)
）の分だけ尤度がずれてしまいます（筆者自身、最初の実行では対数尤度の差が \(1.7 \times 10^{-2}\)
生じ、原因の特定に手間取りました）。ゼロ平均モデルと正しく比較するには、trend="n" を明示する必要があります。
fit = ARIMA(x, order=(1, 0, 1), trend="n").fit()
以下の結果はすべて trend="n" に統一した上での比較です。
| パラメータ | フルスクラッチKF | statsmodels | 差分 |
|---|---|---|---|
| \(\phi\) | 0.676518 | 0.676520 | 1.8×10⁻⁶ |
| \(\theta\) | 0.420873 | 0.420870 | 2.5×10⁻⁶ |
| \(\sigma^2\) | 0.959929 | 0.959929 | 1.2×10⁻⁷ |
| 対数尤度 | -699.899146 | -699.899146 | 1.2×10⁻⁹ |
パラメータ推定値・対数尤度とも6桁以上の精度で一致しました（残差はNelder-Mead法とstatsmodelsの内部最適化アルゴリズムの違いによる数値誤差レベルで、理論的には同一の目的関数を最大化しています）。両者とも真値（\(\phi=0.7\) , \(\theta=0.4\) ）から有意にずれていますが、これは500サンプルという有限サンプルでの推定誤差であり、標準誤差の範囲内です（statsmodelsの出力では \(\phi\) の標準誤差は \(0.0392\) 、\(\theta\) は \(0.0492\) ）。
さらに、推定したパラメータを使って5期先までフィルタ状態から予測した値と、fit.forecast(5) の出力を比較すると、最大絶対誤差は \(5.8 \times 10^{-6}\)
でした。尤度・パラメータ・予測のすべてが一致しており、状態空間表現とカルマンフィルタがARIMA推定の正しい基盤であることが数値的に裏付けられます。
図示：一期先予測と初期共分散 \(P_0\) の選び方が尤度に与える影響
上図(A)は、推定パラメータで計算した一期先予測 \(Z a_{t\vert t-1}\) とその予測区間 \(\pm 2\sqrt{F_{t\vert t-1}}\) を観測系列に重ねたものです。予測分散は定常状態で \(\sqrt{F_\infty} \approx 0.980\) に収束しますが、\(t=1\) の時点（初期共分散 \(P_0\) をそのまま使う最初の予測）だけ \(\sqrt{F_1} \approx 1.758\) とやや広く、その後急速に定常値へ落ち着きます。
上図(B)は、本記事冒頭で「ゼロ行列や大きな値で初期化すると尤度が歪む」と述べた内容を実際の数値で確認したものです。3通りの \(P_0\) で計算した対数尤度の総和は、リアプノフ方程式の解（本記事の手法）で \(-699.899146\) 、ゼロ行列初期化で \(-699.343158\) 、拡散事前分布（\(10^6 I\) ）初期化で \(-706.422969\) となり、最大で \(6.5\) 程度の差が生じます。\(t=1\) の1時点だけを見ると差はさらに顕著で、対数尤度の項はそれぞれ \(-1.523\) （リアプノフ）、\(-1.027\) （ゼロ行列）、\(-8.015\) （拡散事前分布）とばらつきますが、\(t=2\) 以降は3通りともほぼ同じ値に収束します（カルマンフィルタが観測を重ねるにつれて誤った初期値の影響を「忘れていく」ため）。サンプルサイズが小さいほど、この初期化誤差が尤度・パラメータ推定に無視できない影響を残すため、定常ARMAでは式(5)のリアプノフ方程式による厳密な初期化が推奨されます。
なぜこの理解が実務で役立つか
- 欠測値への対応：カルマンフィルタは観測が欠測した時刻でも更新ステップをスキップするだけで自然に扱えます。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など）。
- 多変量への一般化：本記事は単変量ARMAを扱いましたが、同じ状態空間表現とカルマンフィルタによる最尤推定は多変量のVARMAモデルにもそのまま拡張されます。Düker et al. (2024) のVARMAモデルに関するレビュー論文は、識別性・推定・診断の各手法を整理する章の中で「状態空間表現とカルマンフィルタ」を主要な推定手法の一つとして位置づけており、本記事の議論が単変量の特殊ケースに留まらないことを裏付けています。
関連記事
- 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.
- Düker, M.-C., Matteson, D. S., Tsay, R. S., & Wilms, I. (2024). Vector Autoregressive Moving Average Models: A Review. arXiv:2406.19702.