偏自己相関(PACF)とARモデルの次数同定:Yule-Walker・Levinson-DurbinのPython実装

偏自己相関(PACF)とYule-Walker方程式・Levinson-Durbin再帰の理論とPython実装を体系解説。statsmodels の pacf・plot_acf/plot_pacf・levinson_durbin・AutoReg・ar_select_order・arma_order_select_ic・ARIMA を使い、AR(2)過程でPACFがラグ2で切断すること(PACF(3)以降が95%信頼区間±0.088内)と、AIC/BICによる次数選択が(p,q)=(2,0)でPACFの判定と一致することを数値実験で検証します。

はじめに

自己相関(ACF)を計算すると、時系列が「過去とどれだけ似ているか」がラグごとに分かります。しかし ACF だけでは次の問いに答えられません——「ラグ 2 の相関は、本当にラグ 2 の直接的な影響なのか? それともラグ 1 の影響が 2 ステップ伝播しただけなのか?」

この問いに答えるのが偏自己相関関数(PACF: Partial Autocorrelation Function)です。PACF は中間ラグの影響を取り除いた「純粋な」相関を測るため、ARモデルの次数 \(p\) がそのままPACF が非ゼロになる最大ラグとして現れます。この性質が、ARIMAモデルの次数決定で ACF/PACF プロットを眺める理由の理論的根拠です。

本記事は自己相関の基礎編に対する発展編として、PACF の定義(線形予測残差の相関)、Yule-Walker 方程式、Levinson-Durbin 再帰を数式で導出し、AR/MA の次数同定の判定表と AIC/BIC による情報量規準を整理します。そのうえで statsmodels(pacf, plot_acf/plot_pacf, levinson_durbin, AutoReg, ar_select_order, arma_order_select_ic, ARIMA)を使い、既知の AR(2) 過程で「PACF の切断」と「AIC/BIC の最小化」が同じ次数を指すことを数値実験で確認します。

偏自己相関(PACF)の定義

線形予測残差の相関としての定義

定常時系列 \(\{x_t\}\) のラグ \(k\) の偏自己相関は、中間の値 \(x_{t-1}, \ldots, x_{t-k+1}\) の線形効果を両側から除去した後の相関として定義されます。

\[\alpha(k) = \operatorname{Corr}\left(x_t - \hat{x}_t^{(k-1)},\; x_{t-k} - \hat{x}_{t-k}^{(k-1)}\right) \tag{1}\]

ここで \(\hat{x}_t^{(k-1)}\) は \(x_{t-1}, \ldots, x_{t-k+1}\) による \(x_t\) の最良線形予測、\(\hat{x}_{t-k}^{(k-1)}\) は同じ中間値による \(x_{t-k}\) の最良線形予測(後ろ向き予測)です。つまり PACF は「中間ラグで説明できる部分を差し引いた残差同士の相関」であり、回帰分析の偏相関係数の時系列版です。

AR回帰の最終係数としての定義

式 \((1)\) と等価なもう一つの定義があります。次数 \(k\) の自己回帰

\[x_t = \phi_{k1} x_{t-1} + \phi_{k2} x_{t-2} + \cdots + \phi_{kk} x_{t-k} + \varepsilon_t \tag{2}\]

を当てはめたときの最終係数 \(\phi_{kk}\) が、ちょうどラグ \(k\) の偏自己相関になります。

\[\alpha(k) = \phi_{kk} \tag{3}\]

「ラグ \(1\) から \(k-1\) までを説明変数に入れてもなお、\(x_{t-k}\) が追加でどれだけ説明力を持つか」を測っているので、直観的にも式 \((1)\) と同じものです。なお \(\alpha(1)\) は除去すべき中間ラグがないため、通常の自己相関 \(\rho_1\) と一致します。

なぜAR過程でPACFが切断するのか

AR(1) 過程 \(x_t = \phi x_{t-1} + \varepsilon_t\) を考えます。ACF は \(\rho_k = \phi^k\) と指数的に減衰し、ラグ 2 でも \(\rho_2 = \phi^2 \neq 0\) です。しかしこの相関は「\(x_{t-2} \to x_{t-1} \to x_t\) 」と伝播した間接効果に過ぎません。\(x_{t-1}\) の効果を除去すると \(x_{t-2}\) の直接の寄与は残らないため、\(\alpha(2) = 0\) になります。

一般に AR(p) 過程では、\(x_t\) は直近 \(p\) 個の値だけで決まる(式 \((2)\) で \(k > p\) なら \(\phi_{kk} = 0\) )ため、

\[\alpha(k) = 0 \quad (k > p) \tag{4}\]

が成り立ちます。PACF がラグ \(p\) で切断(cut off)する——これが AR モデルの次数同定の核心です。

Yule-Walker 方程式

PACF を実際に計算するには式 \((2)\) の係数 \(\phi_{kj}\) が必要です。AR(p) 過程

\[x_t = \phi_1 x_{t-1} + \cdots + \phi_p x_{t-p} + \varepsilon_t, \quad \varepsilon_t \sim \mathrm{WN}(0, \sigma^2) \tag{5}\]

の両辺に \(x_{t-k}\) (\(k \geq 1\) )を掛けて期待値を取ると、\(E[\varepsilon_t x_{t-k}] = 0\) より自己相関 \(\rho_k\) に関する線形方程式が得られます。

\[\rho_k = \phi_1 \rho_{k-1} + \phi_2 \rho_{k-2} + \cdots + \phi_p \rho_{k-p}, \quad k = 1, \ldots, p \tag{6}\]

これが Yule-Walker 方程式です。\(k = 1, \ldots, p\) をまとめて行列形式で書くと

\[ \begin{pmatrix} 1 & \rho_1 & \cdots & \rho_{p-1} \\ \rho_1 & 1 & \cdots & \rho_{p-2} \\ \vdots & \vdots & \ddots & \vdots \\ \rho_{p-1} & \rho_{p-2} & \cdots & 1 \end{pmatrix} \begin{pmatrix} \phi_1 \\ \phi_2 \\ \vdots \\ \phi_p \end{pmatrix} = \begin{pmatrix} \rho_1 \\ \rho_2 \\ \vdots \\ \rho_p \end{pmatrix} \tag{7} \]

となります。左辺の行列は対角線に沿って同じ値が並ぶテプリッツ行列です。標本 ACF を代入して式 \((7)\) を解けば AR 係数の推定値(Yule-Walker 推定)が得られ、これを \(k = 1, 2, 3, \ldots\) と次数を増やしながら解いたときの各最終係数 \(\phi_{kk}\) の列が標本 PACF になります。

Levinson-Durbin 再帰

次数 \(k\) ごとに式 \((7)\) を独立に解くと \(O(k^3)\) の逆行列計算を繰り返すことになります。テプリッツ構造を利用して、次数 \(k-1\) の解から次数 \(k\) の解を \(O(k)\) で更新するのが Levinson-Durbin 再帰です。

初期値を \(\phi_{11} = \rho_1\) , \(\sigma_1^2 = \gamma_0(1 - \rho_1^2)\) として、\(k = 2, 3, \ldots\) について次を繰り返します。

\[\phi_{kk} = \frac{\rho_k - \sum_{j=1}^{k-1} \phi_{k-1,j}\, \rho_{k-j}}{1 - \sum_{j=1}^{k-1} \phi_{k-1,j}\, \rho_j} \tag{8}\] \[\phi_{kj} = \phi_{k-1,j} - \phi_{kk}\, \phi_{k-1,k-j}, \quad j = 1, \ldots, k-1 \tag{9}\] \[\sigma_k^2 = \sigma_{k-1}^2 \left(1 - \phi_{kk}^2\right) \tag{10}\]

ここで \(\sigma_k^2\) は次数 \(k\) の線形予測の残差分散(イノベーション分散)です。3 つの式にはそれぞれ明確な意味があります。

  • 式 \((8)\) :分子は「現在のモデルで説明しきれないラグ \(k\) の相関」、分母は正規化項。この \(\phi_{kk}\) が**反射係数(reflection coefficient)**と呼ばれ、そのまま PACF になる
  • 式 \((9)\) :次数を 1 つ上げたときの既存係数の補正
  • 式 \((10)\) :\(|\phi_{kk}| < 1\) である限り残差分散は単調に減少し、真の次数 \(p\) を超えると \(\phi_{kk} \approx 0\) となって減少が止まる

最大次数 \(p\) まで回して合計 \(O(p^2)\) で、全次数の AR 係数と PACF が一度に得られます。この再帰は音声符号化の線形予測分析(LPC)やラティスフィルタの設計と同じ数学であり、時系列解析と信号処理をつなぐ蝶番になっています。

次数同定の判定表

AR と MA では ACF と PACF の役割がちょうど入れ替わります。MA(q) 過程 \(x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}\) は \(q\) ステップより離れた項と共有するノイズを持たないため、ACF がラグ \(q\) で切断します。一方 MA 過程を AR 表現に直すと無限次数になるため、PACF は切断せず減衰します。

過程ACFPACF
AR(p)減衰(指数・振動)ラグ \(p\) で切断
MA(q)ラグ \(q\) で切断減衰(指数・振動)
ARMA(p, q)減衰(ラグ \(q\) 以降)減衰(ラグ \(p\) 以降)
白色雑音全ラグで \(\approx 0\)全ラグで \(\approx 0\)

「切断」の判定には信頼区間を使います。真の値が 0 のラグでは、標本 PACF は漸近的に \(\mathcal{N}(0, 1/N)\) に従うため、

\[|\hat{\alpha}(k)| > \frac{1.96}{\sqrt{N}} \tag{11}\]

を有意とみなします(95%信頼区間)。ただし約 5% のラグは偶然この区間を超えるので、高ラグの単発の有意は無視するのが実務的です。

情報量規準(AIC / BIC)による次数選択

ACF/PACF プロットによる判定は視覚的で分かりやすい一方、切断が曖昧なケース(ARMA 混合や小標本)では主観が入ります。これを補うのが情報量規準です。パラメータ数 \(k\) 、標本サイズ \(N\) 、最大対数尤度 \(\ln \hat{L}\) に対して

\[\mathrm{AIC} = -2 \ln \hat{L} + 2k \tag{12}\] \[\mathrm{BIC} = -2 \ln \hat{L} + k \ln N \tag{13}\]

を候補次数ごとに計算し、最小の次数を選びます。第 1 項(当てはまりの良さ)と第 2 項(複雑さへの罰則)のトレードオフで、\(\ln N > 2\) (\(N \geq 8\) )では BIC の罰則の方が重いため、BIC は AIC より小さい次数を選ぶ傾向があります。理論的には、真のモデルが候補に含まれるとき BIC は一致性(\(N \to \infty\) で真の次数を選ぶ確率が 1)を持ち、AIC は予測誤差最小化の意味で効率的です。実務では両方計算し、一致すればその次数、割れたら「予測目的なら AIC、構造解釈なら BIC」を目安にします。

注意点として、モデル比較は同じデータ点に対する尤度で行う必要があります。AR(p) は先頭 \(p\) 点を初期値に使うため、次数ごとに素朴に当てはめると標本が変わってしまいます。後述の ar_select_order はこの点を揃えて比較してくれます。

Python実装:AR(2)過程での数値実験

AR(2)の生成とACF/PACFの計算

真のモデルが分かっている AR(2) 過程 \(x_t = 0.7 x_{t-1} - 0.3 x_{t-2} + \varepsilon_t\) (\(\varepsilon_t \sim \mathcal{N}(0, 1)\) , \(N = 500\) )を生成し、PACF がラグ 2 で切断することを確認します。

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

np.random.seed(42)
N = 500
phi = np.array([0.7, -0.3])      # 真のAR係数
ar = np.r_[1, -phi]              # statsmodels の符号規約: [1, -phi1, -phi2]
ma = np.r_[1]
x = ArmaProcess(ar, ma).generate_sample(nsample=N)

nlags = 20
acf_vals = acf(x, nlags=nlags, fft=True)
pacf_vals = pacf(x, nlags=nlags, method="ywadjusted")
ci = 1.96 / np.sqrt(N)           # 95%信頼区間

print(f"95%信頼区間: ±{ci:.3f}")
print(f"{'lag':>4} {'ACF':>8} {'PACF':>8}")
for k in range(1, 6):
    line = f"{k:>4} {acf_vals[k]:>8.3f} {pacf_vals[k]:>8.3f}"
    if abs(pacf_vals[k]) > ci:
        line += " *"
    print(line)

# 定番の可視化: plot_acf / plot_pacf
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(x, lags=nlags, ax=axes[0], title="ACF: gradual decay")
plot_pacf(x, lags=nlags, ax=axes[1], method="ywm",
          title="PACF: cuts off after lag 2")
plt.tight_layout()
plt.show()

実行結果:

95%信頼区間: ±0.088
 lag      ACF     PACF
   1    0.532    0.533 *
   2    0.060   -0.312 *
   3   -0.136   -0.020
   4   -0.135   -0.020
   5   -0.013    0.079

判定表どおりの結果です。ACF は \(0.532 \to 0.060 \to -0.136\) と符号を変えながら減衰し(\(\phi_2 < 0\) による振動減衰)、明確な切断を示しません。一方 PACF は \(\hat{\alpha}(1) = 0.533\) , \(\hat{\alpha}(2) = -0.312\) が有意で、\(\hat{\alpha}(3)\) 以降(\(-0.020, -0.020, 0.079\) )はすべて信頼区間 \(\pm 0.088\) の内側——ラグ 2 での切断が確認でき、AR(2) と同定できます。\(\hat{\alpha}(2) = -0.312\) が真の \(\phi_2 = -0.3\) の推定になっている点にも注目してください(式 \((3)\) で \(k = p\) のとき \(\phi_{pp} = \phi_p\) )。

Levinson-Durbin再帰でPACFとイノベーション分散を確認

levinson_durbin は再帰の副産物(全次数の PACF・イノベーション分散)をまとめて返します。

from statsmodels.tsa.stattools import levinson_durbin

sigma_v, arcoefs, pacf_ld, sigma, _ = levinson_durbin(x, nlags=5, isacov=False)
print("PACF (Levinson-Durbin):", np.round(pacf_ld[1:6], 3))
print("innovation variance:", np.round(sigma[1:6], 3))

実行結果:

PACF (Levinson-Durbin): [ 0.532 -0.311 -0.02  -0.02   0.078]
innovation variance: [1.067 0.964 0.963 0.963 0.957]

pacf(method="ywadjusted") とは ACF の推定法(不偏補正の有無)が異なるため 3 桁目がわずかにずれますが、実質同じ値です。式 \((10)\) のとおりイノベーション分散は次数 1 → 2 で \(1.067 \to 0.964\) と大きく下がり、次数 3 以降はほぼ横ばい(\(0.963, 0.963, 0.957\) )。「次数を 2 から増やしても予測は良くならない」ことが残差分散の側からも読み取れます。

AIC/BICによる次数選択(AutoReg / ar_select_order)

次数 \(p = 0, \ldots, 6\) の AR モデルを同一標本(hold_back=10 で先頭 10 点を除外して条件を揃える)に当てはめ、AIC/BIC を比較します。

from statsmodels.tsa.ar_model import AutoReg, ar_select_order

print(f"{'p':>2} {'AIC':>9} {'BIC':>9}")
for p in range(0, 7):
    res = AutoReg(x, lags=p, hold_back=10, old_names=False).fit()
    print(f"{p:>2} {res.aic:>9.2f} {res.bic:>9.2f}")

sel_aic = ar_select_order(x, maxlag=10, ic="aic", old_names=False)
sel_bic = ar_select_order(x, maxlag=10, ic="bic", old_names=False)
print("AIC選択ラグ:", sel_aic.ar_lags)
print("BIC選択ラグ:", sel_bic.ar_lags)

res2 = AutoReg(x, lags=2, old_names=False).fit()
print("AutoReg(2) 係数:", np.round(res2.params, 3))

実行結果:

 p       AIC       BIC
 0   1592.91   1601.30
 1   1430.99   1443.57
 2   1381.90   1398.68
 3   1383.54   1404.51
 4   1385.38   1410.55
 5   1384.62   1413.98
 6   1385.15   1418.71
AIC選択ラグ: [1, 2]
BIC選択ラグ: [1, 2]
AutoReg(2) 係数: [ 0.006  0.702 -0.313]

AIC・BIC とも \(p = 2\) (AIC \(= 1381.90\) , BIC \(= 1398.68\) )で最小となり、ar_select_order もラグ \(\{1, 2\}\) を選択。PACF の切断(ラグ 2)と情報量規準の選択(\(p = 2\) )が一致しました。推定係数 \((\hat{\phi}_1, \hat{\phi}_2) = (0.702, -0.313)\) も真値 \((0.7, -0.3)\) を高精度に回復しています。

ARMA全体での次数探索(arma_order_select_ic / ARIMA)

MA 項も候補に含めた場合でも AR(2) が選ばれるかを arma_order_select_ic で確認し、最終モデルを ARIMA で最尤推定します。

from statsmodels.tsa.stattools import arma_order_select_ic
from statsmodels.tsa.arima.model import ARIMA

order_sel = arma_order_select_ic(x, max_ar=4, max_ma=2, ic=["aic", "bic"])
print("AIC最小の次数 (p, q):", tuple(map(int, order_sel.aic_min_order)))
print("BIC最小の次数 (p, q):", tuple(map(int, order_sel.bic_min_order)))

arima_res = ARIMA(x, order=(2, 0, 0)).fit()
print("ARIMA(2,0,0) AR係数:", np.round(arima_res.arparams, 3))
print("ARIMA(2,0,0) 残差分散:", round(arima_res.params[-1], 3))

実行結果:

AIC最小の次数 (p, q): (2, 0)
BIC最小の次数 (p, q): (2, 0)
ARIMA(2,0,0) AR係数: [ 0.7   -0.312]
ARIMA(2,0,0) 残差分散: 0.96

\((p, q)\) の 2 次元グリッド(\(5 \times 3 = 15\) モデル)を探索しても、AIC・BIC とも \((2, 0)\) を選びます。最尤推定の AR 係数 \((0.700, -0.312)\) 、残差分散 \(0.96\) (真値 \(1.0\) )も整合的です。PACF・AIC・BIC という 3 つの独立な手がかりがすべて AR(2) を指す——これが次数同定の理想形です。

対比実験:MA(2)過程と白色雑音

判定表の残り 2 行も数値で確認します。MA(2) 過程 \(x_t = \varepsilon_t + 0.6\varepsilon_{t-1} + 0.4\varepsilon_{t-2}\) と純粋な白色雑音を生成し、同じ解析にかけます。

# MA(2): ACFがラグ2で切断するはず
np.random.seed(29)
theta = np.array([0.6, 0.4])
y = ArmaProcess(np.r_[1], np.r_[1, theta]).generate_sample(nsample=N)
acf_y = acf(y, nlags=10, fft=True)
pacf_y = pacf(y, nlags=10, method="ywadjusted")
print("MA(2) ACF  1-5:", np.round(acf_y[1:6], 3))
print("MA(2) PACF 1-5:", np.round(pacf_y[1:6], 3))
print("MA(2) ACF 有意ラグ:", (np.where(np.abs(acf_y[1:]) > ci)[0] + 1).tolist())

# 白色雑音: どちらも全ラグで非有意のはず
np.random.seed(2)
w = np.random.randn(N)
acf_w = acf(w, nlags=10, fft=True)
pacf_w = pacf(w, nlags=10, method="ywadjusted")
print("WN ACF  1-5:", np.round(acf_w[1:6], 3))
print("WN PACF 1-5:", np.round(pacf_w[1:6], 3))
print("WN 有意ラグ(ACF/PACF):",
      (np.where(np.abs(acf_w[1:]) > ci)[0] + 1).tolist(),
      (np.where(np.abs(pacf_w[1:]) > ci)[0] + 1).tolist())

実行結果:

MA(2) ACF  1-5: [ 0.598  0.31   0.065  0.003 -0.01 ]
MA(2) PACF 1-5: [ 0.599 -0.074 -0.142  0.065  0.009]
MA(2) ACF 有意ラグ: [1, 2]
WN ACF  1-5: [-0.034  0.002  0.039 -0.027  0.039]
WN PACF 1-5: [-0.034  0.001  0.04  -0.025  0.038]
WN 有意ラグ(ACF/PACF): [] []

AR(2) とちょうど鏡写しの結果になりました。MA(2) では ACF がラグ 2 で切断(有意ラグは \(\{1, 2\}\) のみ、\(\hat{\rho}_3 = 0.065\) 以降は区間内)し、PACF の側は \(0.599 \to -0.074 \to -0.142\) と切断せずに減衰します。理論値 \(\rho_1 = 0.553\) , \(\rho_2 = 0.263\) に対して標本値 \(0.598\) , \(0.310\) も妥当な範囲です。白色雑音では ACF・PACF とも全ラグが信頼区間内に収まり、モデル化すべき構造がないことが分かります。

実務上の注意点

  • 定常性が前提:PACF・Yule-Walker はいずれも定常過程を仮定します。トレンドや単位根がある系列は、まずARIMAの枠組みで差分を取って(\(d\) を決めて)から \(p, q\) を同定します
  • 切断は「統計的な」切断:約 5% のラグは偶然有意になります。実際、本記事の AR(2) 実験でもラグを 20 まで見るとラグ 19 が偶然区間を超えました。低ラグからの連続した有意性を重視し、高ラグの単発は無視します
  • ARMA混合は判定表だけでは決まらない:ACF・PACF がどちらも減衰する場合は候補を絞る程度にとどめ、arma_order_select_ic のような情報量規準のグリッド探索に委ねます
  • 同定後の残差検定を忘れない:選んだモデルの残差が白色雑音になっているか(残差の ACF が全ラグ非有意か、Ljung-Box 検定が棄却されないか)の確認までがワンセットです。残差に構造が残る場合の検出は時系列異常検知の残差分析と同じ考え方です
  • 線形モデルの限界:PACF が捉えるのはあくまで線形の依存構造です。非線形な依存が支配的な系列では、LSTMのような非線形モデルが必要になります。とはいえその場合でも、ACF/PACF はラグ特徴量の設計や入力ウィンドウ長の目安として有効で、機械学習×時系列ガイドで整理したパイプラインの入り口に位置します

まとめ

  • **偏自己相関(PACF)**は中間ラグの線形効果を除去した相関であり、AR(k) 回帰の最終係数 \(\phi_{kk}\) と一致する
  • Yule-Walker 方程式(テプリッツ系)を次数を上げながら解けば PACF が得られ、Levinson-Durbin 再帰なら \(O(p^2)\) で全次数の係数・PACF・イノベーション分散が一度に求まる
  • 次数同定の判定表:AR(p) は PACF がラグ \(p\) で切断、MA(q) は ACF がラグ \(q\) で切断、ARMA は両方減衰
  • AR(2) の数値実験では、PACF の切断(\(\hat{\alpha}(2) = -0.312\) が有意、\(\hat{\alpha}(3)\) 以降は \(\pm 0.088\) 内)、ar_select_order の AIC/BIC(ともに \(p=2\) )、arma_order_select_ic(AIC/BIC とも \((2,0)\) )がすべて一致した
  • MA(2) では逆に ACF が切断(有意ラグ \(\{1,2\}\) のみ)、白色雑音では両者とも全ラグ非有意——判定表の全行を数値で確認できた
  • 判定表と情報量規準は競合ではなく相補的。視覚的判定で候補を絞り、AIC/BIC で客観的に決め、残差の白色性検定で締める

おすすめ書籍

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

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

関連記事

参考文献

  • Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time Series Analysis: Forecasting and Control (5th ed.). Wiley.
  • Brockwell, P. J., & Davis, R. A. (2016). Introduction to Time Series and Forecasting (3rd ed.). Springer.
  • Durbin, J. (1960). “The fitting of time-series models.” Revue de l’Institut International de Statistique, 28(3), 233-244.
  • Akaike, H. (1974). “A new look at the statistical model identification.” IEEE Transactions on Automatic Control, 19(6), 716-723.
  • statsmodels: statsmodels.tsa.stattools.pacf
  • statsmodels: statsmodels.tsa.ar_model.AutoReg