ARIMA・SARIMAの理論とPython実装：statsmodels.tsa.arima.model/ACF/PACF/自己相関による時系列予測

はじめに

時系列データの予測は、需要予測、株価分析、気象予報など幅広い分野で必要とされます。**ARIMA（AutoRegressive Integrated Moving Average）**モデルは、時系列予測の古典的かつ強力な手法です。

本記事は、AR・MA・ARIMA・SARIMA の数理を天下り的にではなく単位円上の根から導出し、Box-Jenkins法（同定 → 推定 → 診断的検証）の枠組みに沿ってPython実装を通しで解説する「基礎理論と実践」編です。カルマンフィルタによる正確な最尤推定の理論的背景はhttps://yuhi-sa.github.io/posts/20260715_arima_state_space_kalman/1/に、PACFによる次数同定の詳細な導出はhttps://yuhi-sa.github.io/posts/20260710_pacf_ar_identification/1/にそれぞれ譲り、本記事ではBox-Jenkins法全体の実践的な流れ（定常性検定→次数同定→条件付き最尤推定→残差診断→予測）を一つの具体例で最初から最後まで動かすことに焦点を絞ります。

AR（自己回帰）モデル

AR(\(p\) )モデルは、過去 \(p\) 時点の値を用いて現在の値を表現します。

ここで \(\phi_i\) は自己回帰係数、\(\varepsilon_t \sim \mathcal{N}(0, \sigma^2)\) はホワイトノイズです。

定常性条件の導出

AR(1)モデル \(X_t = \phi X_{t-1} + \varepsilon_t\) を例に、なぜ「係数の大きさ」が定常性を決めるのかを導出します。式を過去に向かって \(k\) 回反復代入すると

が得られます。\(k \to \infty\) の極限で、\(|\phi| < 1\) ならば \(\phi^k X_{t-k} \to 0\) （平均二乗収束）となり、

という過去の有限の分散を持つショックだけで書ける因果的（causal）なMA(∞)表現が得られます。このとき分散は等比級数の和として

と有限確定します。逆に \(|\phi| > 1\) なら式(2)の \(\phi^k X_{t-k}\) は発散し、\(|\phi| = 1\) （単位根、ランダムウォーク）では分散が時刻とともに際限なく増加し、どちらも時不変な有限分散を持つ定常過程にはなりません。

AR(1)の特性方程式は \(1 - \phi z = 0\) で、根は \(z = 1/\phi\) です。したがって

——根が単位円の外側にあることと定常性は同値です。この対応は一般の AR(\(p\) ) にも拡張できます。特性方程式

の根をすべて \(z_1, \ldots, z_p\) とすると、部分分数分解により \(1/\phi(z)\) は \(z\) のべき級数 \(\sum_j \psi_j z^j\) として展開でき、その収束半径はちょうど \(\min_i |z_i|\) になります。係数 \(\{\psi_j\}\) が絶対総和可能（\(\sum_j |\psi_j| < \infty\) 、すなわち因果的MA(∞)表現が存在する条件）であるためには、収束半径が \(1\) より大きくなければなりません。つまり

すべての根が単位円の外側にあることが、AR(\(p\) )が定常なMA(∞)表現を持つための必要十分条件です。

これを本記事で使う数値実験（後述）のAR(2)過程 \(\phi_1=0.5, \phi_2=-0.3\) で検算すると、特性方程式 \(1 - 0.5z + 0.3z^2 = 0\) の根は \(z = 0.8333 \pm 1.6245i\) 、絶対値は \(|z| = 1.8257 > 1\) となり、確かに定常であることが確認できます（後述のPythonコードで数値的に再現します）。

MA（移動平均）モデル

MA(\(q\) )モデルは、過去 \(q\) 時点のノイズの線形結合で表現します。

MA モデルは有限個の定常な項の和なので、\(\theta_j\) の値によらず常に定常です。

反転可能性（invertibility）の導出

定常性は自動的に満たされますが、パラメータを一意に識別してデータから推定するには別の条件——反転可能性——が必要です。MA(1) \(X_t = \varepsilon_t + \theta \varepsilon_{t-1}\) を考えます。バックシフト演算子 \(B\) （\(B\varepsilon_t = \varepsilon_{t-1}\) ）を使うと \(X_t = (1+\theta B)\varepsilon_t\) と書け、形式的に逆演算子を展開すると

というAR(∞)表現が得られます。この級数が収束する（現在のショック \(\varepsilon_t\) を過去の観測値だけから復元できる）条件は \(|\theta| < 1\) です。特性方程式 \(1 + \theta z = 0\) の根は \(z = -1/\theta\) なので、これも「根が単位円の外側」という AR の定常性条件と同じ形の条件になります。

反転可能性が必要な理由は、パラメータの非一意性にあります。MA(1)の理論ACFは

であり、\(\theta \to 1/\theta\) と置き換えても \(\rho_1\) の値は変わりません（例えば \(\theta=2\) と \(\theta=0.5\) は同じ自己相関構造を持つ）。ACFだけからは \(\theta\) と \(1/\theta\) を区別できないため、慣例として \(|\theta|<1\) を満たす方を採用し、一意な表現に固定します。一般の MA(\(q\) ) では、\(\theta(z) = 1 + \theta_1 z + \cdots + \theta_q z^q\) の根がすべて単位円の外側にあることが反転可能性の条件です。

ARMA モデル

AR と MA を組み合わせた ARMA(\(p, q\) ) モデルです。

定常性はAR部分の特性方程式の根で、反転可能性はMA部分の特性方程式の根でそれぞれ独立に判定します。両方を満たして初めて、パラメータが一意に識別可能な定常過程になります。

ARIMA モデル

非定常な時系列に対して、\(d\) 回の差分をとって定常化した後に ARMA モデルを適用するのが ARIMA(\(p, d, q\) ) です。

ここで \(\Delta X_t = X_t - X_{t-1}\) は1次差分です。

• \(p\) : 自己回帰の次数
• \(d\) : 差分の回数（通常0, 1, 2）
• \(q\) : 移動平均の次数

条件付き最尤法によるパラメータ推定

Box-Jenkins法は「同定（identification）→ 推定（estimation）→ 診断的検証（diagnostic checking）」の3段階を反復するフレームワークです。ここでは推定段階の核心である**条件付き最尤法（Conditional Sum of Squares, CSS）**を導出します。

AR(\(p\) )モデルで \(\varepsilon_t \sim \mathcal{N}(0, \sigma^2)\) を仮定すると、式(1)より

というマルコフ性が成り立ちます。最初の \(p\) 個の観測値 \(X_1, \ldots, X_p\) を所与として条件付けると、同時密度はこの条件付き分布の積に分解できます。

対数を取ると、条件付き対数尤度は

となります。\(\sigma^2\) を固定して \(\phi_i\) について式(7)を最大化することは、右辺第2項の残差平方和を最小化することと同値——つまり**\(x_t\) を \(x_{t-1}, \ldots, x_{t-p}\) に回帰する通常の最小二乗法（OLS）と完全に一致**します。これが「条件付き最小二乗法（CSS）」と呼ばれる所以です。\(\hat{\phi}\) を代入して \(\partial \ell_c/\partial \sigma^2 = 0\) を解くと、MLE推定量は残差平方和を自由度で割った

になります。ARMA一般の場合も同様の考え方が使えますが、MA項があると \(\varepsilon_{t-j}\) が観測されないため、通常は \(\varepsilon_0 = \varepsilon_{-1} = \cdots = 0\) と初期化するか、逆向きに予測（バックキャスト）してから条件付き尤度を組み立てます。

CSS法は計算が単純な反面、最初の \(p\) 個の観測値が持つ周辺分布の情報を捨てているという近似です。これに対し、\(X_1, \ldots, X_p\) の定常分布を初期条件として取り込み、カルマンフィルタの予測誤差分解によって観測列全体の正確な（近似なしの）ガウス対数尤度を計算する方法があり、statsmodelsのARIMAクラスはデフォルトでこちらを内部的に使っています。この正確な尤度の導出と、CSS法との数値的な一致・乖離の検証はhttps://yuhi-sa.github.io/posts/20260715_arima_state_space_kalman/1/で詳しく扱っています。本記事でも後述のPython実装で、CSS法を手動実装した結果とstatsmodelsの厳密推定を数値的に比較します。

パラメータの選択（同定ステップ）

ACF（自己相関関数）と PACF（偏自己相関関数）

なぜAR過程でPACFが切断するのか（中間ラグの線形効果を除去した相関であるため）、Yule-Walker方程式やLevinson-Durbin再帰による具体的な計算方法はhttps://yuhi-sa.github.io/posts/20260710_pacf_ar_identification/1/で数式と数値実験の両方から掘り下げています。本記事ではこの判定表を前提に、実際にPythonでプロットして次数を読み取るところから話を進めます。

情報量規準

AIC（赤池情報量規準）やBICを使って、複数の \((p, d, q)\) の組み合わせから最適なモデルを選択します。

ここで \(L\) は最大尤度、\(k\) はパラメータ数です。

Python実装：定常性検定からモデル推定まで

データ生成とADF検定

真のAR(2)過程（\(\phi_1=0.5, \phi_2=-0.3\) ）に線形トレンドを加えたデータを生成し、拡張Dickey-Fuller（ADF）検定で定常性を確認します。

import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller

# --- データ生成（AR(2)プロセス + トレンド） ---
np.random.seed(42)
n = 300
noise = np.random.normal(0, 1, n)
data = np.zeros(n)
for t in range(2, n):
    data[t] = 0.5 * data[t-1] - 0.3 * data[t-2] + noise[t]
data += np.linspace(0, 10, n)  # トレンドを追加

ts = pd.Series(data, index=pd.date_range('2024-01-01', periods=n, freq='D'))

# --- AR(2)部分の特性方程式の根を検算 ---
roots = np.roots([-(-0.3), -(0.5), 1])  # 1 - 0.5z + 0.3z^2 = 0
print("特性方程式の根:", roots, "絶対値:", np.abs(roots))

# --- 定常性の検定（ADF検定） ---
result = adfuller(ts)
print(f"ADF統計量: {result[0]:.4f}, p値: {result[1]:.4f}")
print("臨界値:", result[4])

ts_diff = ts.diff().dropna()
result_diff = adfuller(ts_diff)
print(f"1次差分後 - ADF統計量: {result_diff[0]:.4f}, p値: {result_diff[1]:.4f}")

実行結果:

特性方程式の根: [0.83333333+1.62446572j 0.83333333-1.62446572j] 絶対値: [1.82574186 1.82574186]
ADF統計量: -0.1789, p値: 0.9410
臨界値: {'1%': -3.4532, '5%': -2.8716, '10%': -2.5721}
1次差分後 - ADF統計量: -7.6752, p値: 0.0000
臨界値(diff): {'1%': -3.4538, '5%': -2.8718, '10%': -2.5723}

特性方程式の根の絶対値 \(1.8257\) は前節で導出した理論値と一致し、AR(2)部分は定常です。しかし元系列は線形トレンドを含むためADF統計量 \(-0.1789\) （臨界値をすべて上回り、\(p=0.9410\) ）で単位根の帰無仮説を棄却できず——見かけ上非定常と判定されます。1次差分後はADF統計量 \(-7.6752\) （1%臨界値 \(-3.4538\) を大きく下回り、\(p \approx 0\) ）で明確に定常と判定され、\(d=1\) の妥当性が数値的に確認できます。

ACF/PACFプロットによる次数同定

import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(ts_diff, ax=axes[0], lags=20, title='ACF (1st Difference)')
plot_pacf(ts_diff, ax=axes[1], lags=20, method='ywm', title='PACF (1st Difference)')
plt.tight_layout()
plt.show()

PACFがラグ2で切断し、判定表どおりAR(2)成分が示唆されます（差分後もAR構造が残る理由は次項の「過差分」で説明します）。

条件付きMLE(CSS)とexact MLEの数値比較

前節で導出した条件付き最小二乗法を手動実装し、statsmodelsの厳密推定（カルマンフィルタベース）と比較します。

train = ts[:250]
test = ts[250:]

# --- 条件付きMLE(CSS)を手動実装 ---
y = train.diff().dropna().values
X = np.column_stack([y[1:-1], y[0:-2]])  # ラグ1, ラグ2
yt = y[2:]
beta = np.linalg.lstsq(X, yt, rcond=None)[0]
resid_css = yt - X @ beta
sigma2_css = np.sum(resid_css**2) / len(resid_css)
print("[CSS法] phi1, phi2:", beta, "sigma2:", sigma2_css)

# --- statsmodelsの厳密推定(カルマンフィルタベースMLE) ---
fitted = ARIMA(train, order=(2, 1, 0)).fit()
print("[exact MLE] phi1, phi2:", fitted.arparams, "sigma2:", fitted.params.iloc[-1])
print(fitted.summary())

実行結果:

[CSS法] phi1, phi2: [-0.14256314 -0.3405605 ] sigma2: 1.3407 (n=247)
[exact MLE] phi1, phi2: [-0.14198861 -0.33822858] sigma2: 1.3316 (llf=-389.10, N=249)

                 coef    std err          z      P>|z|
ar.L1         -0.1420      0.067     -2.125      0.034
ar.L2         -0.3382      0.064     -5.307      0.000
sigma2         1.3316      0.120     11.059      0.000
AIC: 784.19   BIC: 794.74

CSS法（\(\hat{\phi}_1=-0.1426, \hat{\phi}_2=-0.3406\) ）とstatsmodelsの厳密MLE（\(\hat{\phi}_1=-0.1420, \hat{\phi}_2=-0.3382\) ）は小数点以下2〜3桁まで一致します。わずかな差は、CSS法が最初の2点を条件付けて捨てる（有効データ数 \(N=247\) ）のに対し、厳密MLEは定常分布を初期条件として全 \(N=249\) 点の情報を使うことに由来します——https://yuhi-sa.github.io/posts/20260715_arima_state_space_kalman/1/で導出したカルマンフィルタの予測誤差分解が、この差を埋める正確な尤度計算に対応します。

残差診断とモデル選択の落とし穴：過差分（overdifferencing）

Box-Jenkins法の3段階目「診断的検証」として、Ljung-Box検定で残差が白色雑音とみなせるかを確認します。

from statsmodels.stats.diagnostic import acorr_ljungbox

lb = acorr_ljungbox(fitted.resid, lags=[10, 20], return_df=True)
print(lb)

実行結果:

      lb_stat     lb_pvalue
10  48.436410  5.165530e-07
20  60.495144  5.968126e-06

両ラグとも \(p < 0.001\) で棄却され、残差に有意な自己相関が残っています。 ADF検定が1次差分後の定常性を支持し、PACFもラグ2での切断を示していたにもかかわらず、ARIMA(2,1,0)は誤指定だったことになります。

原因は**過差分（overdifferencing）**です。本記事のデータ生成過程は「定常なAR(2) + 決定論的な線形トレンド」であり、確率的な単位根（stochastic trend）を持ちません。にもかかわらず1次差分 \(d=1\) を取ると、定常AR(2)過程を差分したことで人為的に非可逆なMA成分が混入します（AR(1)の例で言えば、\(x_t - x_{t-1} = \phi(x_{t-1}-x_{t-2}) + (\varepsilon_t - \varepsilon_{t-1})\) という形になり、右辺の \(\varepsilon_t - \varepsilon_{t-1}\) 項がMA(1)成分を生む）。この人為的なMA成分をARIMA(2,1,0)（MA次数0）でモデル化しようとした結果、残差に構造が残ったのです。

これを2通りの方法で修正し、AIC・残差診断・予測誤差を比較します。

# 修正案1: MA項を追加してMA成分を吸収
m211 = ARIMA(train, order=(2, 1, 1)).fit()
lb211 = acorr_ljungbox(m211.resid, lags=[10, 20], return_df=True)

# 修正案2: 差分せず、決定論的トレンドを明示的にモデル化（真のDGPと一致）
m200t = ARIMA(train, order=(2, 0, 0), trend='t').fit()
lb200t = acorr_ljungbox(m200t.resid, lags=[10, 20], return_df=True)

for name, m, lb in [("ARIMA(2,1,0)", fitted, lb), ("ARIMA(2,1,1)", m211, lb211), ("ARIMA(2,0,0)+trend", m200t, lb200t)]:
    fc = m.forecast(steps=len(test))
    rmse = np.sqrt(np.mean((test.values - fc.values)**2))
    mae = np.mean(np.abs(test.values - fc.values))
    print(f"{name}: AIC={m.aic:.2f}, LB(10) p={lb['lb_pvalue'].iloc[0]:.3f}, RMSE={rmse:.4f}, MAE={mae:.4f}")

実行結果:

真のデータ生成過程と一致する ARIMA(2,0,0)+trend が AIC・残差診断・予測誤差のすべてで最良という結果になりました。一方、AIC単体で見れば ARIMA(2,1,0) は ARIMA(2,1,1) より劣っているのに対し、予測RMSEだけを見るとARIMA(2,1,0)（1.2494）はARIMA(2,1,1)（1.5652）より良く見えてしまいます。これは特定のテスト窓での偶然の当たり外れであり、残差が白色雑音でない時点で予測区間の信頼性は保証されません。AICや点予測のRMSEだけで判断せず、必ず残差診断（Ljung-Box検定）とセットで評価することがBox-Jenkins法の要諦です。

なお、単純な「直近値を維持する」ナイーブ予測のRMSEは \(1.4017\) 、MAEは \(1.1301\) でした。ARIMA(2,1,0)はこのナイーブ基準は上回っていますが、真の仕様である ARIMA(2,0,0)+trend にはさらに劣ります。

過差分ARIMA(2,1,0)（赤破線）はテスト期間のトレンドの伸びに追従できずほぼ横ばいの予測になっているのに対し、正しい仕様のARIMA(2,0,0)+trend（青破線）はトレンドに沿った予測ができていることが視覚的にも確認できます。

SARIMA（季節性ARIMA）

季節性を持つデータには SARIMA(\(p, d, q\) )(\(P, D, Q\) )\(_s\) を使います。

ここで \(s\) は季節周期（月次データなら12）、\(B\) は後退オペレータ（\(BX_t = X_{t-1}\) ）です。

数値実験：季節性を無視するとどれだけ悪化するか

真のSARIMA(1,1,1)(1,1,1,12)過程（SARIMAX.simulateで生成、\(\phi=0.5, \theta=-0.4, \Phi=0.3, \Theta=-0.5\) ）から月次データ180点を生成し、季節差分後のADF検定、SARIMAXでの推定、季節性を無視した非季節ARIMA(1,1,1)との予測比較を行います。

from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller
from statsmodels.stats.diagnostic import acorr_ljungbox

n = 180
dummy = pd.Series(np.zeros(n), index=pd.date_range('2010-01-01', periods=n, freq='MS'))
mod = SARIMAX(dummy, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12), trend='n')
true_params = [0.5, -0.4, 0.3, -0.5, 1.0]  # ar, ma, ar.S, ma.S, sigma2
sim = mod.simulate(true_params, n, initial_state=np.zeros(mod.k_states),
                    random_state=9) + 50  # random_state指定で再現性を確保
ts_s = pd.Series(sim.values, index=dummy.index)
train_s, test_s = ts_s[:150], ts_s[150:]

# 季節差分後のADF検定
res_sd = adfuller(ts_s.diff(12).dropna())
print(f"季節差分後ADF: {res_sd[0]:.4f}, p値: {res_sd[1]:.4f}")

# SARIMAモデル
sarima_fit = SARIMAX(train_s, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12)).fit(disp=False)
lb_s = acorr_ljungbox(sarima_fit.resid.iloc[13:], lags=[12, 24], return_df=True)
fc_s = sarima_fit.get_forecast(steps=len(test_s)).predicted_mean
rmse_s = np.sqrt(np.mean((test_s.values - fc_s.values)**2))
mae_s = np.mean(np.abs(test_s.values - fc_s.values))

# 季節性を無視した非季節ARIMA
nonseason_fit = ARIMA(train_s, order=(1, 1, 1)).fit()
fc_ns = nonseason_fit.forecast(steps=len(test_s))
rmse_ns = np.sqrt(np.mean((test_s.values - fc_ns.values)**2))
mae_ns = np.mean(np.abs(test_s.values - fc_ns.values))

print(f"SARIMA: AIC={sarima_fit.aic:.2f}, RMSE={rmse_s:.3f}, MAE={mae_s:.3f}")
print(f"非季節ARIMA: AIC={nonseason_fit.aic:.2f}, RMSE={rmse_ns:.3f}, MAE={mae_ns:.3f}")
print(lb_s)

実行結果:

季節差分後ADF: -3.4639, p値: 0.0090

SARIMA(1,1,1)(1,1,1,12): AIC=384.24, RMSE=3.872, MAE=3.595
      lb_stat  lb_pvalue
12   8.024676   0.783199
24  14.797854   0.926534

非季節ARIMA(1,1,1): AIC=685.03, RMSE=7.720, MAE=6.421

季節差分後のADF統計量は \(-3.4639\) （\(p=0.0090\) ）で定常性が確認され、SARIMAの残差はLjung-Box検定で白色雑音と判定されます（\(p=0.783, 0.927\) ）。季節性を無視した非季節ARIMA(1,1,1)は、AICが \(685.03\) とSARIMAの \(384.24\) より大幅に悪化し、RMSE/MAEもおよそ2.0倍/1.8倍に悪化しました（RMSE: \(3.872 \to 7.720\) 、MAE: \(3.595 \to 6.421\) ）。季節周期を明示的にモデルへ組み込むことの効果が数値で裏付けられます。

なお、SARIMAX.simulate は乱数生成に独自の内部状態を持ち、np.random.seed() によるグローバルなシード指定だけでは結果が再現されません。上記コードのように random_state 引数へ明示的に整数を渡すことで、初めて実行のたびに同じ系列が得られます（これも季節性モデルを扱う際の見落としやすい落とし穴です）。

モデル選択の自動化（auto_arima）とその限界

pmdarima の auto_arima で最適なパラメータを自動探索できますが、探索アルゴリズムの性質上、常に最良のモデルが選ばれるとは限りません。

import pmdarima as pm

# デフォルト: 段階的探索（stepwise Hyndman-Khandakar法）
auto_stepwise = pm.auto_arima(train, seasonal=False, stepwise=True, suppress_warnings=True)
print("stepwise選択:", auto_stepwise.order, "AIC:", auto_stepwise.aic())

# 全探索（stepwise=False）
auto_full = pm.auto_arima(train, seasonal=False, stepwise=False, suppress_warnings=True,
                           max_p=5, max_q=5, max_d=2, n_jobs=-1)
print("全探索選択:", auto_full.order, "AIC:", auto_full.aic())

実行結果:

stepwise選択: (0, 1, 0) AIC: 813.02
全探索選択: (5, 1, 0) AIC: 732.96

デフォルトの段階的探索は ARIMA(0,1,0)（AIC \(813.02\) ）というほぼ何もしないモデルを選び、前節で手動比較した ARIMA(2,1,0)（AIC \(784.19\) ）よりも悪化しています。stepwise=False で全探索に切り替えると ARIMA(5,1,0)（AIC \(732.96\) ）とやや改善しますが、それでも真の仕様である ARIMA(2,0,0)+trend（AIC \(698.94\) ）には及びません。これは2つの独立した問題を示しています。

• 段階的探索は局所最適に陥りうる: stepwise=True は候補空間全体を探索しないため、今回のように最適から大きく外れることがあります。
• 探索空間自体が真のモデルを含んでいない: auto_arima(seasonal=False) は差分次数 \(d\) とAR/MA次数だけを探索し、決定論的トレンド項を候補に含めません。今回のようにトレンドが確率的単位根ではなく決定論的な場合、AR/MA次数をどれだけ増やしても真の仕様には到達できず、次数を無駄に上げること（\(p=5\) ）で近似するしかありません。

自動化ツールはあくまで候補空間内でのヒューリスティック探索であり、データ生成過程についてのドメイン知識（トレンドが決定論的か確率的か、季節周期の有無など）を踏まえた探索空間の設計が依然として必要です。

実務上の注意点

• 過差分に注意する: ADF検定・PACFの判定表だけに頼ると、決定論的トレンドを持つ定常過程を誤って差分してしまうことがあります。差分後は必ずLjung-Box検定で残差の白色性を確認してください。
• 反転可能性の境界に注意する: MA係数の推定値が \(\pm 1\) 付近（非可逆の境界）になる場合、標準誤差が極端に大きくなることがあります（本記事の予備実験でも、季節差分・非季節差分を両方適用した過差分ケースで ma.L1 の標準誤差が数百倍に膨張する現象を確認しました）。次数や差分回数を見直すサインです。
• AICと点予測のRMSEは別の指標: 本記事の数値実験でも、残差診断で棄却されたモデルの方が特定のテスト窓ではRMSEが良く見える逆転が起きました。AIC・残差診断・予測誤差は三位一体で確認します。
• 自動化ツールは万能ではない: auto_arima の段階的探索は局所最適に陥ることがあり、探索空間（トレンド項の有無、季節性の有無）の設計はユーザーの責任です。

最新研究動向

古典的なBox-Jenkins法（ARIMA）と機械学習手法の実務上の使い分けは、現在も活発に研究されています。Tjøstheim (2025) の Selected Topics in Time Series Forecasting: Statistical Models vs. Machine Learning（Entropy 27(3), 279, DOI: 10.3390/e27030279）は、Makridakis予測コンペティション（M1〜M6）の結果を横断的にレビューし、ARIMAや指数平滑法などの古典的な線形パラメトリックモデルが、低頻度（月次など）でスペクトルエントロピーの低い「素直な」系列では依然として機械学習に匹敵する精度を示す一方、ボラティリティ予測や高頻度金融データのような非線形性・分散変動が支配的な領域では機械学習が優位に立つことを整理しています。M3では指数平滑法の拡張であるThetaメソッドが優勝し、M4では指数平滑法と深層学習を組み合わせたハイブリッド手法が優勝するなど、「古典的統計モデル vs 機械学習」ではなく「両者の組み合わせ」が近年の潮流であるという指摘は、本記事のARIMAが単体の万能モデルではなく、データの性質（線形性・エントロピー）を見極めて使い分けるべき手法の一つであることを裏付けています。

Python実装まとめ

# 上記の全コードをまとめた最終モデル（真の仕様に近いARIMA）
final_model = ARIMA(train, order=(2, 0, 0), trend='t').fit()
print(final_model.summary())
forecast = final_model.get_forecast(steps=len(test))
mean_forecast = forecast.predicted_mean
conf_int = forecast.conf_int()

おすすめ書籍

AR・MA・ARMA モデルから単位根・見せかけの回帰まで、時系列分析の標準理論を丁寧に解説した定番教科書です。

※ 上記は Amazon アソシエイトのリンクです。

関連記事

• ARMA/ARIMAの状態空間表現とカルマンフィルタによる最尤推定 - statsmodelsのARIMA推定が内部で使っているカルマンフィルタベースの正確な尤度計算を導出し、本記事の条件付きMLE(CSS)との数値的な一致・乖離を検証した発展編です。
• PACF・ARモデル次数同定 - 本記事の次数決定ステップを深掘りし、Yule-Walker方程式・Levinson-Durbin再帰によるPACFの数理と、AR次数 \(p\) の同定を解説しています。
• GARCH(1,1)モデルの理論とPython実装 - 本記事が捉えられない分散のクラスタリング（ボラティリティクラスタリング）を補完するモデルです。
• 機械学習による時系列予測・分類・異常検知の選び方ハブ - ARIMA を含む 7 手法（k-means/GMM/RF/GBDT/LSTM/Kalman/IsolationForest）を比較し、データ特性に応じた選択指針を提供します。
• LSTMによる時系列予測 - ARIMA の線形仮定に対し、LSTM は非線形・長期依存を扱えます。両者の使い分けが整理できます。
• デジタル信号処理と機械学習の学習ロードマップ - ARIMA を含む 5 大ハブの学習動線を体系化したメタロードマップです。
• 時系列データの異常検知：統計的手法からカルマンフィルタまで - ARIMAによる予測値からの逸脱を異常検知に利用する手法を解説しています。
• 移動平均フィルタの種類と比較 - 移動平均フィルタの種類を解説しており、MAモデルの直感的な理解に役立ちます。
• カルマンフィルタの理論とPython実装 - 状態空間モデルに基づく時系列推定手法であり、ARIMAとの比較が有益です。
• Matplotlib実践Tips：論文品質のグラフを作る - 時系列グラフの可視化品質を高めるTipsを紹介しています。
• 指数移動平均（EMA）フィルタの周波数特性 - EMAは時系列の平滑化手法であり、ARIMAのMA成分と関連します。
• ベイズ線形回帰の基礎 - 時系列の回帰的な側面をベイズ的に捉える手法を解説しています。

参考文献

• Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
• Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and Practice (3rd ed.).
• Tjøstheim, D. (2025). “Selected Topics in Time Series Forecasting: Statistical Models vs. Machine Learning.” Entropy, 27(3), 279. DOI: 10.3390/e27030279
• statsmodels ARIMA documentation