昨日の気温が高ければ、今日の気温も高い確率が大きい。先週ぐっと売れた商品は、今週もまだ余熱で売れているかもしれない。先月の株価が大きく動いたなら、今月の値動きにもその記憶がうっすら残っているかもしれない。私たちのまわりにある時系列データは、独立した点の寄せ集めではなく、過去が現在に影響を残しながら進んでいきます。
この「過去と現在のつながり」を1つの数値で定量化するのが、自己共分散(auto-covariance) と、それを正規化した 自己相関係数(auto-correlation function, ACF) です。自己共分散は、いま観測した値とラグ $k$ だけ前の値との「足並みの揃い方」を測る統計量で、時系列解析のあらゆる手法の出発点に座っています。
応用範囲はとても広く、たとえば次のような場面でフル活用されます。
- 金融: 株価や為替リターンに「ボラティリティクラスタリング」と呼ばれる相関構造があるか、ACFで確認する。アルゴリズムトレード戦略の基礎情報になります。
- 気象・気候: 気温・降水量・海面水温の周期成分(日周期・季節周期・ENSO周期)を、ACFのスパイクから検出する。
- 産業センサ/IoT: 振動・温度ログから機械の正常な周期パターンを学習し、ACFが崩れたら異常検知のトリガにする。
- 計量経済学: ARMA・ARIMA・SARIMAなどの古典モデルでは、ACFとPACFの形を見て次数 $(p, q)$ を決める「モデル同定」が出発点です。
- 信号処理: ノイズ除去・予測・推定の最適性は、信号の自己相関構造に依存します。ウィーナーフィルタなどはACFをそのまま使って設計されます。
本記事の内容
- 自己共分散・自己相関の直感と数学的定義 (なぜ「自分自身との共分散」なのか)
- 偏自己相関(PACF)で「途中の媒介を除いた直接の相関」を取り出す方法
- 弱定常性とエルゴード性 — 自己共分散が意味を持つための前提
- 周期性検出: ACFのスパイクから季節成分・周期成分を読み取る
- AR/MA/ARMAモデルとの関係 — ACF/PACFの形でモデルを同定する
- Ljung-Box検定で「相関が残っていないか」を統計的に検定する
- numpy + statsmodels によるPython実装と、シミュレーション・実データへの応用
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
通常の2変数の共分散 $\mathrm{Cov}(X, Y)$ がしっかり頭に入っていると、自己共分散は「同じ系列の過去と未来をペアにした共分散」だと素直に受け入れられます。
自己共分散の直感 — 「自分の過去とペアを組む共分散」
自己共分散の数式に行く前に、まずイメージを掴みましょう。通常の共分散 $\mathrm{Cov}(X, Y)$ は、別々の確率変数 $X$ と $Y$ がどれだけ一緒に動くかを測る量でした。たとえば「気温と冷たい飲料の売上」、「広告費と売上」のように、2つの違う変数を並べてペアを作りました。
ところが時系列では、わざわざ別の変数を持ち出さなくても「ペア」が作れます。同じ系列 $y_t$ を、自分自身からラグ $k$ だけずらしたコピー $y_{t-k}$ と並べればよいのです。すると、
- $y_t$ と $y_{t-1}$ のペア (ラグ1)
- $y_t$ と $y_{t-2}$ のペア (ラグ2)
- …
という、自分自身を「過去にずらしたコピー」とのペアが無限に作れます。このずれたペアの共分散こそが自己共分散です。「auto-(自分自身)」という接頭辞は、ここから来ています。
ペアを作って共分散を取るというこの操作の意味は、後で典型的な4系列のACFを並べた図 (本記事の後ろ) を見るとさらに鮮明になります。まずは直感を押さえましょう。
- ラグ1の自己共分散が大きい: 昨日が高ければ今日も高い、昨日が低ければ今日も低い。時系列が「滑らか」で、急激には変わらない。
- ラグ1がほぼゼロ: 過去の値と今の値は無関係 (例: 公正なコイン投げの記録)。
- ラグ1が負: 昨日が高ければ今日は低い、昨日が低ければ今日は高い、というジグザグした変動。
これだけでも「時系列の性格」が大雑把に読み取れます。ラグを変えながら自己共分散を並べたものを 自己共分散関数 (autocovariance function) と呼び、時系列の構造を測る基本指標になります。
ここまでで「自己共分散がペアの共分散である」というイメージが掴めました。次に、これを数式できちんと書き下しましょう。
自己共分散の定義
時系列 $\{y_t\}_{t=1}^{T}$ に対する 標本自己共分散 (sample autocovariance) は、ラグ $k$ について次式で定義されます。
$$ \begin{equation} \hat{\gamma}_k = \frac{1}{T} \sum_{t=k+1}^{T} (y_t – \bar{y})(y_{t-k} – \bar{y}) \end{equation} $$
ここで、
- $k \in \{0, 1, 2, \dots\}$ は ラグ (lag): 何ステップ遡るかを表す非負整数。
- $\bar{y} = \frac{1}{T}\sum_{t=1}^{T} y_t$ は時系列全体の標本平均。
- $T$ は系列の長さ (観測点数)。
通常の共分散 $\frac{1}{n}\sum_i (x_i – \bar{x})(y_i – \bar{y})$ と比べると、$x$ と $y$ がそれぞれ「いまの値」と「$k$ ステップ前の値」に置き換わっただけだと分かります。
$k=0$ のときは分散になる
ラグ0、つまり $k = 0$ の場合を $(1)$ 式に代入すると、
$$ \hat{\gamma}_0 = \frac{1}{T}\sum_{t=1}^{T} (y_t – \bar{y})^2 $$
となり、これはまさに $y_t$ の標本分散そのものです。自己共分散はラグ0で分散と一致し、ラグが増えるにつれて「時間的に離れた値の連動」を測るように一般化された量だと考えられます。
分母を $T$ にする理由
(1) 式の分母は $T – k$ ではなく $T$ になっています。データのペア数は $T – k$ 個しかないのに、なぜ $T$ で割るのでしょうか? これは「ラグが大きいほど推定の信頼性が下がる」ことを陽に表現するための、バイアスを犠牲にして分散を抑える設計です。
- 分母を $T – k$ にする方式 (不偏推定量寄り) → 大きなラグで値が暴れやすい
- 分母を $T$ にする方式 (本記事の式) → 大きなラグで値が0に向けて減衰し、推定が安定する
時系列解析のソフトウェア (statsmodels, R の acf など) はほぼ全て後者を採用しています。本記事もそれに合わせます。
計算のイメージ
数式だけだと味気ないので、簡単な例で手を動かしてみましょう。データ点 $y = [4, 2, 6, 1, 9, 7, 1, 2, 7, 6]$ を考えると、平均は $\bar{y} = 4.5$ です。
$k = 1$ の自己共分散を計算するなら、
$$ \hat{\gamma}_1 = \frac{1}{10}\sum_{t=2}^{10} (y_t – 4.5)(y_{t-1} – 4.5) $$
を計算します。式の意味するところは「いまの値からの偏差」と「1つ前の値からの偏差」を掛け合わせ、それを全時刻で足し合わせて平均するということです。
上の図では、青色の実線が時系列 $y_t$、赤色の線が平均 $\bar{y}$ で、$(y_t – \bar{y})(y_{t-k} – \bar{y})$ という積に対応する量が、図中で同じ色の領域同士の掛け算として表現されています。同色の組を全部掛けて足し、$T$ で割ったものが自己共分散の値です。
$k = 2$ では1つ多く遡るので、ペアの数が1つ減ります。
このように、ラグを変えながら $\hat{\gamma}_k$ を計算していくと、$\hat{\gamma}_0, \hat{\gamma}_1, \hat{\gamma}_2, \dots$ という数列ができます。これが自己共分散関数の標本推定値です。
ここまでで、自己共分散が「自分自身の過去とのペアで作った共分散」であることが定義式と図の両方で確認できました。しかし $\hat{\gamma}_k$ にはまだ不便な点があります。値の大きさがデータのスケールに依存するのです。1m単位で測った身長の時系列と、1mm単位で測ったものでは、共分散の値が $10^6$ 倍も違ってきます。系列を横断して比較したり、シンプルな尺度で「強い相関か弱い相関か」を判断したい場合は困ります。そこで次に、これを正規化した自己相関係数を導入します。
自己相関係数 (ACF)
定義
自己相関係数 (autocorrelation function, ACF) は、自己共分散をラグ0の自己共分散 (= 分散) で割って正規化したものです。
$$ \begin{equation} \hat{\rho}_k = \frac{\hat{\gamma}_k}{\hat{\gamma}_0} \end{equation} $$
通常の相関係数 $\mathrm{Corr}(X, Y) = \mathrm{Cov}(X, Y) / \sqrt{\mathrm{Var}(X)\mathrm{Var}(Y)}$ で、$X$ と $Y$ が同じ系列の異なるラグである場合、両者の分散はともに $\hat{\gamma}_0$ なので、分母が $\sqrt{\hat{\gamma}_0 \cdot \hat{\gamma}_0} = \hat{\gamma}_0$ となり、(2) 式が自然に出てきます。
ACFの性質
ACFには以下の重要な性質があります。
- $\hat{\rho}_0 = 1$: $k = 0$ では同じ値同士の比較なので、相関は必ず1。
- $|\hat{\rho}_k| \leq 1$: 通常の相関係数と同じく、$-1$ から $1$ の間に収まる。
- 対称性: 母集団のACFは $\rho_k = \rho_{-k}$ (定常な実数値系列では時刻の方向を入れ替えても同じ)。標本版でもラグの符号は無視できます。
- 正規化済み: スケールに依存しないので、別の時系列と直接比較できる。
ACFを縦軸に取り、横軸にラグ $k$ を並べた図を コレログラム (correlogram) と呼びます。慣れてくると、コレログラムを一目見るだけで「この系列は AR(1) っぽい」「これは周期12の季節性がある」と当たりが付くようになります。
95%信頼区間と「相関なし」の判定
白色雑音 (i.i.d. なノイズ) の標本ACFは、$T$ が十分大きいとき近似的に平均0・分散 $1/T$ の正規分布に従うことが知られています (Bartlett の公式)。そこから、
$$ \text{95% 信頼区間: } \quad \pm \frac{1.96}{\sqrt{T}} $$
の範囲に $\hat{\rho}_k$ が収まっていれば「ラグ $k$ で統計的に有意な相関は無い」と判断されます。コレログラムでよく見る灰色の帯はこの信頼区間です。
上の図は、白色雑音・AR(1)・MA(1)・AR(2) の4種類の系列について、それぞれの標本ACFとPACFを並べたものです。次の3点が読み取れます。
- 白色雑音は、ラグ0の1本だけが立ち、それ以降はすべて灰色の信頼区間に収まっている。
- AR(1) (係数0.7) のACFは $0.7^k$ で指数的に減衰し、長く尾を引く。一方でPACFはラグ1だけが立つ「鋭い切断 (cut-off)」。
- MA(1) (係数0.7) はちょうど逆で、ACFがラグ1で切れ、PACFが減衰する。
この「ACFは尾を引きPACFは切れる/ACFは切れPACFは尾を引く」のパターン認識が、後で ARMA モデルの次数を決めるときに直接効いてきます。
ここまでで、自己共分散と自己相関の数学的な定義と性質、コレログラムの読み方が掴めました。しかしACFには重大な弱点があります。ラグ $k$ のACF $\hat{\rho}_k$ は、間にある $y_{t-1}, y_{t-2}, \dots, y_{t-k+1}$ を経由した間接的な伝播まで含めて測ってしまうのです。「$y_t$ と $y_{t-k}$ の直接の関係だけを取り出したい」と思ったときに困ります。それを解決するのが、次に紹介する偏自己相関 (PACF) です。
偏自己相関 (PACF) — 媒介を除いた直接の相関
モチベーション
AR(1) モデル $y_t = \phi y_{t-1} + \varepsilon_t$ を考えると、$y_t$ は $y_{t-1}$ にしか直接依存しません。ところが $y_{t-1}$ 自身が $y_{t-2}$ に依存しているので、$y_t$ と $y_{t-2}$ の間にも見かけの相関が生じます。
$$ y_t = \phi(\phi y_{t-2} + \varepsilon_{t-1}) + \varepsilon_t = \phi^2 y_{t-2} + \phi \varepsilon_{t-1} + \varepsilon_t $$
なので、ラグ2のACFは $\phi^2$ という形で表れます。ですが直接 $y_{t-2}$ に依存しているわけではありません。媒介変数 $y_{t-1}$ を経由した「お下がり」の相関です。
「$y_t$ と $y_{t-k}$ の関係から、間の $y_{t-1}, \dots, y_{t-k+1}$ の影響を取り除いた、直接の相関」を測るのが偏自己相関 (partial autocorrelation function, PACF) です。
定義
PACF は、$y_t$ と $y_{t-k}$ をそれぞれ $\{y_{t-1}, \dots, y_{t-k+1}\}$ で線形回帰した残差同士の相関係数として定義されます。形式的には、ラグ $k$ の PACF $\phi_{kk}$ は、AR($k$) モデル
$$ y_t = \phi_{k1} y_{t-1} + \phi_{k2} y_{t-2} + \dots + \phi_{kk} y_{t-k} + \varepsilon_t $$
を最小二乗法で当てはめたときの、最大ラグの係数 $\phi_{kk}$ として与えられます。実用的には、Yule-Walker方程式を Durbin-Levinson 漸化式で解くことで、ACF $\{\hat{\rho}_k\}$ から効率的に PACF $\{\phi_{kk}\}$ を計算できます。
ACFとPACFの「対の役割」
ACFとPACFは時系列モデル同定の道具立てとして対をなしています。
| モデル | ACF | PACF |
|---|---|---|
| AR($p$) | 指数 / 正弦的に減衰 | ラグ $p$ で切れる |
| MA($q$) | ラグ $q$ で切れる | 指数 / 正弦的に減衰 |
| ARMA($p, q$) | 減衰 | 減衰 |
つまり、
- PACF が切れたラグ = AR の次数 $p$ の候補
- ACF が切れたラグ = MA の次数 $q$ の候補
このパターン認識を Box-Jenkins のモデル同定 と呼びます。1976年以来、伝統的な時系列解析の中核を担ってきた方法論で、現代でも統計コンサルや量的金融のスタンダードな出発点になっています。
ここまででACFとPACFという「時系列の構造を測る2つの基本指標」が揃いました。しかし、これらが意味のある量になるためには、時系列がある種の安定性を持っている必要があります。具体的には、平均や分散、自己共分散が時刻の関数になっていない、つまり「時系列の統計的性質が時間で変わらない」という性質です。この性質を定常性と言います。
定常性とエルゴード性 — 自己共分散が意味を持つ前提
弱定常性の定義
実用上、最もよく使われるのが 弱定常性 (weak / covariance stationarity) です。次の3条件を満たすとき、系列 $\{y_t\}$ は弱定常であると言います。
- 平均が一定: $\mathbb{E}[y_t] = \mu$ (時刻 $t$ に依存しない)
- 分散が有限かつ一定: $\mathrm{Var}(y_t) = \sigma^2 < \infty$
- 自己共分散が時刻差のみに依存: $\mathrm{Cov}(y_t, y_{t-k}) = \gamma_k$ (時刻 $t$ には依存しない)
条件3が決定的に重要です。弱定常な系列では、「ラグ $k$ の自己共分散」は時刻 $t$ がどこであっても同じ値を取るので、1つの数値 $\gamma_k$ で代表させることができます。これがあるから、(1) 式のような標本平均で $\gamma_k$ を推定する操作が意味を持つのです。
なぜ定常性が必要か
非定常な系列(平均がトレンドで動いていたり、分散が時間とともに膨らんだりするもの)に対して (1) 式を機械的に計算しても、得られる $\hat{\gamma}_k$ は「時刻と無関係な真の $\gamma_k$」を推定したものにはなりません。むしろ平均が動いているせいで、$y_t – \bar{y}$ がほとんど常に同じ符号になり、$\hat{\gamma}_k$ が緩やかに減衰するだけのほぼ無意味なグラフが出てきてしまいます。
上の図には、定常な AR(1) と、3種類の非定常系列 (トレンドあり / ランダムウォーク / 分散変動) を並べました。次の特徴が一目で見て取れます。
- 定常な AR(1) では、移動平均 (太い黒線) が水平に近く、振幅も時刻に依らずほぼ一定。
- トレンドありの系列では、移動平均が右肩上がりに動き、平均が時刻の関数になっている。
- ランダムウォークでは、ある瞬間にはトレンドに見える挙動が、別の時期には逆方向に動く。一般に分散が $t$ に比例して発散する。
- 分散変動の系列では、後半で振幅が明らかに大きくなり、分散が時刻で変わっている。
トレンドがある場合は、差分を取る ($\Delta y_t = y_t – y_{t-1}$) ことで定常化することが多く、これが ARIMA モデルの「I (Integrated)」の正体です。分散変動の場合は対数変換や GARCH 系のモデルで対処します。
エルゴード性
弱定常性は「真の $\gamma_k$ が時刻に依らず存在する」ことを保証してくれますが、標本 $\hat{\gamma}_k$ を計算するときには、もう1つ別の前提が必要です。それがエルゴード性 (ergodicity) です。
弱定常性は「同じ確率実験を繰り返したアンサンブル平均」での話で、エルゴード性は「1本の長い実現値の時間平均が、アンサンブル平均に収束する」という性質です。私たちが手元に持つのはたいてい1本の系列なので、その時間平均で母集団の量 ($\mu, \gamma_k$) を推定できる根拠がエルゴード性です。
実用上は、弱定常 + 自己共分散が十分速く減衰する ($\sum_k |\gamma_k| < \infty$ など) ならエルゴード性が成り立つ十分条件になります。AR(1) で $|\phi| < 1$ ならOK、というのが典型例です。
ここまでで「自己共分散が意味を持つには弱定常性とエルゴード性が必要」という土台が固まりました。次に、自己相関の代表的な応用として、周期性 (季節性) の検出を見ていきます。日次データに潜む週・年の周期、機械の振動データに潜む回転周期など、ACFはこうした繰り返しパターンを浮き彫りにする強力な道具です。
周期性検出 — 季節成分はACFに何を残すか
周期性とACFの関係
周期 $P$ で振動する成分が時系列に含まれている場合、$y_t$ と $y_{t-P}$ は近い値を取りやすくなります。これは ACF にラグ $P$ で大きなスパイクとして現れます。さらにラグ $2P, 3P, \dots$ でも、減衰しながらスパイクが続きます。
たとえば月次の小売売上データなら、12カ月ごとの季節性 (クリスマス商戦、ボーナス商戦など) がラグ12, 24, 36, … に山として現れます。週次の電力消費データなら、ラグ7, 14, 21, … に強いピークが出ます。日次気温データなら、ラグ365, 730, … に年周期のスパイクが立ちます。
多重周期の検出例
上の図は、年周期 (365日) と週周期 (7日) の2つが重なり合った合成時系列に対する標本ACFです。上段は元の系列(2年分)、下段はラグ0~400のACFです。
次の構造が読み取れます。
- ラグ7, 14, 21, … の位置に小さなスパイクが等間隔で並んでおり、週周期の存在を示している。
- ラグ365 には大きな1本の山がぽつんと立ち、年周期の存在を示している。
- ACF全体は緩やかに減衰しており、トレンド成分も少しは含まれていることが示唆される。
このように、ACFは複数の周期が混在している系列でも、それぞれの周期を別々のラグに「分けて」表示してくれます。フーリエ変換 (周波数領域の分析) が「振動成分の振幅を周波数別に並べる」のと相補的な役割を果たします。
Ljung-Box検定との合わせ技
「目で見て周期がありそう」と「統計的に有意に相関がある」は別物です。後者を裏付けるのが、後述する Ljung-Box 検定です。複数のラグをまとめて検定することで、「ACFのどこかにシグナルが残っているか」を1つの p 値に集約できます。
ここまでで、周期成分がACFにどんなパターンを残すかを見てきました。次に、ACFとPACFの形を見て「データ生成過程のモデル」を当てる、Box-Jenkins モデル同定の中身を詳しく追っていきましょう。
AR / MA / ARMA モデルとの関係
AR (自己回帰) モデル
AR($p$) モデルは、現在の値を過去 $p$ 個の自分自身の線形結合で表現する確率過程です。
$$ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p} + \varepsilon_t, \quad \varepsilon_t \sim \mathrm{WN}(0, \sigma^2) $$
ここで $\varepsilon_t$ は平均0・分散 $\sigma^2$ の白色雑音 (互いに無相関) です。
定常性のためには、特性方程式 $1 – \phi_1 z – \dots – \phi_p z^p = 0$ の根がすべて単位円の外にある必要があります。AR(1) なら単純に $|\phi_1| < 1$。
AR(1) の理論ACF
AR(1) モデル $y_t = \phi y_{t-1} + \varepsilon_t$ (定数項は0、$|\phi| < 1$) の自己共分散を求めてみましょう。両辺に $y_{t-k}$ を掛けて期待値を取ります。
$$ \mathbb{E}[y_t y_{t-k}] = \phi \mathbb{E}[y_{t-1} y_{t-k}] + \mathbb{E}[\varepsilon_t y_{t-k}] $$
$\varepsilon_t$ は $y_{t-k}$ ($k \geq 1$) と無相関なので右第2項は0。よって、
$$ \gamma_k = \phi \gamma_{k-1} \quad (k \geq 1) $$
という漸化式が得られます。ここから $\gamma_k = \phi^k \gamma_0$、すなわち
$$ \boxed{\rho_k = \phi^k} $$
がAR(1) の理論ACFです。指数関数的な減衰が AR(1) の指紋です。$\phi > 0$ なら全部正で滑らかに減衰、$\phi < 0$ なら符号が交互に入れ替わるジグザグな減衰になります。
MA (移動平均) モデル
MA($q$) モデルは、現在の値を過去 $q$ 個のノイズの線形結合で表現します。
$$ y_t = \mu + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \dots + \theta_q \varepsilon_{t-q} $$
MA は常に定常で、ACFはちょうどラグ $q$ で切れます(0になる)。MA(1) なら、
$$ \rho_k = \begin{cases} 1 & k = 0 \\ \theta_1 / (1 + \theta_1^2) & k = 1 \\ 0 & k \geq 2 \end{cases} $$
これも導出は同様で、$y_t y_{t-k}$ の期待値を取り、ノイズ項同士の独立性を使って計算するだけです。
ARMA / ARIMA / SARIMA
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} $$
- ARIMA($p, d, q$): 元系列を $d$ 階差分してから ARMA を当てはめる。トレンドのある非定常系列に対応。
- SARIMA($p, d, q$)($P, D, Q$)$_s$: 季節性 (周期 $s$) を別の AR/MA で扱う ARIMA の拡張。月次・週次データの定番。
Box-Jenkins 手順のまとめ
- 可視化: 時系列をプロットして、トレンド・季節性・異常値を目視。
- 定常化: 必要なら差分・対数変換で定常な系列にする。
- ACF/PACFのプロット: コレログラムから $(p, q)$ の候補を絞り込む。
- 当てはめ: 候補の $(p, q)$ で最尤推定。
- 残差診断: 残差が白色雑音か Ljung-Box などで検定。OK でなければ次数を変えて再挑戦。
- 予測: 通った最終モデルで先の値を予測する。
ここまでで、自己共分散・ACF・PACF が ARMA 系モデルの設計に直接結びついていることが分かりました。続いて、ステップ5の「残差が白色雑音かどうか」を統計的に判定する道具、Ljung-Box 検定を見ていきます。
Ljung-Box 検定 — 「相関が残っていないか」を一度に検定する
検定統計量
各ラグの $\hat{\rho}_k$ が個別に信頼区間に収まっているかを見るだけでは、「複数のラグで小さな相関がまとめて残っている」状況を見逃しがちです。そこで、最初の $h$ 個のラグをまとめて検定する統計量が Ljung-Box (1978) によって提案されました。
$$ \begin{equation} Q(h) = T(T + 2) \sum_{k=1}^{h} \frac{\hat{\rho}_k^2}{T – k} \end{equation} $$
帰無仮説 $H_0$ は「$\rho_1 = \rho_2 = \dots = \rho_h = 0$」、つまり「最初の $h$ ラグまで自己相関は無い (白色雑音)」です。$H_0$ のもとで、$Q(h)$ は漸近的に自由度 $h$ のカイ二乗分布 $\chi^2_h$ に従います。よって、
- $Q(h) > \chi^2_{h, 1-\alpha}$ ならば、有意水準 $\alpha$ で $H_0$ を棄却
- そうでなければ「白色雑音と矛盾しない」と判断
ARMA モデルの残差に対して適用する場合は、推定したパラメータ数 $p + q$ を引いて自由度を $h – p – q$ にします。
振る舞いの可視化
上の図は、白色雑音と AR(1) ($\phi = 0.6$) の2系列について、$h = 1, \dots, 20$ の範囲で $Q(h)$ をプロットしたものです。
- 白色雑音 (緑) の $Q(h)$ は自由度 $h$ のカイ二乗の上側5%臨界値 (灰破線) を一度も超えず、帰無仮説は棄却されない。
- AR(1) (赤) の $Q(h)$ は最初のラグから既に臨界値を大きく超え、「白色雑音」という帰無仮説は強く棄却される。
このように Ljung-Box 検定は、「ACFを目で見て判断する」を統計的な p 値に変える重要な道具です。とくに ARMA モデルを当てはめた後の残差診断で必ず使われます。
Box-Pierce 検定との違い
似た統計量に Box-Pierce 検定 $Q^{*}(h) = T \sum_{k=1}^h \hat{\rho}_k^2$ がありますが、これは小標本で漸近分布への近似が悪いため、Ljung-Box の係数 $T(T+2)/(T-k)$ で補正したものが事実上の標準です。
ここまでで、時系列に相関が残っているかを統計的に検定する道具まで揃いました。一方、私たちは標本ACFの値を「真の ACF の推定値」として使ってきましたが、標本サイズが小さいとき推定値はどれだけ信頼できるのでしょうか? 次に、シミュレーションでその挙動を確認します。
標本ACFは理論値に収束するか — シミュレーション
ACFの定義から、$\hat{\rho}_k$ は母集団の $\rho_k$ の標本推定量です。エルゴード性が成り立つなら $T \to \infty$ で $\hat{\rho}_k \to \rho_k$ に収束するはずですが、有限の $T$ ではどれだけ理論値に近いのでしょう?
AR(1) ($\phi = 0.75$) を例にとり、$T = 50, 200, 1000, 10000$ の4通りで標本ACFを計算し、理論値 $\phi^k = 0.75^k$ と比較してみます。
このシミュレーションから次の3点が読み取れます。
- $T = 50$ では、標本ACFは理論値からそれなりにバラついており、ラグ10付近では負の値さえ出ている。短い系列ではACFの推定は不確実。
- $T = 200$ で大幅に改善するが、大きなラグでは依然として理論値からずれが残る。
- $T = 1000, 10000$ になると、標本ACFは理論値の指数曲線にきれいに乗り、エルゴード性による収束が視覚的に確認できる。
このことは実務上、短いデータでACFの細かい構造を信用しすぎてはいけないという強い警告でもあります。経験則として、各ラグについて少なくとも50点、できれば100点以上のペアが取れるよう $T \gg k_{\max}$ にしておくのが目安です。
理論との一致が確認できたので、次に「多次元時系列」の場合に拡張します。1つの系列の自分との相関ではなく、2つの系列のお互いとの相関を測る量が 交差相関 (cross-correlation) です。
多次元時系列の交差相関 (CCF)
定義
2つの時系列 $\{x_t\}$, $\{y_t\}$ に対し、ラグ $k$ の交差共分散と交差相関は次のように定義されます。
$$ \gamma_{xy}(k) = \mathrm{Cov}(x_t, y_{t+k}), \quad \rho_{xy}(k) = \frac{\gamma_{xy}(k)}{\sqrt{\gamma_{xx}(0)\,\gamma_{yy}(0)}} $$
自己共分散と違う点は、
- ラグの符号に意味がある: $\rho_{xy}(k)$ と $\rho_{xy}(-k)$ は別物。
- 正のラグは $x$ が $y$ より早く動く (リーディング) ことに対応、負のラグは逆。
リーディング・インジケータの検出
たとえば「広告費が売上に2週間遅れて効く」「PMI (景気先行指数) が GDP に1四半期先行する」のような関係は、CCFのピーク位置として現れます。
上の図では、先行系列 $x_t$ (青) を5期遅らせてノイズを足したものを追従系列 $y_t$ (赤) とし、CCF を計算しました。下段のCCFを見ると、ラグ $k = 5$ に明瞭なピークが立っています。この位置の検出が、リーディング・インジケータ分析や、入力と出力の時間遅れを持つシステム同定の基礎になります。
注意点 — 偽相関
CCFは便利ですが、両系列が独立に非定常 (例: トレンドあり) の場合、見かけ上の高い相関が出てしまう「偽相関 (spurious correlation)」が発生しやすい点には注意が必要です。CCFを取る前に、まず差分や対数差分で定常化することが推奨されます。
ここまでで、自己相関を多次元に拡張する道筋まで掴めました。最後に、ラグごとの相関を散布図として直接見てみる、シンプルだけど直感的な可視化を紹介します。
ラグ別の散布図 (Lag plot)
ACFが「全ラグの相関を1枚の棒グラフにまとめたもの」なら、ラグ散布図は「特定のラグの相関を、生の散布図で見せる」ものです。
上の図は AR(1) ($\phi = 0.7$) の系列に対して、ラグ1, 2, 5, 10 の散布図を並べたものです。横軸が $y_t$、縦軸が $y_{t+k}$。次のことが見て取れます。
- ラグ1 では強い右肩上がりの関係 ($r_1 \approx 0.69$)。直前の値が次の値をよく予測する。
- ラグ2 ($r_2 \approx 0.49$)、ラグ5 ($r_5 \approx 0.17$) と進むにつれ、点群が円形に広がり相関が弱まる。
- ラグ10 ではほぼ無相関。
ACFの値 $\hat{\rho}_k$ は、こうした散布図の傾き (回帰直線の傾き) を正規化したものに他なりません。コレログラムの数値だけでなく、たまにはラグ散布図で「生の関係」を覗いてみると、外れ値・非線形性・分布の歪みなど、ACFには現れない情報も拾えます。
ここまでで概念と可視化の全ツールが揃いました。最後に、これらを Python で実装していきましょう。
Python実装
自前で計算してみる (numpy のみ)
まずは定義式 (1) を素直に実装した自己共分散と、(2) を実装した自己相関係数を書いてみましょう。
import numpy as np
import matplotlib.pyplot as plt
def auto_covariance(y, max_lag):
"""標本自己共分散 gamma_k を k=0..max_lag まで返す."""
y = np.asarray(y, dtype=float)
T = len(y)
mean = y.mean()
centered = y - mean
gamma = np.empty(max_lag + 1)
for k in range(max_lag + 1):
gamma[k] = np.dot(centered[k:], centered[:T - k]) / T
return gamma
def auto_correlation(y, max_lag):
"""標本自己相関係数 rho_k を返す."""
gamma = auto_covariance(y, max_lag)
return gamma / gamma[0]
ポイントは2つあります。1つ目は、np.dot(centered[k:], centered[:T-k]) で内積を使ってループを潰している点。for t in range(k+1, T+1) を素朴に書くより数十倍速くなります。2つ目は、分母を T 固定にしている点で、これにより大きなラグでも値が安定して0に向かいます。
ダミー時系列でテスト
サイン波にノイズを乗せた合成系列を作り、ACFを描いてみます。
np.random.seed(5656)
x = np.linspace(0, 50, 50)
y = np.sin(x) + np.random.normal(0, 1, size=len(x)) + 10
fig, ax = plt.subplots(figsize=(10, 2), dpi=100)
ax.plot(x, y, color="red", linestyle="-")
ax.set_xlabel("時刻 t"); ax.set_ylabel("y(t)")
plt.show()
生成した時系列はおおよそ周期 $2\pi$ で振動しており、これにガウスノイズが重畳しています。振動の山と山がほぼ周期的に並んでいるので、ラグが $2\pi$ の整数倍のところでACFが大きくなりそうだ、と予想できます。
gamma = auto_covariance(y, max_lag=49)
result = gamma[1:] # k=1..49
fig, ax = plt.subplots(figsize=(10, 2), dpi=100)
ax.set_xlabel("ラグ k")
ax.set_xticks(range(1, 50), minor=True)
ax.stem(range(1, 50), result)
plt.show()
このコレログラムから、$2\pi \approx 6.28$ の整数倍 (約6, 12, 19, 25, …) のラグで自己共分散が大きな正の値を取っているのが読み取れます。これは元のサイン波の周期性が忠実に検出されているということ。逆に $\pi$ の奇数倍 ($\approx 3, 9, \dots$) のラグでは負の値を取り、「半周期ずらすと逆位相になる」性質も再現されています。
statsmodels で本格的なACF/PACF
実務では statsmodels の関数を使うと、信頼区間・PACF・Ljung-Box検定までまとめて扱えます。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
# AR(1) 系列を生成
np.random.seed(0)
T = 1000
phi = 0.7
eps = np.random.normal(0, 1, T + 100)
y = np.zeros(T + 100)
for t in range(1, T + 100):
y[t] = phi * y[t - 1] + eps[t]
y = y[100:] # 過渡部分を捨てる
# ACF / PACF を計算
acf_vals = acf(y, nlags=30, fft=True)
pacf_vals = pacf(y, nlags=30, method="ywm")
print(f"理論ACF (phi^k): {phi**1:.3f}, {phi**2:.3f}, {phi**3:.3f}")
print(f"標本ACF: {acf_vals[1]:.3f}, {acf_vals[2]:.3f}, {acf_vals[3]:.3f}")
print(f"標本PACF: {pacf_vals[1]:.3f}, {pacf_vals[2]:.3f}, {pacf_vals[3]:.3f}")
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(y, lags=30, ax=axes[0])
plot_pacf(y, lags=30, ax=axes[1], method="ywm")
plt.tight_layout()
plt.show()
出力例: – 理論ACF (phi^k): 0.700, 0.490, 0.343 – 標本ACF: 0.696, 0.481, 0.330 – 標本PACF: 0.696, -0.005, 0.014
理論値 $\phi^k = 0.7^k$ と標本ACFがほぼ一致しているのが確認できます。PACF はラグ1だけが大きく、ラグ2以降はほぼゼロ — AR(1) の指紋通りです。実務でも、plot_acf / plot_pacf を眺めて「PACFがラグ1で切れているから AR(1) で試そう」と次数を決めるのが定石です。
Ljung-Box検定の実行
import numpy as np
from statsmodels.stats.diagnostic import acorr_ljungbox
np.random.seed(1)
wn = np.random.normal(0, 1, 500) # 白色雑音
ar1 = np.zeros(500)
eps = np.random.normal(0, 1, 500)
for t in range(1, 500):
ar1[t] = 0.6 * ar1[t - 1] + eps[t] # AR(1)
# h = 1, 5, 10, 20 で検定
for h in [1, 5, 10, 20]:
res_wn = acorr_ljungbox(wn, lags=[h], return_df=True)
res_ar = acorr_ljungbox(ar1, lags=[h], return_df=True)
print(f"h={h:2d} 白色雑音: Q={res_wn['lb_stat'].iloc[0]:7.2f}, "
f"p={res_wn['lb_pvalue'].iloc[0]:.3f} | "
f"AR(1): Q={res_ar['lb_stat'].iloc[0]:7.2f}, "
f"p={res_ar['lb_pvalue'].iloc[0]:.3f}")
実行結果例:
h= 1 白色雑音: Q= 0.62, p=0.430 | AR(1): Q= 183.71, p=0.000
h= 5 白色雑音: Q= 3.94, p=0.558 | AR(1): Q= 267.34, p=0.000
h=10 白色雑音: Q= 8.95, p=0.537 | AR(1): Q= 287.16, p=0.000
h=20 白色雑音: Q= 18.66, p=0.543 | AR(1): Q= 295.41, p=0.000
白色雑音では全ての $h$ で p 値が大きく (0.4 ~ 0.6)、帰無仮説 $H_0$「相関なし」は棄却されません。AR(1) ではどの $h$ でも p 値が事実上0になり、強く棄却されます。実務的には、ARMA を当てはめた後の残差に対して同じ検定を行い、p > 0.05 ならモデルは適合していると判断する流れになります。
周期性検出 — 月次小売売上に応用
実データを使った例として、月次データに対する季節性検出を見てみましょう。シミュレーションですが、現実の月次小売データに近い特性を持たせています。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
# 12カ月周期 + トレンド + ノイズの月次データ
np.random.seed(42)
T = 120 # 10年分
t = np.arange(T)
trend = 0.5 * t
seasonal = 30 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 5, T)
sales = 100 + trend + seasonal + noise
# 差分して定常化
diff = np.diff(sales)
# 差分系列のACF
acf_vals = acf(diff, nlags=36, fft=True)
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
axes[0].plot(t, sales, color="steelblue", lw=1.2)
axes[0].set_title("月次売上 (合成データ)"); axes[0].set_xlabel("月"); axes[0].set_ylabel("売上")
axes[0].grid(alpha=0.3)
axes[1].stem(range(37), acf_vals)
axes[1].axhline(1.96/np.sqrt(len(diff)), color="gray", ls="--")
axes[1].axhline(-1.96/np.sqrt(len(diff)), color="gray", ls="--")
axes[1].set_title("1階差分系列の標本ACF — ラグ12, 24に季節性スパイク")
axes[1].set_xlabel("ラグ"); axes[1].set_ylabel("自己相関")
axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
このグラフでは、元の売上に右肩上がりのトレンドと12カ月周期の波がはっきり見えます。差分系列のACFでは、トレンドが除かれた一方でラグ12と24に明確なスパイクが残っており、年次の季節性は差分では消えないことが分かります。実務では、これを見て「季節差分も併用する SARIMA で試そう」とモデル選択が進みます。
株価リターンへの応用
最後に、株価のような金融データへの応用も触れておきましょう。一般に、株価そのものはランダムウォークに近いので非定常ですが、対数差分 (log-return) を取ると弱定常になります。さらに、log-return 自身のACFはほぼ無相関 (効率的市場仮説と整合) ですが、二乗 log-return のACFには強い相関 (ボラティリティクラスタリング) が残ることが知られています。
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
# GARCH(1,1) っぽい合成データ
np.random.seed(7)
T = 2000
omega, alpha, beta = 0.05, 0.1, 0.85
sigma2 = np.zeros(T)
r = np.zeros(T)
sigma2[0] = omega / (1 - alpha - beta)
for t in range(1, T):
sigma2[t] = omega + alpha * r[t-1]**2 + beta * sigma2[t-1]
r[t] = np.sqrt(sigma2[t]) * np.random.normal()
acf_r = acf(r, nlags=40, fft=True)
acf_r2 = acf(r**2, nlags=40, fft=True)
fig, axes = plt.subplots(2, 1, figsize=(12, 6))
axes[0].stem(range(41), acf_r)
axes[0].axhline(1.96/np.sqrt(T), color="gray", ls="--")
axes[0].axhline(-1.96/np.sqrt(T), color="gray", ls="--")
axes[0].set_title("log-return r_t のACF — ほぼ無相関")
axes[0].set_xlabel("ラグ"); axes[0].set_ylabel("ACF"); axes[0].grid(alpha=0.3)
axes[1].stem(range(41), acf_r2)
axes[1].axhline(1.96/np.sqrt(T), color="gray", ls="--")
axes[1].axhline(-1.96/np.sqrt(T), color="gray", ls="--")
axes[1].set_title("二乗 log-return r_t^2 のACF — 強い相関 (ボラクラスタリング)")
axes[1].set_xlabel("ラグ"); axes[1].set_ylabel("ACF"); axes[1].grid(alpha=0.3)
plt.tight_layout(); plt.show()
この2つのコレログラムから、金融時系列の重要な性質が読み取れます。リターンそのもののACFは信頼区間内に収まる「白色雑音っぽい」振る舞いを見せる一方、二乗リターンには長く尾を引く強い相関が残っています。これは「大きな変動の後には大きな変動が続きやすい」というボラティリティクラスタリングの定量的表現です。
この「平均はランダムだが分散は予測可能」という構造を捉えるのが GARCH 系のモデルで、自己共分散 → ACF → 二乗系列のACF → GARCH という流れは、量的金融の基本ルートになっています。
自己共分散の応用領域まとめ
最後に、実際にどんな分野で自己共分散が使われているかを駆け足でまとめます。
| 分野 | 使われ方 |
|---|---|
| 計量経済学 | GDP、消費者物価指数、失業率などのマクロ時系列の ARIMA モデリング |
| 量的金融 | リターンのACF (効率的市場仮説の検証)、二乗リターンのACF (ボラティリティモデリング) |
| 信号処理 | ウィーナーフィルタ、最尤検出器、レーダーの整合フィルタの設計 |
| 通信工学 | チャネル推定、雑音の色 (white / colored noise) の同定 |
| 制御工学 | システム同定、入出力相互相関による伝達関数推定 |
| 気象学 | 季節成分・ENSO・年々変動などの周期検出 |
| 産業 IoT | 機械振動からの異常検知、回転周期の同定 |
| 生物学 | 脳波 (EEG)・心拍 (HRV) の周期解析、概日リズム |
| 天文学 | 周期的変光星の検出、ガンマ線バーストの統計 |
| 地震学 | 余震系列のクラスタリング、地震活動の時間相関 |
このように、自己共分散と自己相関は時系列が現れるすべての分野で基本ツールとして使われています。先のステップとしては、ARMAモデルの最尤推定、状態空間モデル (カルマンフィルタ)、スペクトル解析 (ピリオドグラム) などへ進むと、自己共分散と表裏一体の周波数領域の見方が手に入ります。
まとめ
本記事では、時系列解析の基本指標である自己共分散と自己相関について、定義・直感・実装・応用までを通して解説しました。
- 自己共分散 $\gamma_k$ は「自分自身をラグ $k$ ずらしたコピーとの共分散」。ラグ0で分散と一致する。
- 自己相関係数 $\rho_k$ は自己共分散を $\gamma_0$ で正規化したもの。スケール不変で系列横断比較に便利。
- 偏自己相関 $\phi_{kk}$ は、媒介ラグを除いた「$y_t$ と $y_{t-k}$ の直接の相関」。AR の次数同定に使う。
- 弱定常性とエルゴード性 が、$\hat{\gamma}_k$ を母集団 $\gamma_k$ の推定値として使う前提。
- ACF/PACFの形 からARMA次数を選ぶ Box-Jenkins 法は古典時系列の中核。
- Ljung-Box 検定 は複数ラグの相関をまとめて検定する道具。残差診断の必須項目。
- 周期成分 はACFのスパイクとして現れ、季節性検出のシグナルになる。
- 交差相関 で2変量時系列のリードラグ関係を測れる。
次のステップとして、以下の記事も参考にしてください。