時系列データの分析・予測において、ARIMA（AutoRegressive Integrated Moving Average）モデルは最も基本的かつ広く使われている手法の一つです。株価、気温、売上高など、時間とともに変化するデータのパターンを捉え、将来の値を予測するための数学的な枠組みを提供します。
ARIMA は AR（自己回帰）、I（差分）、MA（移動平均）の3つの構成要素を組み合わせたモデルであり、Box-Jenkins 法と呼ばれる体系的なモデル構築手法によって運用されます。本記事では、各構成要素の数学的定式化を省略なく導出し、Python の statsmodels を用いた実践的な実装を行います。
本記事の内容
- AR(p) モデルの定式化と定常性条件（特性方程式の根が単位円外）
- MA(q) モデルの定式化と可逆性条件
- ARMA(p,q) モデル
- ARIMA(p,d,q) モデル（非定常時系列への拡張）
- Box-Jenkins 法（同定→推定→診断のサイクル）
- ACF / PACF による次数決定
- AIC / BIC によるモデル選択
- SARIMA（季節性 ARIMA の拡張）
- Python での実装（ADF 検定→ACF/PACF→モデル選択→予測→残差診断）
前提知識
この記事を読む前に、以下の概念を理解しておくと理解が深まります。
- 確率変数、期待値、分散、共分散の基本概念
- 回帰分析の基礎
時系列データと定常性
時系列データの定義
時系列データとは、時間の順序に沿って観測されたデータの列 $\{x_t\}_{t=1}^{T}$ のことです。ここで $x_t$ は時刻 $t$ における観測値であり、確率変数と見なします。
弱定常性（Wide-Sense Stationarity）
時系列 $\{x_t\}$ が弱定常であるとは、以下の条件を満たすことです。
- 期待値が一定: $\mathbb{E}[x_t] = \mu$ （全ての $t$ に対して同じ値）
- 分散が有限かつ一定: $\text{Var}[x_t] = \gamma(0) < \infty$
- 自己共分散が時間差のみに依存: $\text{Cov}[x_t, x_{t+h}] = \gamma(h)$ （$t$ に依存しない）
ここで $\gamma(h)$ は自己共分散関数（autocovariance function）と呼ばれます。これを正規化した自己相関関数（ACF: Autocorrelation Function）は
$$ \rho(h) = \frac{\gamma(h)}{\gamma(0)} = \frac{\text{Cov}[x_t, x_{t+h}]}{\text{Var}[x_t]} $$
で定義されます。$\rho(0) = 1$, $|\rho(h)| \leq 1$ です。
ホワイトノイズ
ホワイトノイズ $\{\varepsilon_t\}$ は最も基本的な定常過程で、以下を満たします。
$$ \mathbb{E}[\varepsilon_t] = 0, \quad \text{Var}[\varepsilon_t] = \sigma^2, \quad \text{Cov}[\varepsilon_t, \varepsilon_s] = 0 \quad (t \neq s) $$
つまり、各時点の値は無相関で、同一の分散を持ちます。ACF は $\rho(0) = 1$, $\rho(h) = 0$ ($h \neq 0$) です。
AR(p) モデル（自己回帰モデル）
定式化
AR(p) モデルは、現在の値が過去 $p$ 個の値の線形結合とホワイトノイズで表されるモデルです。
$$ x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_p x_{t-p} + \varepsilon_t $$
$$ x_t = \sum_{i=1}^{p} \phi_i x_{t-i} + \varepsilon_t $$
ここで $\phi_1, \phi_2, \dots, \phi_p$ は AR 係数、$\varepsilon_t \sim \text{WN}(0, \sigma^2)$ はホワイトノイズです。
後退シフト演算子（backshift operator） $B$ を用いると簡潔に書けます。$Bx_t = x_{t-1}$, $B^k x_t = x_{t-k}$ として
$$ (1 – \phi_1 B – \phi_2 B^2 – \cdots – \phi_p B^p) x_t = \varepsilon_t $$
$$ \phi(B) x_t = \varepsilon_t $$
ここで $\phi(B) = 1 – \phi_1 B – \phi_2 B^2 – \cdots – \phi_p B^p$ はAR 多項式です。
AR(1) の具体例
最も単純な AR(1) モデルは
$$ x_t = \phi x_{t-1} + \varepsilon_t $$
再帰的に展開すると
$$ \begin{align} x_t &= \phi x_{t-1} + \varepsilon_t \\ &= \phi(\phi x_{t-2} + \varepsilon_{t-1}) + \varepsilon_t \\ &= \phi^2 x_{t-2} + \phi \varepsilon_{t-1} + \varepsilon_t \\ &= \cdots \\ &= \phi^n x_{t-n} + \sum_{j=0}^{n-1} \phi^j \varepsilon_{t-j} \end{align} $$
$|\phi| < 1$ のとき、$n \to \infty$ で $\phi^n x_{t-n} \to 0$ となり
$$ x_t = \sum_{j=0}^{\infty} \phi^j \varepsilon_{t-j} $$
これは因果的な MA($\infty$) 表現であり、過去のホワイトノイズの加重和として $x_t$ が表現できることを意味します。
定常性条件
AR(p) モデルが定常であるための条件を導出します。
特性方程式: AR 多項式 $\phi(z) = 1 – \phi_1 z – \phi_2 z^2 – \cdots – \phi_p z^p = 0$ の根を $z_1, z_2, \dots, z_p$ とします。
定常性条件: AR(p) モデルが弱定常であるためには、特性方程式の全ての根が単位円の外側にある必要があります。
$$ |z_i| > 1 \quad \text{for all } i = 1, 2, \dots, p $$
AR(1) の場合の導出: $\phi(z) = 1 – \phi z = 0$ より $z = 1/\phi$。$|z| > 1$ は $|1/\phi| > 1$、すなわち $|\phi| < 1$ と同値です。
AR(2) の場合の導出: $\phi(z) = 1 – \phi_1 z – \phi_2 z^2 = 0$ の2根 $z_1, z_2$ が $|z_i| > 1$ を満たす条件は
$$ \phi_1 + \phi_2 < 1, \quad \phi_2 - \phi_1 < 1, \quad |\phi_2| < 1 $$
の3つの不等式で与えられます。これは $(1, 0)$, $(-1, 0)$, $(0, -1)$ を頂点とする三角形の内部に $(\phi_1, \phi_2)$ があることに対応します。
AR(p) の ACF と PACF
AR(1) の自己相関関数は
$$ \rho(h) = \phi^{|h|} $$
で、$|\phi| < 1$ のとき指数的に減衰します。
AR(p) の ACF は一般に指数的・減衰振動的に減衰しますが、具体的な形は $\phi_1, \dots, \phi_p$ に依存します。
偏自己相関関数（PACF: Partial Autocorrelation Function）は、AR モデルの次数同定において重要な役割を果たします。$\alpha(h)$ は $x_t$ と $x_{t-h}$ の間の、$x_{t-1}, \dots, x_{t-h+1}$ の影響を除いた相関です。
AR(p) モデルの PACF は
$$ \alpha(h) = \begin{cases} \neq 0 & h \leq p \\ 0 & h > p \end{cases} $$
つまり、PACF がラグ $p$ で切断（cutoff）されることが AR(p) の特徴です。
MA(q) モデル（移動平均モデル）
定式化
MA(q) モデルは、現在の値が現在と過去 $q$ 個のホワイトノイズの線形結合で表されるモデルです。
$$ x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \cdots + \theta_q \varepsilon_{t-q} $$
$$ x_t = \sum_{j=0}^{q} \theta_j \varepsilon_{t-j} \quad (\theta_0 = 1) $$
後退シフト演算子を用いると
$$ x_t = \theta(B) \varepsilon_t $$
ここで $\theta(B) = 1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q$ はMA 多項式です。
MA(q) の自己共分散
MA(q) の自己共分散を導出します。
$$ \begin{align} \gamma(h) &= \text{Cov}[x_t, x_{t+h}] \\ &= \text{Cov}\left[\sum_{i=0}^{q}\theta_i\varepsilon_{t-i},\;\sum_{j=0}^{q}\theta_j\varepsilon_{t+h-j}\right] \\ &= \sum_{i=0}^{q}\sum_{j=0}^{q}\theta_i\theta_j\,\text{Cov}[\varepsilon_{t-i}, \varepsilon_{t+h-j}] \end{align} $$
$\text{Cov}[\varepsilon_{t-i}, \varepsilon_{t+h-j}] = \sigma^2$ only if $t – i = t + h – j$、すなわち $j = i + h$ のときのみ非零です。したがって
$$ \gamma(h) = \begin{cases} \sigma^2 \sum_{i=0}^{q-h} \theta_i \theta_{i+h} & 0 \leq h \leq q \\ 0 & h > q \end{cases} $$
つまり、MA(q) の ACF はラグ $q$ で切断される（$\rho(h) = 0$ for $h > q$）ことが特徴です。
MA(1) の具体例
MA(1): $x_t = \varepsilon_t + \theta \varepsilon_{t-1}$
$$ \gamma(0) = \sigma^2(1 + \theta^2), \quad \gamma(1) = \sigma^2 \theta, \quad \gamma(h) = 0 \; (h \geq 2) $$
$$ \rho(1) = \frac{\theta}{1 + \theta^2}, \quad \rho(h) = 0 \; (h \geq 2) $$
可逆性条件
MA(q) モデルの可逆性とは、MA 表現を AR($\infty$) 表現に変換できることです。
MA(1) の場合: $x_t = \varepsilon_t + \theta \varepsilon_{t-1}$ から $\varepsilon_t = x_t – \theta \varepsilon_{t-1}$ を再帰的に展開すると
$$ \varepsilon_t = x_t – \theta x_{t-1} + \theta^2 x_{t-2} – \theta^3 x_{t-3} + \cdots = \sum_{j=0}^{\infty} (-\theta)^j x_{t-j} $$
これが収束するためには $|\theta| < 1$ が必要です。
一般に、MA(q) モデルが可逆であるための条件は、MA 多項式 $\theta(z) = 0$ の全ての根が単位円の外側にあることです。
$$ |z_i| > 1 \quad \text{for all roots } z_i \text{ of } \theta(z) = 0 $$
ARMA(p,q) モデル
定式化
AR と MA を組み合わせた 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} $$
$$ \phi(B) x_t = \theta(B) \varepsilon_t $$
定常性条件: $\phi(z) = 0$ の全根が $|z| > 1$ 可逆性条件: $\theta(z) = 0$ の全根が $|z| > 1$
ARMA の ACF / PACF
| モデル | ACF | PACF |
|---|---|---|
| AR(p) | 指数/振動減衰 | ラグ $p$ で切断 |
| MA(q) | ラグ $q$ で切断 | 指数/振動減衰 |
| ARMA(p,q) | ラグ $q-p$ 以降で指数/振動減衰 | ラグ $p-q$ 以降で指数/振動減衰 |
この特徴が、ACF/PACF プロットからモデルの次数を推定する際の手がかりとなります。
ARIMA(p,d,q) モデル
非定常時系列と差分
現実のデータの多くは非定常です（例: トレンドを持つ時系列）。非定常時系列を定常にする最も一般的な方法は差分を取ることです。
1階差分: $\nabla x_t = x_t – x_{t-1} = (1 – B) x_t$
$d$ 階差分: $\nabla^d x_t = (1 – B)^d x_t$
$d$ 階差分を取って定常になった時系列に ARMA(p,q) を当てはめるモデルが ARIMA(p,d,q) です。
$$ \phi(B)(1 – B)^d x_t = \theta(B) \varepsilon_t $$
$y_t = \nabla^d x_t = (1 – B)^d x_t$ とおくと
$$ \phi(B) y_t = \theta(B) \varepsilon_t $$
すなわち、$y_t$ は ARMA(p,q) に従います。
ARIMA のパラメータの意味
| パラメータ | 意味 | 決定方法 |
|---|---|---|
| $p$ | AR の次数（過去何個の値に依存するか） | PACF プロット |
| $d$ | 差分の階数（何回差分すれば定常になるか） | ADF 検定 |
| $q$ | MA の次数（過去何個のノイズに依存するか） | ACF プロット |
Box-Jenkins 法
Box と Jenkins が体系化したモデル構築の方法論は、以下の3ステップで構成されます。
ステップ1: 同定（Identification）
- 時系列プロットを確認し、トレンド・季節性の有無を判断
- ADF（Augmented Dickey-Fuller）検定で定常性を検定
- 必要に応じて差分 $d$ を決定（$d = 0, 1, 2$ のいずれか）
- 差分後の ACF / PACF プロットから $p$, $q$ の候補を決定
ステップ2: 推定（Estimation）
候補モデルのパラメータ $\{\phi_i, \theta_j, \sigma^2\}$ を最尤推定法で推定します。
ARMA(p,q) の対数尤度は、データ $\bm{x} = (x_1, \dots, x_T)^T$ に対して
$$ \ell(\bm{\phi}, \bm{\theta}, \sigma^2) = -\frac{T}{2}\log(2\pi\sigma^2) – \frac{1}{2\sigma^2}\sum_{t=1}^{T}\varepsilon_t^2 $$
ここで $\varepsilon_t$ は $x_t$ から予測値を引いた残差です。パラメータの推定には数値最適化（条件付き最尤法、厳密最尤法等）が用いられます。
ステップ3: 診断（Diagnostic Checking）
推定したモデルの妥当性を以下の方法で検証します。
- 残差の ACF: 残差がホワイトノイズに近いか（有意な自己相関がないか）
- Ljung-Box 検定: 残差の自己相関の一括検定
$$ Q(m) = T(T+2)\sum_{h=1}^{m}\frac{\hat{\rho}_\varepsilon(h)^2}{T-h} $$
帰無仮説「残差に自己相関がない」のもとで $Q(m) \sim \chi^2(m – p – q)$ に従います。
- 残差の正規性: Q-Qプロットやシャピロ-ウィルク検定
AIC / BIC によるモデル選択
複数の候補モデルを比較する際、情報量基準を使います。
AIC（Akaike Information Criterion）:
$$ \text{AIC} = -2\ell(\hat{\bm{\theta}}) + 2k $$
BIC（Bayesian Information Criterion）:
$$ \text{BIC} = -2\ell(\hat{\bm{\theta}}) + k \log T $$
ここで $\ell(\hat{\bm{\theta}})$ は最大対数尤度、$k$ はパラメータ数（$k = p + q + 1$）、$T$ はサンプルサイズです。値が小さいほど良いモデルです。
BIC は AIC よりもパラメータ数に対するペナルティが大きく（$T > 8$ で $\log T > 2$）、よりパーシモニアス（倹約的）なモデルを選びやすい傾向があります。
SARIMA（季節性 ARIMA）
定式化
季節性を持つ時系列に対応するため、ARIMA を拡張した SARIMA モデルは ARIMA$(p,d,q) \times (P,D,Q)_s$ と表記されます。
$$ \Phi(B^s)\phi(B)(1-B)^d(1-B^s)^D x_t = \Theta(B^s)\theta(B)\varepsilon_t $$
ここで
- $\phi(B)$: 非季節 AR 多項式（次数 $p$）
- $\theta(B)$: 非季節 MA 多項式（次数 $q$）
- $\Phi(B^s)$: 季節 AR 多項式（次数 $P$）
- $\Theta(B^s)$: 季節 MA 多項式（次数 $Q$）
- $(1-B)^d$: 非季節差分
- $(1-B^s)^D$: 季節差分
- $s$: 季節周期（月次データなら $s = 12$）
季節差分
季節差分 $(1 – B^s)x_t = x_t – x_{t-s}$ は、$s$ 期前との差を取ります。例えば月次データで $s = 12$ の場合、$x_t – x_{t-12}$ は前年同月比に相当します。
Python での実装
データの準備と可視化
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller, acf, pacf
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings
warnings.filterwarnings('ignore')
# AirPassengers データセット（月次航空旅客数）
# データを手動生成（典型的な AirPassengers データ）
np.random.seed(42)
months = pd.date_range('1949-01', periods=144, freq='MS')
# トレンド + 季節性 + ノイズのシミュレーション
t = np.arange(144)
trend = 100 + 2.5 * t
seasonal = 30 * np.sin(2 * np.pi * t / 12) + 15 * np.cos(2 * np.pi * t / 12)
noise = np.random.randn(144) * 10
passengers = trend + seasonal + noise
passengers = np.maximum(passengers, 50) # 負にならないように
ts = pd.Series(passengers, index=months, name='Passengers')
# 時系列プロットの可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 原系列
axes[0, 0].plot(ts, 'b-', linewidth=1)
axes[0, 0].set_title('Original Time Series', fontsize=13)
axes[0, 0].set_xlabel('Date')
axes[0, 0].set_ylabel('Passengers')
axes[0, 0].grid(True, alpha=0.3)
# (2) 1階差分
ts_diff1 = ts.diff().dropna()
axes[0, 1].plot(ts_diff1, 'r-', linewidth=1)
axes[0, 1].set_title('First Difference', fontsize=13)
axes[0, 1].set_xlabel('Date')
axes[0, 1].set_ylabel('Diff(Passengers)')
axes[0, 1].grid(True, alpha=0.3)
# (3) 季節差分（12ヶ月）
ts_seasonal_diff = ts.diff(12).dropna()
axes[1, 0].plot(ts_seasonal_diff, 'g-', linewidth=1)
axes[1, 0].set_title('Seasonal Difference (lag=12)', fontsize=13)
axes[1, 0].set_xlabel('Date')
axes[1, 0].set_ylabel('Seasonal Diff')
axes[1, 0].grid(True, alpha=0.3)
# (4) 1階差分 + 季節差分
ts_both_diff = ts.diff().diff(12).dropna()
axes[1, 1].plot(ts_both_diff, 'm-', linewidth=1)
axes[1, 1].set_title('First + Seasonal Difference', fontsize=13)
axes[1, 1].set_xlabel('Date')
axes[1, 1].set_ylabel('Double Diff')
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
ADF 検定による定常性の確認
import numpy as np
import pandas as pd
from statsmodels.tsa.stattools import adfuller
# データの再生成
np.random.seed(42)
months = pd.date_range('1949-01', periods=144, freq='MS')
t = np.arange(144)
trend = 100 + 2.5 * t
seasonal = 30 * np.sin(2 * np.pi * t / 12) + 15 * np.cos(2 * np.pi * t / 12)
noise = np.random.randn(144) * 10
passengers = trend + seasonal + noise
passengers = np.maximum(passengers, 50)
ts = pd.Series(passengers, index=months, name='Passengers')
def adf_test(series, title=''):
"""ADF検定を実行し結果を表示する"""
result = adfuller(series, autolag='AIC')
print(f'--- ADF検定: {title} ---')
print(f'ADF統計量: {result[0]:.4f}')
print(f'p値: {result[1]:.4f}')
print(f'ラグ数: {result[2]}')
print(f'観測数: {result[3]}')
for key, val in result[4].items():
print(f'臨界値 ({key}): {val:.4f}')
if result[1] < 0.05:
print('→ 帰無仮説（単位根あり）を棄却。定常と判断。')
else:
print('→ 帰無仮説を棄却できない。非定常の可能性。')
print()
# 原系列の ADF 検定
adf_test(ts, '原系列')
# 1階差分の ADF 検定
adf_test(ts.diff().dropna(), '1階差分')
# 季節差分 + 1階差分
adf_test(ts.diff().diff(12).dropna(), '1階差分 + 季節差分')
ACF / PACF プロットによる次数の同定
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# データの再生成
np.random.seed(42)
months = pd.date_range('1949-01', periods=144, freq='MS')
t = np.arange(144)
trend = 100 + 2.5 * t
seasonal = 30 * np.sin(2 * np.pi * t / 12) + 15 * np.cos(2 * np.pi * t / 12)
noise = np.random.randn(144) * 10
passengers = trend + seasonal + noise
passengers = np.maximum(passengers, 50)
ts = pd.Series(passengers, index=months, name='Passengers')
# 1階差分系列の ACF/PACF
ts_diff = ts.diff().dropna()
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
plot_acf(ts_diff, lags=40, ax=axes[0], title='ACF of First Difference')
plot_pacf(ts_diff, lags=40, ax=axes[1], title='PACF of First Difference',
method='ywm')
plt.tight_layout()
plt.show()
ACF プロットでは、ラグ 12, 24 付近にスパイクが見えれば季節性が残っていることを示します。PACF プロットでは、AR 次数の候補がスパイクの切断位置から読み取れます。
AIC / BIC によるモデル選択
import numpy as np
import pandas as pd
from statsmodels.tsa.arima.model import ARIMA
import warnings
warnings.filterwarnings('ignore')
# データの再生成
np.random.seed(42)
months = pd.date_range('1949-01', periods=144, freq='MS')
t = np.arange(144)
trend = 100 + 2.5 * t
seasonal = 30 * np.sin(2 * np.pi * t / 12) + 15 * np.cos(2 * np.pi * t / 12)
noise = np.random.randn(144) * 10
passengers = trend + seasonal + noise
passengers = np.maximum(passengers, 50)
ts = pd.Series(passengers, index=months, name='Passengers')
# グリッドサーチで最適な ARIMA(p,d,q) を探索
results = []
for p in range(4):
for d in range(2):
for q in range(4):
try:
model = ARIMA(ts, order=(p, d, q))
fitted = model.fit()
results.append({
'order': (p, d, q),
'AIC': fitted.aic,
'BIC': fitted.bic,
'Log-Likelihood': fitted.llf
})
except:
continue
results_df = pd.DataFrame(results)
results_df = results_df.sort_values('AIC')
print("=== ARIMA モデル選択（AIC順） ===")
print(results_df.head(10).to_string(index=False))
print()
# 最良モデル
best_order = results_df.iloc[0]['order']
print(f"AIC 最小モデル: ARIMA{best_order}")
SARIMA モデルの推定と予測
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings
warnings.filterwarnings('ignore')
# データの再生成
np.random.seed(42)
months = pd.date_range('1949-01', periods=144, freq='MS')
t = np.arange(144)
trend = 100 + 2.5 * t
seasonal = 30 * np.sin(2 * np.pi * t / 12) + 15 * np.cos(2 * np.pi * t / 12)
noise = np.random.randn(144) * 10
passengers = trend + seasonal + noise
passengers = np.maximum(passengers, 50)
ts = pd.Series(passengers, index=months, name='Passengers')
# 学習データとテストデータの分割
train = ts[:'1959-12']
test = ts['1960-01':]
# SARIMA モデルの推定
model = SARIMAX(train,
order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False)
fitted = model.fit(disp=False)
print(fitted.summary())
# 予測
forecast = fitted.get_forecast(steps=len(test))
forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int(alpha=0.05)
# 予測結果の可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 予測と実測値の比較
axes[0, 0].plot(train, 'b-', label='Train', linewidth=1)
axes[0, 0].plot(test, 'g-', label='Test (Actual)', linewidth=1.5)
axes[0, 0].plot(forecast_mean, 'r--', label='Forecast', linewidth=1.5)
axes[0, 0].fill_between(forecast_ci.index,
forecast_ci.iloc[:, 0],
forecast_ci.iloc[:, 1],
color='red', alpha=0.1, label='95% CI')
axes[0, 0].set_title('SARIMA(1,1,1)(1,1,1)[12] Forecast', fontsize=13)
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
# (2) 残差の時系列プロット
residuals = fitted.resid
axes[0, 1].plot(residuals, 'k-', linewidth=0.8)
axes[0, 1].axhline(y=0, color='r', linestyle='--')
axes[0, 1].set_title('Residuals', fontsize=13)
axes[0, 1].grid(True, alpha=0.3)
# (3) 残差の ACF
plot_acf(residuals, lags=30, ax=axes[1, 0],
title='ACF of Residuals')
# (4) 残差のヒストグラム
axes[1, 1].hist(residuals, bins=25, density=True, alpha=0.7, color='steelblue')
x_range = np.linspace(residuals.min(), residuals.max(), 100)
from scipy.stats import norm
mu, sigma = norm.fit(residuals)
axes[1, 1].plot(x_range, norm.pdf(x_range, mu, sigma), 'r-', linewidth=2,
label=f'N({mu:.1f}, {sigma:.1f}²)')
axes[1, 1].set_title('Residual Distribution', fontsize=13)
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Ljung-Box 検定
lb_test = acorr_ljungbox(residuals, lags=[10, 20, 30], return_df=True)
print("\n=== Ljung-Box 検定 ===")
print(lb_test)
AR(p) の定常性条件の可視化
import numpy as np
import matplotlib.pyplot as plt
# AR(2) の定常性領域を可視化
phi1 = np.linspace(-2.5, 2.5, 500)
phi2 = np.linspace(-1.5, 1.5, 500)
PHI1, PHI2 = np.meshgrid(phi1, phi2)
# 定常性条件: phi1 + phi2 < 1, phi2 - phi1 < 1, |phi2| < 1
cond1 = PHI1 + PHI2 < 1
cond2 = PHI2 - PHI1 < 1
cond3 = np.abs(PHI2) < 1
stationary = cond1 & cond2 & cond3
plt.figure(figsize=(8, 6))
plt.contourf(PHI1, PHI2, stationary.astype(float),
levels=[0.5, 1.5], colors=['lightblue'], alpha=0.5)
plt.contour(PHI1, PHI2, stationary.astype(float),
levels=[0.5], colors=['blue'], linewidths=2)
# 頂点を表示
vertices = [(1, 0), (-1, 0), (0, -1)]
for v in vertices:
plt.plot(*v, 'ro', markersize=8)
plt.annotate(f'({v[0]}, {v[1]})', xy=v, fontsize=10,
xytext=(v[0]+0.15, v[1]+0.08))
# 境界線の方程式
plt.plot(phi1, 1 - phi1, 'b--', alpha=0.5, label=r'$\phi_1 + \phi_2 = 1$')
plt.plot(phi1, 1 + phi1, 'g--', alpha=0.5, label=r'$\phi_2 - \phi_1 = 1$')
plt.axhline(y=1, color='r', linestyle='--', alpha=0.5, label=r'$\phi_2 = \pm 1$')
plt.axhline(y=-1, color='r', linestyle='--', alpha=0.5)
plt.xlabel(r'$\phi_1$', fontsize=14)
plt.ylabel(r'$\phi_2$', fontsize=14)
plt.title('Stationarity Region of AR(2)', fontsize=14)
plt.legend(loc='upper left', fontsize=9)
plt.grid(True, alpha=0.3)
plt.xlim(-2.5, 2.5)
plt.ylim(-1.5, 1.5)
plt.tight_layout()
plt.show()
この図は、AR(2) モデル $x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \varepsilon_t$ のパラメータ空間において、定常性が保証される三角形の領域を示しています。3つの定常性条件（$\phi_1 + \phi_2 < 1$, $\phi_2 - \phi_1 < 1$, $|\phi_2| < 1$）の交差部分が定常領域です。
まとめ
本記事では、ARIMA モデルの理論と実装について解説しました。
- AR(p) モデル は過去 $p$ 個の値に依存し、特性方程式の根が単位円外にあるとき定常である
- MA(q) モデル は過去 $q$ 個のホワイトノイズに依存し、ACF がラグ $q$ で切断される
- ARIMA(p,d,q) は非定常時系列に $d$ 階差分を適用して ARMA でモデル化する
- Box-Jenkins 法 により同定→推定→診断のサイクルで体系的にモデルを構築できる
- ACF/PACF は次数同定の基本ツールであり、AIC/BIC が複数モデルの比較に有効
- SARIMA は季節性を持つ時系列に対応した拡張であり、多くの実データに適用可能
次のステップとして、以下の記事も参考にしてください。