ARMA・ARIMAモデル｜時系列分析

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：AR・MAモデル　|　関連：ランダムウォークと単位根・予測の評価指標と時系列CV

要点（BLUF）

ARMA( $p,q$ )＝AR（過去の値）と MA（過去のショック）を両方載せたモデル。少ないパラメータで複雑な自己相関を表せます（倹約性）。
ARIMA( $p,d,q$ )＝非定常系列を $d$ 回差分して定常化し（I＝和分、ランダムウォークと単位根）、その差分系列に ARMA を当てるモデル。 $d$ は差分回数。
予測は点予測だけでなく予測区間を出します。和分系列では予測の不確実性が累積するので、区間は先のホライズンほど広がります。最後はウォークフォワードで素朴予測と比べて評価します。

1. ARMA( $p,q$ )：AR と MA を統合する

ARMA( $p,q$ ) は、AR 項と MA 項を同時に持ちます。

y_t = c + \sum_{i=1}^{p}\phi_i y_{t-i} + \varepsilon_t + \sum_{j=1}^{q}\theta_j\varepsilon_{t-j} \quad\Longleftrightarrow\quad \phi(L)\,y_t = c + \theta(L)\,\varepsilon_t

ここで $\phi(L)=1-\phi_1L-\cdots-\phi_pL^p$ 、 $\theta(L)=1+\theta_1L+\cdots+\theta_qL^q$ （AR・MAモデル）。定常性は AR 部分（ $\phi(z)=0$ の根が単位円の外）、反転可能性は MA 部分（ $\theta(z)=0$ の根が単位円の外）が担います。

ARMA の利点は倹約性（parsimony）。純粋な AR や MA だと高次が必要な自己相関でも、ARMA なら低次の $(p,q)$ で表せることがあります。代償として、ARMA は ACF も PACF も尾を引いて切れないため、 $(p,q)$ を図だけで読みにくい——そこで AIC などの情報量規準で選びます（モデル選択と残差診断で扱います）。

コード①：ARMA(1,1) の φ と θ を同時に復元

真の $\phi=0.6,\ \theta=0.4$ を仕込んだ ARMA(1,1) を生成し、ARIMA(order=(1,0,1)) で AR 係数と MA 係数の両方を復元します。

import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima_process import arma_generate_sample

# 真の ARMA(1,1): y_t = 0.6 y_{t-1} + ε_t + 0.4 ε_{t-1}
phi_true, theta_true = 0.6, 0.4
np.random.seed(7)
y = arma_generate_sample(np.array([1, -phi_true]), np.array([1, theta_true]),
                         nsample=3000, scale=1.0)

res = ARIMA(y, order=(1, 0, 1), trend="c").fit()
phi_hat = res.params[res.param_names.index("ar.L1")]
theta_hat = res.params[res.param_names.index("ma.L1")]
print(f"{'係数':>6}{'真値':>10}{'推定値':>10}")
print(f"{'φ(AR)':>6}{phi_true:>10.3f}{phi_hat:>10.3f}")
print(f"{'θ(MA)':>6}{theta_true:>10.3f}{theta_hat:>10.3f}")
print(f"AIC={res.aic:.1f}")

出力：

    係数        真値       推定値
 φ(AR)     0.600     0.593
 θ(MA)     0.400     0.419
AIC=8483.8

出力の意味：AR 係数 $\hat\phi=0.593$ 、MA 係数 $\hat\theta=0.419$ と、仕込んだ $0.6,\ 0.4$ を両方とも復元しました。ARMA は AR と MA が混ざっても、最尤推定で分離して当てられます。AIC は次数選択の指標——候補モデル間で比べ、小さいものを選びます（絶対値そのものに意味はなく、相対比較に使う）。

2. ARIMA( $p,d,q$ )：差分（I）で非定常を取り込む

現実の系列はトレンドや単位根で非定常なことが多い（ランダムウォークと単位根）。ARMA は定常系列のモデルなので、そのままでは当てられません。そこで $d$ 階差分で定常化してから ARMA を当てるのが ARIMA( $p,d,q$ ) です。

\phi(L)\,(1-L)^d\,y_t = c + \theta(L)\,\varepsilon_t

$(1-L)y_t=y_t-y_{t-1}$ が1階差分。 $d$ は差分回数で、 $(1-L)^d y_t$ が定常な ARMA( $p,q$ ) になる、という構造です。たとえばランダムウォークは ARIMA(0,1,0)（1階差分するとホワイトノイズ）。 $d$ は ADF 検定や、差分前後の ACF の減衰（定常性と自己相関）で決めます。予測の流れは「差分系列で ARMA 予測 → 逆差分（累積）で水準へ戻す」。逆差分のとき不確実性も累積するので、予測区間が先のホライズンほど広がります。

flowchart LR
    A["原系列 y_t（非定常・トレンド）"] --> B["d 階差分 (1-L)^d で定常化"]
    B --> C["定常系列に ARMA(p,q) を推定"]
    C --> D["差分系列を予測 + 予測分散"]
    D --> E["逆差分（累積）で水準へ戻す"]
    E --> F["点予測 + 予測区間（先ほど広がる）"]

コード②：トレンド付き非定常系列を ARIMA(1,1,1) で予測（予測区間つき）

ドリフト付きランダムウォーク（＋AR(1) 的な揺れ）という非定常な上昇トレンド系列を作り、ARIMA(order=(1,1,1))（ $d=1$ ）で当てて 30 期先を予測します。点予測だけでなく 95% 予測区間（get_forecast().conf_int()）を出し、区間幅がホライズンとともに広がることを数値と図で確かめます。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib  # 日本語ラベル用
from statsmodels.tsa.arima.model import ARIMA

# トレンドを持つ非定常系列：ドリフト付きランダムウォーク + AR(1)的な揺れ
# Δy_t = 0.3(ドリフト) + u_t,  u_t = 0.5 u_{t-1} + ε_t  → y は I(1)
np.random.seed(11)
n = 200
eps = np.random.normal(0, 1.0, n)
u = np.zeros(n)
for t in range(1, n):
    u[t] = 0.5 * u[t-1] + eps[t]
drift = 0.3
dy = drift + u                 # 1階差分（定常）
y = 10 + np.cumsum(dy)         # 原系列（非定常・上昇トレンド）

# ARIMA(1,1,1)：d=1 で差分、trend='t' でドリフトを推定
model = ARIMA(y, order=(1, 1, 1), trend="t")
res = model.fit()

# 30 期先を予測（点予測 + 95% 予測区間）
fc = res.get_forecast(steps=30)
mean = fc.predicted_mean
ci = fc.conf_int(alpha=0.05)   # shape (30, 2): [下限, 上限]
print(f"1期先 予測区間幅 = {ci[0,1]-ci[0,0]:.2f}")
print(f"30期先 予測区間幅 = {ci[29,1]-ci[29,0]:.2f}")
print(f"区間は先のホライズンほど {'広がる' if (ci[29,1]-ci[29,0]) > (ci[0,1]-ci[0,0]) else '広がらない'}")

future_idx = np.arange(n, n+30)
plt.figure(figsize=(9, 4.5))
plt.plot(np.arange(n), y, color="gray", lw=1, label="実測")
plt.plot(future_idx, mean, color="C1", lw=2, label="点予測")
plt.fill_between(future_idx, ci[:, 0], ci[:, 1], color="C1", alpha=0.25, label="95%予測区間")
plt.axvline(n-1, ls=":", color="k")
plt.xlabel("時点 t"); plt.ylabel("y"); plt.legend()
plt.title("ARIMA(1,1,1) による予測と予測区間（先ほど広がる）")
plt.tight_layout(); plt.show()

出力：

1期先 予測区間幅 = 3.66
30期先 予測区間幅 = 35.44
区間は先のホライズンほど 広がる

出力の意味：1期先の区間幅は $3.66$ 、30期先では $35.44$ ——約10倍に広がりました。点予測（橙線）はドリフトを引き継いで上昇を続けますが、その周りの不確実性は急速に膨らむ。図では橙の帯がラッパ状に開きます。これは「和分系列では各ステップのショック分散が累積するから」（次節）。先の予測ほど自信を持てない——点予測だけ見て将来を断言するのは危険で、必ず区間で語るべきだ、というのが時系列予測の鉄則です。

3. なぜ予測区間は先ほど広がるのか

最も単純な ARIMA(0,1,0)＝ランダムウォーク $y_t=y_{t-1}+\varepsilon_t$ で考えると直観的です。現在 $t$ から $h$ 期先は、

y_{t+h}=y_t+\varepsilon_{t+1}+\varepsilon_{t+2}+\cdots+\varepsilon_{t+h}

なので、予測の分散は独立ショックの和で $\mathrm{Var}[y_{t+h}\mid y_t]=h\,\sigma^2$ 。標準偏差は $\sigma\sqrt{h}$ 、95% 区間幅は $2\times1.96\,\sigma\sqrt{h}$ と** $\sqrt{h}$ に比例して広がります**。一般の ARIMA でも、差分系列の予測誤差を逆差分（累積）する過程で誤差が積み上がり、ホライズンとともに区間が広がります。一方、定常な ARMA（ $d=0$ ）の予測は平均回帰し、区間幅は有限値（無条件分散）に収束します——非定常（ $d\ge1$ ）だけが青天井に広がる点が対照的です。

コード③：ウォークフォワードで1期先予測を評価（RMSE・区間カバレッジ）

真の AR(2)（定常・実構造あり）を生成し、後半をウォークフォワード（過去だけで1期先を予測 → 1歩進める）でバックテストします（予測の評価指標と時系列CV）。RMSE を素朴予測（前回値）と比較し、95% 予測区間が実際に約95%をカバーするか（較正）も確かめます。

import numpy as np
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.arima_process import arma_generate_sample

# 真の AR(2)（定常・実構造あり）を生成
np.random.seed(21)
y = arma_generate_sample(np.array([1, -0.5, 0.3]), np.array([1.0]),
                         nsample=300, scale=1.0)

start = 250                     # ここから1期ずつ前進検証（ウォークフォワード）
err_arima, err_naive, covered = [], [], 0
for i in range(start, len(y)):
    hist = y[:i]                # 過去だけ（未来は見ない）
    res = ARIMA(hist, order=(2, 0, 0), trend="c").fit()
    fc = res.get_forecast(steps=1)
    pred = fc.predicted_mean[0]
    lo, hi = fc.conf_int(alpha=0.05)[0]
    actual = y[i]
    err_arima.append(pred - actual)
    err_naive.append(hist[-1] - actual)   # 素朴予測：前回値
    if lo <= actual <= hi:
        covered += 1

err_arima = np.array(err_arima); err_naive = np.array(err_naive)
rmse_arima = np.sqrt(np.mean(err_arima**2))
rmse_naive = np.sqrt(np.mean(err_naive**2))
nstep = len(err_arima)
print(f"1期先 RMSE: ARIMA(2,0,0)={rmse_arima:.3f}  素朴(前回値)={rmse_naive:.3f}")
print(f"改善率 = {(1 - rmse_arima/rmse_naive)*100:.1f}%（素朴比）")
print(f"95%予測区間カバレッジ = {covered}/{nstep} = {covered/nstep*100:.0f}%")

出力：

1期先 RMSE: ARIMA(2,0,0)=0.949  素朴(前回値)=1.244
改善率 = 23.7%（素朴比）
95%予測区間カバレッジ = 48/50 = 96%

出力の意味：ARIMA(2,0,0) の RMSE は $0.949$ で、素朴予測（前回値）の $1.244$ より23.7% 良い。系列に本物の AR(2) 構造があるので、それを学んだモデルがベースラインを上回りました——これで初めて「モデルに価値がある」と言えます（ベースライン未満なら無意味、予測の評価指標と時系列CV）。さらに 95% 予測区間が 50 回中 48 回（96%）実測を捕らえ、名目の 95% とよく一致——区間が較正されている証拠です。ウォークフォワードでは毎回「過去だけ」で再推定するので、未来情報の漏洩がありません。

4. 数式の直観

ARMA は AR の慣性と MA のショックの尾の合成： $\phi(L)y_t=\theta(L)\varepsilon_t$ を $y_t=\frac{\theta(L)}{\phi(L)}\varepsilon_t$ と見ると、低次の有理多項式 $\theta/\phi$ で豊かな MA( $\infty$ ) を近似している——これが倹約性の正体です。
I（差分）はトレンドの除去装置： $(1-L)^d$ は「水準の情報」を落として「変化の情報」だけにする操作。差分系列が定常な ARMA になれば、ARIMA の枠で予測でき、最後に累積して水準へ戻します。
予測区間＝累積したショックの不確実性：点予測は条件付き期待値、区間は条件付き分散から。和分系列ではこの分散が $h$ とともに増えるので、区間が広がります。区間はモデルが正しく、誤差が（おおむね）正規という前提のもとでの目安です。

⚠️ よくある誤解

「点予測だけ出せば十分」ではない：時系列予測は不確実性が本質。必ず予測区間を添え、先ほど広がることを前提に意思決定します。
「 $d$ は大きいほど安全」ではない：過剰差分は人工的な負の自己相関（MA 的構造）を生み、予測を悪化させます。 $d$ は最小限（ふつう 0〜2）に。ADF と ACF で見極めます（ランダムウォークと単位根）。
「ARIMA は遠い未来も当たる」ではない：定常 ARMA は数ステップ先で平均（や差分系列の平均＝ドリフト）へ回帰します。長期は外挿で、区間も急速に広がる——遠未来の点予測を過信しないこと。
「区間に入る割合は気にしなくてよい」ではない：95% 区間が実測の 95% を捕らえているか（カバレッジ）で較正を点検します。大きく外れるなら誤差分布の仮定（正規・等分散）を疑います。
「ARMA の次数は ACF/PACF だけで決まる」ではない：ARMA は ACF も PACF も切れないので、図に加えて AIC など情報量規準で選びます（モデル選択と残差診断）。