広義定常過程とは?定常性の定義と検証を理解する

センサーから流れてくる電圧の時系列を100サンプル受け取ったとします。「このセンサーのノイズの標準偏差はいくつですか」と聞かれたら、あなたは迷わず100個の標準偏差を計算して答えるでしょう。では、株価の日次終値を100日ぶん受け取って「この株価の平均はいくつですか」と聞かれたら、同じように100日ぶんの平均を計算して答えていいのでしょうか。

答えは「センサーはよくて、株価はまずい」です。この2つを分けているものが定常性(stationarity)です。センサーのノイズは、昨日測ろうが今日測ろうが同じ確率法則から出てきます。だから100個のサンプルは同じ壺から引いたくじで、平均を取ることに意味があります。ところが株価は違います。1年前の株価が従っていた分布と今日の株価が従っている分布は別物で、その100日ぶんの平均は「過去100日にたまたま通過した水準の平均」でしかなく、明日を予測する材料にはなりません。

定常性は、時系列を扱うほぼすべての技術の土台になっています。まず時系列予測(ARモデル・ARMAモデル・状態空間モデル)は、ほぼ例外なく「対象が定常であること」を前提に理論が組まれています。定常でないデータにARモデルをそのまま当てはめると、係数の推定値が意味を失い、予測区間も信用できません。次に信号処理では、パワースペクトル密度という概念そのものが定常性の上に立っています。自己相関関数のフーリエ変換がスペクトルになるというウィーナー・ヒンチンの定理は、自己相関が「時間差だけの関数」でなければ書き下すことすらできません。さらに異常検知では、「平常時の統計量が一定であること」を暗黙の前提にして閾値を引きます。平常時の分散がじわじわ増えているのに固定閾値を使えば、誤警報が時間とともに増えていくだけです。

やっかいなのは、「定常である」を厳密に定義しようとすると、すべての有限次元分布が時間シフトに対して不変、という非常に強い要求になってしまい、実データで確かめることが事実上不可能になる点です。そこで実務では、1次と2次のモーメント(平均と自己共分散)だけに条件を課した「広義定常」という緩い定義を使います。この記事では、その厳密版(強定常)と実務版(広義定常)の関係を整理し、広義定常が課す条件が実データで何を意味するのかを、導出と検証の両面から詰めていきます。

本記事の内容

  • 強定常(狭義定常)と広義定常(弱定常)の定義と、両者の包含関係
  • ガウス過程では強定常と広義定常が一致することの証明
  • 自己共分散関数が満たすべき性質(偶関数性・$|R(\tau)| \le R(0)$・正定値性)の導出
  • AR(1)過程の定常条件 $|\phi| < 1$ とランダムウォークの非定常性
  • トレンド定常と差分定常(単位根)の違い
  • Pythonによる3段構えの検証(AR(1)比較 → 窓分割の目視診断 → ADF検定)

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

定常性とは — 「レシピが変わらない」ということ

定常性を直感的に掴むために、工場のラインを想像してみてください。あるラインでネジを作っています。設計は直径5.00 mm、実際に出てくるネジは加工のばらつきで5.00 mmのまわりに標準偏差0.02 mmで散らばります。このラインが定常であるとは、「今日出てくるネジも、来月出てくるネジも、まったく同じ確率法則(平均5.00、標準偏差0.02)から生まれる」という意味です。刃物が摩耗して直径がじわじわ太くなっていけば、平均が動くので定常ではありません。機械のガタが増えてばらつきだけが広がっていけば、平均は動かなくても分散が動くので、やはり定常ではありません。

ここで重要なのは、定常性はデータそのものの性質ではなく、データを生み出す「レシピ」の性質だという点です。定常な過程から出てくる時系列も、当然ながら上下に揺れます。パッと見て「動いている」からといって非定常とは限りません。逆に、たまたま100サンプルのあいだ平らに見えても、レシピが動いていれば非定常です。私たちが手にしているのは1本の実現値(サンプルパス)だけなので、そこから背後のレシピの不変性を判断しなければならない — ここに定常性の検証の難しさがあります。

もうひとつ、確率過程の見方として「アンサンブル」の視点を持っておくと理解が速くなります。確率過程 $\{X_t\}$ は、同じレシピから生まれる時系列の無限個の束だと思ってください。時刻 $t$ で束をスパッと切ると、そこには無数の値が並んでいて、それが確率変数 $X_t$ の分布です。定常性とは、この「切り口の分布」が $t$ を動かしても変わらない、さらに言えば複数の時刻で同時に切ったときの同時分布も、切る位置を一斉にずらしても変わらない、ということです。

定常な過程と非定常な過程のアンサンブルを比べた模式図。定常では切り口の分布が3時刻とも同じ形、非定常では時刻が後ろになるほど分布が広がる

この図が定常性のイメージそのものです。左は同じレシピから生まれた40本の実現値で、$t=60,150,240$ のどこで縦に切っても、切り口に現れる分布(黒い曲線)はぴたりと同じ形をしています。右は切るたびに分布が横に太っていき、$\pm 2\sigma$ の幅が $\sqrt{t}$ に比例して開いていきます。1本ずつの線は左右どちらも「動いて」見えるのに、定常か非定常かを分けているのは切り口の分布のほうだという点が重要です。

この「切り口の分布が動かない」をどこまで厳密に要求するかで、定常性には強弱2つのバージョンが生まれます。次のセクションで、まず厳密なほう(強定常)から定義します。

強定常(狭義定常)の定義

厳密版の定常性は、「有限個の時刻を選んで同時分布を作り、それを時間方向にまるごと平行移動しても分布が変わらない」と要求します。単に「各時刻の分布が同じ」ではなく、「時刻のを選んだときの同時分布まで同じ」と言っているのがポイントです。

確率過程 $\{X_t\}_{t \in \mathbb{Z}}$ が強定常(strictly stationary、狭義定常)であるとは、任意の正整数 $n$、任意の時刻の組 $t_1 < t_2 < \dots < t_n$、任意のシフト量 $\tau$ に対して、

$$ \begin{equation} F_{X_{t_1}, \dots, X_{t_n}}(x_1, \dots, x_n) = F_{X_{t_1+\tau}, \dots, X_{t_n+\tau}}(x_1, \dots, x_n) \end{equation} $$

が任意の $(x_1, \dots, x_n) \in \mathbb{R}^n$ で成り立つことをいいます。ここで $F$ は同時累積分布関数です。

この定義が何を言っているか、$n$ を小さいほうから追ってみましょう。

$n=1$ とすると、$F_{X_t}(x) = F_{X_{t+\tau}}(x)$ です。つまりすべての時刻で周辺分布がまったく同じ。平均も分散も歪度も尖度も、そもそも分布の形自体が時刻によらない、と要求しています。

$n=2$ とすると、$(X_{t_1}, X_{t_2})$ の同時分布が $(X_{t_1+\tau}, X_{t_2+\tau})$ の同時分布と一致します。$\tau = -t_1$ と選べば、$(X_{t_1}, X_{t_2})$ の同時分布は $(X_0, X_{t_2 – t_1})$ の同時分布と等しい。つまり2点の同時分布は、2点の絶対位置ではなく間隔 $t_2 – t_1$ だけで決まるということです。「1月1日と1月5日のペア」と「7月10日と7月14日のペア」は、間隔が4日で同じなら統計的に区別がつかない、という主張になります。

$n=3, 4, \dots$ と続ければ、3点同時分布も4点同時分布も、すべて相対位置だけで決まることになります。

強定常の定義の図解。n=1では周辺分布がシフト不変、n=2では間隔が同じペアの同時分布が一致することを、3000本のアンサンブルの散布図で確認している

左は $n=1$ の条件で、$t_1$ と $t_1+\tau$ の切り口の分布が同じ形であることを示しています。中央は $n=2$ の条件で、間隔6のペアを2組、時間軸上の離れた場所から取っています。右はその2組の同時分布をアンサンブル3000本で実際に描いたもので、オレンジと緑の雲が完全に重なっています(相関はそれぞれ0.531と0.542、理論値 $\phi^6 = 0.531$)。強定常とは「ペアの絶対位置を変えても、この散布図の形が1ミリも動かない」という要求なのだと押さえてください。

強定常は非常にきれいな性質ですが、実データで検証するとなると絶望的です。$n=1$ の条件、つまり「すべての時刻で周辺分布が同一」だけでも、分布全体を比べる必要があります。$n=2$ 以降になると同時分布の一致を確かめなければならず、有限のサンプルからはほぼ不可能です。しかも、モーメントが存在するとは一言も言っていません。たとえばコーシー分布に従う独立同分布列は、時刻をどうシフトしても同時分布が変わらないので立派な強定常過程ですが、平均すら存在しません。

強すぎて使えない。ならば条件を「実際に推定できる量」だけに絞ればいい — こうして登場するのが広義定常です。

広義定常(弱定常)の定義

私たちが有限のデータから安定して推定できる量は、実質的に平均と共分散、つまり1次と2次のモーメントに限られます。だったら定常性の条件も、そこだけに課してしまおう。これが広義定常の発想です。分布全体の一致という要求を捨てて、「平均が動かない」「共分散が時間差だけで決まる」の2点に絞ります。

確率過程 $\{X_t\}$ が広義定常(wide-sense stationary、WSS)あるいは弱定常(weakly stationary)であるとは、次の3条件を満たすことをいいます。

条件1(2次モーメントの有限性): すべての $t$ で

$$ \begin{equation} E[X_t^2] < \infty \end{equation} $$

条件2(平均一定): ある定数 $\mu$ が存在して、すべての $t$ で

$$ \begin{equation} E[X_t] = \mu \end{equation} $$

条件3(自己共分散が時間差のみに依存): ある関数 $R(\cdot)$ が存在して、すべての $t$ と $\tau$ で

$$ \begin{equation} \mathrm{Cov}(X_t, X_{t+\tau}) = E[(X_t – \mu)(X_{t+\tau} – \mu)] = R(\tau) \end{equation} $$

条件1は、そもそも分散や共分散が定義できるための前提です。条件2と条件3が定常性の本体になります。

条件3の意味を噛み砕いておきましょう。左辺の $\mathrm{Cov}(X_t, X_{t+\tau})$ は本来 $t$ と $\tau$ の2変数関数です。それが $\tau$ だけの1変数関数 $R(\tau)$ に潰れる、と言っています。つまり「4月の値と4月+3日の値の関係」と「11月の値と11月+3日の値の関係」が同じ強さである、ということです。$\tau = 0$ とおけば $R(0) = \mathrm{Var}(X_t)$ となり、分散も $t$ によらず一定という条件が自動的に含まれていることに注意してください。

$R(\tau)$ を自己共分散関数、それを $R(0)$ で割って規格化した

$$ \rho(\tau) = \frac{R(\tau)}{R(0)} $$

自己相関関数と呼びます。$\rho(0) = 1$ です。

広義定常の2つの条件を可視化した図。左は平均が水平かどうか、中央と右は共分散を時刻と時間差の平面にヒートマップで描いたもの

広義定常が過程を覗いている「2つの窓」を絵にしたのがこの図です。窓1(左)は平均の折れ線が水平かどうかだけを見ます。窓2(中央・右)は $\mathrm{Cov}(X_t, X_{t+\tau})$ を横軸 $t$・縦軸 $\tau$ の平面に描いたもので、定常なら横方向($t$ 方向)に動かしても色が変わらず、縦方向($\tau$ 方向)にだけ変化するのが特徴です。右のランダムウォークでは横に進むほど明るくなっており、共分散が $t$ に依存していること、つまり条件3が破れていることが一目でわかります。

広義定常は「1次と2次のモーメントだけ見る」ので、3次以上の構造がどう暴れていても文句を言いません。たとえば分散が時間とともに変わるGARCH型のモデルは、条件付き分散こそ動きますが無条件分散は一定なので、広義定常になり得ます。ここには「見ているものが浅い」という弱点と、「実際に推定できるものだけを条件にしている」という実用上の強みが同居しています。

さて、こうして2つの定常性が出そろいました。名前からして「強」が「弱」を含んでいそうですが、本当にそうでしょうか。次はこの包含関係を丁寧に確かめます。

強定常と広義定常の包含関係

素朴には「強定常なら当然、広義定常だろう」と思えます。実際そのとおりなのですが、ひとつだけ条件が要ります。

主張: $\{X_t\}$ が強定常で、かつ $E[X_t^2] < \infty$ ならば、$\{X_t\}$ は広義定常である。

証明は定義から直接たどれます。強定常の $n=1$ の条件より、$X_t$ と $X_{t+\tau}$ は同一の分布を持ちます。分布が同じなら期待値も同じなので

$$ E[X_t] = E[X_{t+\tau}] \quad \text{(すべての } t, \tau \text{ で)} $$

となり、$E[X_t]$ は $t$ に依存しない定数 $\mu$ です。これで条件2が言えました。

次に $n=2$ の条件を使います。$(X_t, X_{t+\tau})$ の同時分布は、シフト量を $-t$ と選ぶことで $(X_0, X_\tau)$ の同時分布に一致します。共分散は同時分布から計算される量なので

$$ \mathrm{Cov}(X_t, X_{t+\tau}) = \mathrm{Cov}(X_0, X_\tau) =: R(\tau) $$

となり、右辺は $t$ を含みません。これで条件3が言えました。$E[X_t^2] < \infty$ は仮定しているので条件1も満たされ、広義定常が従います。

ここで $E[X_t^2] < \infty$ という但し書きが本質的である点を強調しておきます。2次モーメントが有限でなければ、強定常でも広義定常にはなりません。反例はさきほど触れたコーシー分布です。$X_t \sim \text{Cauchy}(0,1)$ が独立同分布に並んでいるとき、任意の時刻の組の同時分布は積の形になり、シフトしても変わらないので強定常です。しかしコーシー分布は平均も分散も存在しないので、広義定常の条件1・2をそもそも書き下せません。「強いほうが弱いほうを含む」という直感は、モーメントの存在という橋がないと渡れないわけです。

逆は成り立ちません。広義定常でも強定常とは限らない例を作ってみましょう。互いに独立な確率変数の列 $\{X_t\}$ を、次のように定めます。

  • $t$ が偶数のとき: $X_t \sim \mathcal{N}(0, 1)$
  • $t$ が奇数のとき: $X_t$ は $+1$ と $-1$ を確率 $1/2$ ずつ取る

どの $t$ でも $E[X_t] = 0$、$\mathrm{Var}(X_t) = 1$ です($\pm 1$ の分布も分散1)。独立なので $\tau \ne 0$ のとき $\mathrm{Cov}(X_t, X_{t+\tau}) = 0$。したがって

$$ R(\tau) = \begin{cases} 1 & (\tau = 0) \\ 0 & (\tau \ne 0) \end{cases} $$

となり、$t$ に依存しません。文句なしの広義定常です。ところが1次元の周辺分布を見ると、偶数時刻は連続な正規分布、奇数時刻は2点しか取らない離散分布で、まったくの別物です。強定常の $n=1$ の条件からして破れています。

強定常と広義定常の包含関係を示すベン図。広義定常の外側にコーシー分布のiid列、内側に偶数正規・奇数±1の反例が配置されている

2つの定常性の関係を1枚にまとめると上のようになります。青い外側の楕円が広義定常、赤い内側が「強定常かつ2次モーメント有限」で、赤は青に完全に含まれます。しかし強定常そのものは青からはみ出します(反例B:コーシー分布のi.i.d.列)し、青の中には強定常でないもの(反例A)もいます。「強いほうが弱いほうを含む」という直感が、モーメントの存在という条件抜きには成り立たないことを、この図は視覚的に示しています。

偶数時刻が正規分布・奇数時刻が±1の2値をとる反例の数値検証。系列・周辺分布・標本自己相関の3枚のパネル

反例Aを数値で確かめたのがこの図です。左の系列を見ると、青い点(偶数時刻)は連続的な値を取り、赤い点(奇数時刻)は $\pm 1$ にしか現れません。中央の周辺分布は明確に別物ですが、右の標本自己相関はラグ1以上がすべて信頼帯の中に収まり、実測でも平均0・分散1(偶数時刻0.996、奇数時刻1.000)です。広義定常の3条件はすべて満たしているのに、強定常の $n=1$ 条件だけが破れているという状況が、実際に構成できることがわかります。

このように、広義定常は「平均と共分散という2つの窓」からしか過程を見ていないため、その窓に映らない違いは見逃してしまいます。逆に言えば、この2つの窓に映るものがすべてになる状況では、強定常と広義定常は一致するはずです。そんな都合のいい状況が実在します。ガウス過程です。

ガウス過程では両者が一致する

多変量正規分布には、他の分布にはない特別な性質があります。分布が平均ベクトルと共分散行列だけで完全に決まるという性質です。3次モーメント以上に自由度がなく、1次と2次のモーメントを指定した瞬間に分布が一意に定まります。広義定常は1次と2次のモーメントを縛る条件でしたから、ガウス過程においてはそれが分布全体を縛ることになり、強定常に格上げされるはずです。これを確認しましょう。

確率過程 $\{X_t\}$ がガウス過程であるとは、任意の $n$ と任意の時刻の組 $t_1, \dots, t_n$ に対して、ベクトル $\bm{X} = (X_{t_1}, \dots, X_{t_n})^\top$ が多変量正規分布に従うことをいいます。

主張: ガウス過程においては、広義定常 $\iff$ 強定常。

まず $\Leftarrow$ 方向は簡単です。ガウス過程は定義から2次モーメントが有限(正規分布のモーメントはすべて有限)なので、前セクションの主張がそのまま使えて、強定常ならば広義定常です。

問題は $\Rightarrow$ 方向、つまり「広義定常ならば強定常」です。$\{X_t\}$ をガウス過程かつ広義定常とし、平均 $\mu$、自己共分散関数 $R(\tau)$ を持つとします。

任意の $n$、任意の $t_1, \dots, t_n$、任意のシフト $\tau$ を取り、2つのベクトル

$$ \bm{X} = (X_{t_1}, \dots, X_{t_n})^\top, \qquad \bm{X}^{(\tau)} = (X_{t_1 + \tau}, \dots, X_{t_n + \tau})^\top $$

を比べます。どちらもガウス過程の定義から多変量正規分布に従うので、それぞれの分布は平均ベクトルと共分散行列だけで決まります。

$\bm{X}$ の平均ベクトルは、広義定常の条件2(平均一定)より

$$ E[\bm{X}] = (\mu, \mu, \dots, \mu)^\top = \mu \bm{1} $$

です。$\bm{X}^{(\tau)}$ についても各成分の期待値は同じ $\mu$ なので

$$ E[\bm{X}^{(\tau)}] = \mu \bm{1} $$

となり、平均ベクトルは完全に一致します。

次に共分散行列を比べます。$\bm{X}$ の共分散行列を $\bm{\Sigma}$ とすると、その $(i,j)$ 成分は広義定常の条件3より

$$ \Sigma_{ij} = \mathrm{Cov}(X_{t_i}, X_{t_j}) = R(t_j – t_i) $$

です。同じことを $\bm{X}^{(\tau)}$ について計算します。$(i,j)$ 成分は $\mathrm{Cov}(X_{t_i + \tau}, X_{t_j + \tau})$ ですが、条件3は「時間差だけで決まる」と言っているので、この2つの時刻の差を取ると

$$ (t_j + \tau) – (t_i + \tau) = t_j – t_i $$

と $\tau$ が打ち消し合います。したがって

$$ \mathrm{Cov}(X_{t_i + \tau}, X_{t_j + \tau}) = R(t_j – t_i) = \Sigma_{ij} $$

となり、共分散行列も完全に一致します。シフト量 $\tau$ が引き算で消える — これが証明の核心です。

時刻の組を5だけシフトしても共分散行列が完全に一致することを示すヒートマップ2枚と、証明の要点をまとめたパネル

証明の核心を目で見えるようにしたのがこの図です。左は時刻の組 $(0,1,2,4,7)$ の共分散行列、中央は同じ組を $\tau=5$ ずらした $(5,6,7,9,12)$ の共分散行列で、25個の成分がすべて同じ値になっています。時刻の絶対位置は変わっているのに、差 $t_j – t_i$ が変わらないので行列が動かないわけです。多変量正規分布は平均ベクトルと共分散行列だけで決まるので、これで同時分布の一致、すなわち強定常が言えたことになります。

平均ベクトルと共分散行列が一致したので、$\bm{X}$ と $\bm{X}^{(\tau)}$ は同一の多変量正規分布に従います。特性関数で書けば

$$ E[e^{i \bm{u}^\top \bm{X}}] = \exp\left( i \mu \bm{u}^\top \bm{1} – \frac{1}{2} \bm{u}^\top \bm{\Sigma} \bm{u} \right) = E[e^{i \bm{u}^\top \bm{X}^{(\tau)}}] $$

が任意の $\bm{u} \in \mathbb{R}^n$ で成り立ち、特性関数が一致すれば分布も一致します(特性関数の一意性定理)。$n$ と時刻の組と $\tau$ は任意に取っていたので、これは強定常の定義そのものです。証明完了です。

この結果は実務上とても大きな意味を持ちます。ガウス性を仮定できる場面では、「平均一定・自己共分散が時間差のみに依存」を確認するだけで、分布レベルの定常性まで保証されるからです。カルマンフィルタが前提とする線形ガウシアン状態空間モデル、通信系の熱雑音、ガウス過程回帰など、応用の広い領域がまるごとこの恩恵にあずかります。逆に、金融時系列のように裾が厚くガウス性が疑わしい対象では、広義定常が成り立っても分布形が動いている可能性を捨てきれません。この違いは意識しておく価値があります。

ところで、広義定常の条件3では「$R(\tau)$ なる関数が存在する」と書きましたが、$R(\tau)$ はどんな関数でもいいのでしょうか。実は、自己共分散関数を名乗るには満たすべき制約があります。次はそれを導きます。

自己共分散関数が満たすべき3つの性質

$R(\tau)$ は「時間差 $\tau$ 離れた2点がどれくらい連動するか」を表す関数でした。この関数は好き勝手な形を取れるわけではありません。共分散という出自から、3つの制約が自動的に課されます。これを知っていると、推定した自己相関がおかしいときに「そもそも自己共分散関数としてあり得ない形だ」と気づけるようになります。

性質1: 偶関数性 $R(-\tau) = R(\tau)$

まず、負の時間差を考えたときの値です。定義に戻ると

$$ R(-\tau) = \mathrm{Cov}(X_t, X_{t-\tau}) $$

です。共分散は2つの引数の順序を入れ替えても値が変わらない($\mathrm{Cov}(A,B) = \mathrm{Cov}(B,A)$)ので

$$ R(-\tau) = \mathrm{Cov}(X_{t-\tau}, X_t) $$

と書き直せます。ここで $s = t – \tau$ とおくと、右辺は $\mathrm{Cov}(X_s, X_{s+\tau})$ です。広義定常なので、この値は $s$ によらず $R(\tau)$ に等しい。したがって

$$ \begin{equation} R(-\tau) = R(\tau) \end{equation} $$

が示せました。自己共分散関数は原点対称の偶関数です(実数値過程の場合。複素数値過程では $R(-\tau) = \overline{R(\tau)}$ とエルミート対称になります)。この性質のおかげで、自己相関関数のグラフを描くときは $\tau \ge 0$ の側だけ見れば十分です。

性質2: 原点で最大 $|R(\tau)| \le R(0)$

自己共分散の絶対値は、$\tau = 0$ での値、すなわち分散を超えられません。「自分自身との連動より強く連動できる相手はいない」という当たり前の話ですが、きちんと導いておきましょう。

コーシー・シュワルツの不等式を使えば1行です。

$$ |\mathrm{Cov}(X_t, X_{t+\tau})| \le \sqrt{\mathrm{Var}(X_t) \cdot \mathrm{Var}(X_{t+\tau})} = \sqrt{R(0) \cdot R(0)} = R(0) $$

ただ、これだと「なぜそうなるか」が見えにくいので、分散の非負性だけを使う初等的な導出も示します。任意の実数について、和と差の分散を考えます。

$$ \mathrm{Var}(X_{t+\tau} + X_t) \ge 0, \qquad \mathrm{Var}(X_{t+\tau} – X_t) \ge 0 $$

分散の展開公式 $\mathrm{Var}(A \pm B) = \mathrm{Var}(A) + \mathrm{Var}(B) \pm 2\mathrm{Cov}(A,B)$ を適用します。広義定常より $\mathrm{Var}(X_t) = \mathrm{Var}(X_{t+\tau}) = R(0)$、$\mathrm{Cov}(X_t, X_{t+\tau}) = R(\tau)$ を代入すると

$$ 2R(0) + 2R(\tau) \ge 0, \qquad 2R(0) – 2R(\tau) \ge 0 $$

の2本が得られます。それぞれを2で割って整理すると、第1式から $R(\tau) \ge -R(0)$、第2式から $R(\tau) \le R(0)$。両方を合わせて

$$ \begin{equation} -R(0) \le R(\tau) \le R(0) \quad \Longleftrightarrow \quad |R(\tau)| \le R(0) \end{equation} $$

が示せました。自己相関で書けば $|\rho(\tau)| \le 1$ です。相関係数が $[-1, 1]$ に収まるという事実の、時系列版だと思えばわかりやすいはずです。

自己共分散関数の性質1と2を示す図。左は複数のARモデルのR(τ)が原点対称であること、右はR(τ)が±R(0)の帯の中に収まり原点で最大になること

左のパネルでは、$\phi = 0.5, 0.85, -0.6$ という3つのAR(1)の自己共分散関数がいずれも $\tau = 0$ を軸に左右対称になっています。$\phi = -0.6$ のように符号が交互に振れる場合でも対称性は保たれます。右のパネルは性質2で、$R(\tau)$ が水色の帯 $[-R(0), +R(0)]$ から決して出られず、$\tau = 0$(星印)で必ず最大値 $R(0) = 3.60$ を取ることを示しています。「自分自身との連動より強く連動できる相手はいない」が図として現れているわけです。

性質3: 正定値性

3つ目が最も重要で、かつ見落とされがちな性質です。「$\rho(1) = 0.9$、$\rho(2) = 0$」のような自己相関関数は、一見もっともらしく見えますが存在し得ません。その理由が正定値性です。

任意の $n$、任意の時刻 $t_1, \dots, t_n$、任意の実数 $a_1, \dots, a_n$ を取り、線形結合

$$ Y = \sum_{i=1}^{n} a_i X_{t_i} $$

を作ります。$Y$ は確率変数なので、その分散は非負でなければなりません。

$$ \mathrm{Var}(Y) \ge 0 $$

この左辺を展開します。分散の双線形性から

$$ \mathrm{Var}\left( \sum_{i=1}^{n} a_i X_{t_i} \right) = \sum_{i=1}^{n} \sum_{j=1}^{n} a_i a_j \, \mathrm{Cov}(X_{t_i}, X_{t_j}) $$

となります。ここで広義定常の条件3を代入すると、$\mathrm{Cov}(X_{t_i}, X_{t_j}) = R(t_i – t_j)$ なので

$$ \begin{equation} \sum_{i=1}^{n} \sum_{j=1}^{n} a_i a_j \, R(t_i – t_j) \ge 0 \end{equation} $$

が任意の $\{a_i\}$ について成り立ちます。これは行列 $\bm{G}$($G_{ij} = R(t_i – t_j)$、自己共分散行列あるいはトープリッツ行列と呼ばれます)が半正定値であることを意味します。この条件を、$R$ が正定値関数である、と言います。

正定値性がどれだけ強い制約かを、具体例で見ましょう。時刻 $t, t+1, t+2$ の3点を取り、$\rho(0)=1$、$\rho(1)=0.9$、$\rho(2)=0$ という自己相関を仮定します。相関行列は

$$ \bm{G} = \begin{pmatrix} 1 & 0.9 & 0 \\ 0.9 & 1 & 0.9 \\ 0 & 0.9 & 1 \end{pmatrix} $$

です。この固有値を計算すると $-0.273,\ 1.000,\ 2.273$ となり、負の固有値が現れます。負の固有値に対応する固有ベクトルを $\bm{a}$ として線形結合 $Y = \sum a_i X_{t_i}$ を作れば、$\mathrm{Var}(Y) < 0$ という不可能な事態になります。つまりこの自己相関関数を持つ確率過程は存在しません。

なぜ存在しないのか、直感的にも説明できます。$X_t$ と $X_{t+1}$ が相関0.9でほぼ一心同体、$X_{t+1}$ と $X_{t+2}$ も相関0.9でほぼ一心同体なら、$X_t$ と $X_{t+2}$ が無関係でいられるはずがありません。相関係数の三つ組が満たすべき不等式から、$\rho_{13}$ の下限は

$$ \rho_{13} \ge \rho_{12}\rho_{23} – \sqrt{(1-\rho_{12}^2)(1-\rho_{23}^2)} = 0.81 – \sqrt{0.19 \times 0.19} = 0.81 – 0.19 = 0.62 $$

と計算でき、$\rho(2)$ は少なくとも0.62 なければなりません。0 は論外というわけです。

正定値性の検証図。ρ(1)=0.9, ρ(2)=0を仮定した相関行列と、その固有値に-0.273という負の値が現れること、ρ(2)が取れる範囲

左は仮定した相関行列そのもの、中央はその固有値で、$-0.273$ という負の値がはっきり現れています。右のパネルは $\rho(1) = 0.9$ を固定したまま $\rho(2)$ を $-1$ から $1$ まで動かし、相関行列の最小固有値を追跡したものです。最小固有値が0以上になるのは $\rho(2) \ge 0.62$ の緑の領域だけで、$\rho(2) = 0$(赤い×印)は完全に外側にあります。自己相関関数は3つの値を自由に決められるものではない、ということが定量的に見て取れます。

さらに深い事実として、ボホナーの定理(離散時間ではヘルグロッツの定理)があります。$R(\tau)$ が正定値関数であることと、ある非負の有限測度 $F$(スペクトル分布関数)が存在して

$$ R(\tau) = \int_{-\pi}^{\pi} e^{i \omega \tau} \, dF(\omega) $$

と表せることが同値である、という定理です。$F$ が密度 $S(\omega)$ を持つ場合、$S(\omega) \ge 0$ がパワースペクトル密度にほかなりません。つまり「自己共分散関数が正定値である」ことと「パワースペクトル密度が非負である」ことは同じことを言っており、これが自己相関関数とはで扱ったウィーナー・ヒンチンの定理の土台になっています。正定値性を破った自己相関関数は、フーリエ変換したときに負のパワーを持つ周波数が現れてしまい、物理的に意味をなしません。

ここまでで、広義定常の定義と自己共分散関数の制約が揃いました。次は最も身近な定常/非定常の対比、AR(1)過程とランダムウォークを、手で計算して確かめます。

具体例1: AR(1)過程はいつ定常になるか

時系列解析で最初に出会うモデル、1次自己回帰過程を取り上げます。

$$ \begin{equation} X_t = \phi X_{t-1} + \varepsilon_t, \qquad \varepsilon_t \sim \text{WN}(0, \sigma^2) \end{equation} $$

ここで $\text{WN}(0,\sigma^2)$ は平均0・分散 $\sigma^2$ の白色雑音(互いに無相関)です。「前の値を $\phi$ 倍して引き継ぎ、新しい衝撃を足す」という単純な仕組みですが、$\phi$ の値ひとつで定常にも非定常にもなります。

平均を求める

まず定常解が存在すると仮定して、$E[X_t] = \mu$($t$ によらない定数)とおきます。両辺の期待値を取ると

$$ E[X_t] = \phi E[X_{t-1}] + E[\varepsilon_t] $$

$E[\varepsilon_t] = 0$ と $E[X_t] = E[X_{t-1}] = \mu$ を代入すると

$$ \mu = \phi \mu \quad \Longrightarrow \quad \mu(1 – \phi) = 0 $$

$\phi \ne 1$ ならば $\mu = 0$ が唯一の解です。逆に $\phi = 1$ だとこの式は $\mu$ について何の情報も与えず、定数平均が一意に定まりません。すでにここで $\phi = 1$ の特殊性が顔を出しています。

分散を求める

次に $\mathrm{Var}(X_t) = \gamma_0$($t$ によらない)とおきます。$X_{t-1}$ と $\varepsilon_t$ は無相関($\varepsilon_t$ は時刻 $t$ で初めて現れる新しい衝撃なので、過去の $X_{t-1}$ とは独立)であることに注意して、両辺の分散を取ると

$$ \mathrm{Var}(X_t) = \phi^2 \mathrm{Var}(X_{t-1}) + \mathrm{Var}(\varepsilon_t) $$

$\gamma_0 = \phi^2 \gamma_0 + \sigma^2$ を $\gamma_0$ について解くと

$$ \gamma_0 (1 – \phi^2) = \sigma^2 \quad \Longrightarrow \quad \gamma_0 = \frac{\sigma^2}{1 – \phi^2} $$

ここが決定的です。分散が正の有限値であるためには $1 – \phi^2 > 0$、すなわち $|\phi| < 1$ が必要です。$|\phi| \ge 1$ では、この式が負の分散や無限大を返してしまい、定常解が存在しません。

自己共分散を求める

$k \ge 1$ について $\gamma_k = \mathrm{Cov}(X_t, X_{t-k})$ を計算します。定義式 $X_t = \phi X_{t-1} + \varepsilon_t$ の両辺に $X_{t-k}$ を掛けて期待値を取ります($\mu = 0$ なので共分散=積の期待値)。

$$ E[X_t X_{t-k}] = \phi E[X_{t-1} X_{t-k}] + E[\varepsilon_t X_{t-k}] $$

$k \ge 1$ のとき $\varepsilon_t$ は $X_{t-k}$ より未来の衝撃なので $E[\varepsilon_t X_{t-k}] = 0$ です。したがって

$$ \gamma_k = \phi \gamma_{k-1} $$

という漸化式が得られます。$\gamma_0$ から出発して繰り返し適用すると

$$ \begin{equation} \gamma_k = \phi^k \gamma_0 = \frac{\phi^k \sigma^2}{1 – \phi^2}, \qquad \rho_k = \frac{\gamma_k}{\gamma_0} = \phi^k \end{equation} $$

自己相関が $\phi^k$ という指数関数的な減衰を示します。$\phi = 0.5$ なら $\rho_1 = 0.5$, $\rho_2 = 0.25$, $\rho_3 = 0.125$ と半分ずつ減り、10ラグも離れればほぼ0です。「昔のことは忘れる」という性質が、定常性の直感的な中身になっています。

MA($\infty$)表現からも見える

もう一つの見方として、漸化式を過去に向かって展開してみます。

$$ 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 $$

これを無限に続けると

$$ X_t = \sum_{j=0}^{\infty} \phi^j \varepsilon_{t-j} $$

という白色雑音の無限和(MA($\infty$)表現)が得られます。この和が平均二乗収束するには、係数の二乗和が有限、すなわち

$$ \sum_{j=0}^{\infty} \phi^{2j} = \frac{1}{1 – \phi^2} < \infty $$

が必要で、これも $|\phi| < 1$ に帰着します。この表現は「過去の衝撃 $\varepsilon_{t-j}$ の影響が $\phi^j$ という重みで減衰していく」ことを直接示していて、定常性の物理的意味を最もよく表しています。衝撃が忘れられるから定常なのです。

AR(1)の定常条件を3つの角度から示す図。自己相関の指数減衰、φに対する定常分散の発散、MA(∞)係数の減衰

3枚のパネルが同じ結論 $|\phi| < 1$ を別々の側面から示しています。左は $\rho(k) = \phi^k$ で、$\phi = 0.99$ ではラグ20でもまだ0.82と、ほとんど忘れていないことが読み取れます。中央は定常分散 $\sigma^2/(1-\phi^2)$ を対数目盛で描いたもので、$|\phi| \to 1$ で垂直に発散し、赤い破線の外側には定常解が存在しません。右はMA($\infty$)の重みで、$\phi = 0.5$ では10ステップ前の衝撃の影響が0.001以下に潰れるのに対し、$\phi = 1.0$ では何ステップ前だろうと重み1のまま残ります。この右のパネルが、定常か非定常かを分けている構造そのものです。

では $\phi = 1$ のとき、衝撃はどうなるのでしょうか。次のセクションで確かめます。

具体例2: ランダムウォークはなぜ非定常か

$\phi = 1$ とすると、AR(1)は

$$ X_t = X_{t-1} + \varepsilon_t $$

というランダムウォーク(酔歩)になります。$X_0 = 0$ から出発すると、繰り返し代入して

$$ X_t = \varepsilon_1 + \varepsilon_2 + \dots + \varepsilon_t = \sum_{s=1}^{t} \varepsilon_s $$

と、衝撃の累積和になります。MA($\infty$)表現で言えば係数がすべて1、つまり過去の衝撃が一切減衰せず、そのまま残り続ける構造です。

平均を計算すると

$$ E[X_t] = \sum_{s=1}^{t} E[\varepsilon_s] = 0 $$

で、なんと定数です。広義定常の条件2は満たしてしまいます。ここが引っかけどころで、「ランダムウォークは平均が動くから非定常」というよくある説明は正確ではありません。問題は分散にあります。$\varepsilon_s$ は互いに無相関なので、分散は単純に足し算になり

$$ \begin{equation} \mathrm{Var}(X_t) = \sum_{s=1}^{t} \mathrm{Var}(\varepsilon_s) = t \sigma^2 \end{equation} $$

分散が時間に比例して無限に増大します。$R(0) = \mathrm{Var}(X_t)$ が $t$ に依存する時点で、広義定常の条件3が破れています。時間が経つほど、その時点での値は原点から遠くへ散らばっていく。これがランダムウォークの非定常性の正体です。

自己共分散も計算しておきましょう。$k > 0$ として

$$ \mathrm{Cov}(X_t, X_{t+k}) = \mathrm{Cov}\left( \sum_{s=1}^{t} \varepsilon_s, \sum_{u=1}^{t+k} \varepsilon_u \right) $$

共通して現れるのは $s = u \le t$ の項だけで、それ以外は無相関なので消えます。したがって

$$ \mathrm{Cov}(X_t, X_{t+k}) = \sum_{s=1}^{t} \mathrm{Var}(\varepsilon_s) = t\sigma^2 $$

これも $t$ に依存し、時間差 $k$ だけの関数になっていません。相関係数に直すと

$$ \mathrm{Corr}(X_t, X_{t+k}) = \frac{t\sigma^2}{\sqrt{t\sigma^2 \cdot (t+k)\sigma^2}} = \sqrt{\frac{t}{t+k}} $$

という興味深い形になります。$t$ が大きく $k$ が小さければ、この値は1に極めて近い。たとえば $t = 500$, $k = 5$ なら $\sqrt{500/505} = 0.995$ です。十分に時間が経ったランダムウォークでは、隣接する値どうしの相関がほぼ1になるわけです。標本自己相関を計算すると、ラグを増やしてもなかなか0に落ちない、ゆるやかに右下がりの直線状のプロットになります。これは実務で非定常性を疑う代表的なサインです。

ランダムウォークは1階差分を取ると

$$ \Delta X_t = X_t – X_{t-1} = \varepsilon_t $$

と白色雑音になり、これは立派な広義定常です。「差分を取れば定常になる」過程を差分定常(difference stationary)あるいは単位根過程(unit root process)と呼びます。「単位根」という名前は、AR多項式 $1 – \phi z = 0$ の根が $z = 1/\phi = 1$(単位円上)になることに由来します。

ランダムウォークが非定常な理由を示す3枚のパネル。衝撃の累積和、相関が√(t/(t+k))になること、1階差分すると自己相関が消えること

左のパネルでは、灰色の棒が各時刻の衝撃 $\varepsilon_t$、赤い折れ線がその累積和です。一度入った衝撃が差し引かれずに水準に残り続けるため、系列はどこまでも漂っていきます。中央は $\mathrm{Corr}(X_t, X_{t+k}) = \sqrt{t/(t+k)}$ を起点 $t$ ごとに描いたもので、同じ時間差 $k=5$ でも $t=10$ と $t=500$ で値がまるで違う(0.816 と 0.995)ことがわかります。これが「時間差だけの関数になっていない」ということの意味です。右は1階差分の効果で、差分前(赤)はラグ20でも0.6近い相関が残るのに対し、差分後(青)はラグ1以降が信頼帯に収まり、白色雑音=広義定常になっています。

ところで、時間とともに水準が上がっていく時系列には、もうひとつ別のタイプがあります。両者は見た目が似ているのに、扱い方が正反対です。次で区別します。

トレンド定常と差分定常の決定的な違い

右肩上がりのグラフを見たとき、その背後には2つの異なるメカニズムがあり得ます。

タイプA: トレンド定常(trend stationary)

$$ X_t = a + bt + Y_t, \qquad \{Y_t\} \text{ は広義定常} $$

決定論的な直線に、定常な揺らぎが乗っているだけの構造です。$E[X_t] = a + bt$ が $t$ に依存するので広義定常ではありませんが、回帰で直線を引いて残差を取れば定常になります

タイプB: 差分定常(単位根)

$$ X_t = \delta + X_{t-1} + \varepsilon_t $$

ドリフト付きランダムウォークです。累積すると $X_t = X_0 + \delta t + \sum_{s=1}^{t}\varepsilon_s$ となり、こちらも平均は $X_0 + \delta t$ と直線的に増えます。見た目はタイプAとそっくりです。しかし分散は $t\sigma^2$ と増大し続けます。この場合は差分を取らないと定常になりません

両者を分けるのはショックの持続性です。ある時点で $\varepsilon_t$ に大きな正の衝撃が入ったとします。タイプAでは、$Y_t$ が定常なので衝撃の影響はやがて減衰し、時系列は決定論的トレンド $a + bt$ の上に戻ってきます。長期予測は「トレンド線の延長」でよいわけです。ところがタイプBでは、衝撃が累積和にそのまま加わるので、水準が永久にずれたままになります。長期予測をするなら、現在の水準から直線を引き直さなければなりません。

処方箋も正反対です。トレンド定常に対して差分を取ると、$\Delta X_t = b + \Delta Y_t$ となり、$\Delta Y_t$ は「過剰差分」された系列になります。これは自己相関構造をゆがめ、$\rho(1)$ に人工的な負の相関を持ち込みます。逆に、単位根過程に対してトレンド回帰で除去しようとすると、確率的トレンドは直線では取り除けないので、残差にまだ単位根が残ってしまいます。どちらのタイプなのかを見極めてから前処理を選ぶ必要がある、というのが実務上の要点です。

トレンド定常とドリフト付き単位根の比較図。左と中央は似た右肩上がりの系列、右は衝撃が減衰して戻るか水準が永久にずれるかの違い

左(タイプA)と中央(タイプB)を見比べてください。どちらも右肩上がりで、パッと見では区別がつきません。しかしタイプAは黒い破線の決定論的トレンドに終始張りついているのに対し、タイプBは破線から大きく外れて漂い、その外れ方が時間とともに大きくなっています。右のパネルは $t=30$ に大きさ $+5$ の衝撃を1回だけ入れた思考実験で、トレンド定常(青)は数ステップで元の直線に戻るのに対し、単位根(赤)は水準がずれたまま永久に戻らないという決定的な差が現れています。

そして、この見極めをデータから統計的に行う道具が単位根検定です。次のセクションからPythonで実際に手を動かしていきます。

Pythonでの検証: 3段構えの診断手順

定常性の判定は、いきなり検定に飛びつくのではなく、次の3段階で進めるのが確実です。

  1. 系列を眺める — レベルが漂っていないか、ばらつきが広がっていないか、目で見る
  2. 窓に区切って統計量を比べる — 平均・分散・自己相関を区間ごとに計算し、動いていないか確認する
  3. 統計的検定にかける — ADF検定などで単位根の有無を判定する

この順序が大事です。検定は「単位根があるかどうか」という1つの問いにしか答えてくれません。分散だけが増えている、季節性で平均が周期的に動いている、といった非定常性は検定では拾えないので、1と2の目視診断が欠かせません。

まず共通の準備として、AR(1)を $\phi = 0.5$(定常)と $\phi = 1.0$(単位根)の2通りで生成します。

import numpy as np
import matplotlib
import matplotlib.pyplot as plt

# 日本語フォント設定
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
    if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
        plt.rcParams["font.family"] = cand
        break
plt.rcParams["axes.unicode_minus"] = False


def gen_ar1(phi, n, sigma=1.0, seed=0, x0=0.0):
    """AR(1): X_t = phi * X_{t-1} + e_t を生成する"""
    rng = np.random.default_rng(seed)
    e = rng.normal(0.0, sigma, n)
    x = np.empty(n)
    x[0] = x0 + e[0]
    for t in range(1, n):
        x[t] = phi * x[t - 1] + e[t]
    return x


N = 600
x_st = gen_ar1(0.5, N, seed=42)   # 定常: |phi| < 1
x_rw = gen_ar1(1.0, N, seed=42)   # 非定常: 単位根

print("phi=0.5 : 全体平均=%.3f  全体分散=%.3f  (理論分散=%.3f)"
      % (x_st.mean(), x_st.var(), 1 / (1 - 0.5 ** 2)))
print("phi=1.0 : 全体平均=%.3f  全体分散=%.3f" % (x_rw.mean(), x_rw.var()))

実行すると次のように出力されます。

phi=0.5 : 全体平均=-0.050  全体分散=1.329  (理論分散=1.333)
phi=1.0 : 全体平均=-5.784  全体分散=20.691

$\phi = 0.5$ の系列は、標本分散1.329が理論値 $\sigma^2/(1-\phi^2) = 1.333$ とほぼ一致しました。導出した定常分散の式が正しいことの確認になっています。一方 $\phi = 1.0$ の系列は標本分散が20.691と一桁大きく、しかも「全体分散」という量に理論的な意味がありません。ランダムウォークの分散は各時刻で $t\sigma^2$ と違う値を持つので、全期間をひとまとめにした標本分散は「たまたまこの1本がどれくらい遠くまで散歩したか」を測っているにすぎないからです。

ステップ1: 系列を眺める

2つの系列を並べて描いてみます。

fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharex=True)

axes[0].plot(x_st, color="#1f77b4", lw=1.0)
axes[0].axhline(0, color="black", ls="--", lw=1.0, label="理論平均 $\\mu=0$")
axes[0].fill_between(range(N), -2 * np.sqrt(4 / 3), 2 * np.sqrt(4 / 3),
                     color="#1f77b4", alpha=0.12, label="理論の $\\pm 2\\sigma$ 帯(一定幅)")
axes[0].set_ylabel("$X_t$")
axes[0].set_title("定常: AR(1) $\\phi=0.5$ — 一定の幅の中を行き来する")
axes[0].legend(loc="upper right")
axes[0].grid(alpha=0.3)

t = np.arange(N)
axes[1].plot(x_rw, color="#d62728", lw=1.0)
axes[1].axhline(0, color="black", ls="--", lw=1.0, label="出発点")
axes[1].fill_between(t, -2 * np.sqrt(t), 2 * np.sqrt(t),
                     color="#d62728", alpha=0.12, label="理論の $\\pm 2\\sqrt{t}\\sigma$ 帯(広がる)")
axes[1].set_xlabel("時刻 $t$")
axes[1].set_ylabel("$X_t$")
axes[1].set_title("非定常: AR(1) $\\phi=1.0$(ランダムウォーク)— 帯が $\\sqrt{t}$ で開く")
axes[1].legend(loc="upper left")
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

定常なAR(1)と単位根過程のサンプルパスの比較。上段は一定幅の±2σ帯の中を振動し、下段は帯が√tで開いていく

上段の $\phi = 0.5$ は、原点を中心とした一定幅の帯(薄い青)の中を細かく振動しています。どこかへ行こうとしても、$\phi = 0.5$ の引き戻す力(平均回帰)によって毎ステップ半分に縮められるので、遠くへは行けません。下段の $\phi = 1.0$ は、引き戻す力が消えているため、帯そのものが $\sqrt{t}$ に比例して開いていき、系列はその中を漂います。「一定幅の帯」と「開いていく帯」の対比が、定常と非定常の最も視覚的な違いです。

ステップ2: アンサンブルで確かめる

1本のパスだけ見ていると、「たまたま漂っているだけでは?」という疑いが残ります。そこで同じレシピから多数の系列を生成し、各時刻でのアンサンブル平均・アンサンブル分散を測ってみます。定常性の定義はもともとアンサンブルに関する条件なので、これが最も定義に忠実な確認方法です。

import numpy as np
import matplotlib.pyplot as plt

M, N = 4000, 600   # M本の実現値を同時に生成


def gen_ar1_ensemble(phi, M, N, seed=0):
    rng = np.random.default_rng(seed)
    e = rng.normal(0.0, 1.0, (M, N))
    x = np.zeros((M, N))
    for t in range(1, N):
        x[:, t] = phi * x[:, t - 1] + e[:, t]
    return x


ens_st = gen_ar1_ensemble(0.5, M, N, seed=0)
ens_rw = gen_ar1_ensemble(1.0, M, N, seed=1)

for name, ens in [("phi=0.5", ens_st), ("phi=1.0", ens_rw)]:
    for tt in [50, 100, 300, 599]:
        print("%s t=%3d : アンサンブル平均=%+.3f  アンサンブル分散=%.2f"
              % (name, tt, ens[:, tt].mean(), ens[:, tt].var()))
    print()

出力は次のようになります。

phi=0.5 t= 50 : アンサンブル平均=-0.013  アンサンブル分散=1.35
phi=0.5 t=100 : アンサンブル平均=+0.037  アンサンブル分散=1.32
phi=0.5 t=300 : アンサンブル平均=-0.022  アンサンブル分散=1.29
phi=0.5 t=599 : アンサンブル平均=+0.032  アンサンブル分散=1.31

phi=1.0 t= 50 : アンサンブル平均=+0.091  アンサンブル分散=49.59
phi=1.0 t=100 : アンサンブル平均=+0.410  アンサンブル分散=100.94
phi=1.0 t=300 : アンサンブル平均=+0.545  アンサンブル分散=304.32
phi=1.0 t=599 : アンサンブル平均=+0.535  アンサンブル分散=632.01

4000本のアンサンブルから測った平均と分散の時間変化。左は平均がどちらも0付近、右は分散が定常なら一定・単位根なら時刻に比例して増大

この結果を時刻方向に連続的に描いたのが上の図です。左のパネルでは、単位根系列(赤)のアンサンブル平均も $\pm 0.9$ 以内をふらついているだけで、系統的な上昇や下降は見られません。右のパネルでは、定常系列(青)の分散が理論値 $4/3$ の水平線に張りついているのに対し、単位根系列(赤)の分散が理論直線 $t\sigma^2$(赤い点線)にぴたりと重なって直線的に増大しています。非定常性は左(平均)ではなく右(分散)に現れていることが一目でわかります。

この数値は理論を見事に裏づけています。$\phi = 0.5$ ではアンサンブル分散が $t$ によらず 1.3 前後で一定(理論値 $\sigma^2/(1-\phi^2) = 1.333$)です。対して $\phi = 1.0$ では、$t = 50, 100, 300, 599$ に対して分散が 49.6, 100.9, 304.3, 632.0 と、ほぼ $t\sigma^2 = t$ そのものになっています。導出した $\mathrm{Var}(X_t) = t\sigma^2$ が数値実験で確認できました(4000本での推定なので、$t=599$ での相対誤差 5% 程度は推定のばらつきの範囲です)。

さらに注目したいのがアンサンブル平均です。$\phi = 1.0$ でも平均は 0 の近傍($+0.09$ 〜 $+0.55$)に留まっており、標準偏差が $\sqrt{t} \approx 24$ にもなる分布の中心としては、$t$ とともに系統的に動いているとは言えません。実際、アンサンブル平均の推定誤差は $\sqrt{t/M} = \sqrt{599/4000} \approx 0.39$ なので、$+0.535$ はその範囲内です。前のセクションで指摘したとおり、ランダムウォークが非定常なのは平均が動くからではなく、分散が動くからだという事実が、数値からもはっきり読み取れます。

ステップ3: 窓に区切って目視診断する

実データでは実現値は1本しかないので、アンサンブルは作れません。そこで時間方向に窓を切り、区間ごとの統計量を比べるという現実的な代替手段を使います。定常でエルゴード的な過程なら、どの窓で計算しても統計量はほぼ同じ値になるはずです。

import numpy as np
import matplotlib.pyplot as plt


def lag1_acf(w):
    """窓 w の lag-1 標本自己相関"""
    w = w - w.mean()
    return np.sum(w[1:] * w[:-1]) / np.sum(w * w)


W = 100                      # 窓幅
nw = N // W                  # 窓の数
centers = [(k + 0.5) * W for k in range(nw)]

stats = {}
for name, x in [("$\\phi=0.5$(定常)", x_st), ("$\\phi=1.0$(単位根)", x_rw)]:
    ms, vs, a1 = [], [], []
    for k in range(nw):
        w = x[k * W:(k + 1) * W]
        ms.append(w.mean())
        vs.append(w.var())
        a1.append(lag1_acf(w))
    stats[name] = (ms, vs, a1)
    print(name)
    print("  窓ごとの平均 :", np.round(ms, 2), " 最大−最小=%.2f" % (max(ms) - min(ms)))
    print("  窓ごとの分散 :", np.round(vs, 2), " 最大/最小=%.2f" % (max(vs) / min(vs)))
    print("  lag-1 自己相関:", np.round(a1, 3))

出力は次のとおりです。

$\phi=0.5$(定常)
  窓ごとの平均 : [-0.09 -0.03 -0.12  0.2  -0.06 -0.19]  最大−最小=0.40
  窓ごとの分散 : [0.85 1.33 1.35 1.37 1.64 1.35]  最大/最小=1.94
  lag-1 自己相関: [0.539 0.528 0.482 0.515 0.643 0.463]
$\phi=1.0$(単位根)
  窓ごとの平均 : [ 0.99 -7.74 -9.64 -3.76 -5.81 -8.75]  最大−最小=10.63
  窓ごとの分散 : [ 4.9   4.23  6.36 12.86  8.04 10.25]  最大/最小=3.04
  lag-1 自己相関: [0.901 0.878 0.901 0.942 0.931 0.918]

3つの統計量それぞれに、はっきりした差が出ました。

平均が最も雄弁です。$\phi = 0.5$ では6つの窓平均がすべて $-0.19$ から $0.20$ の狭い範囲に収まり、レンジは 0.40 です。標準誤差の目安 $\sqrt{\gamma_0/W}\approx\sqrt{1.33/100}=0.115$ から見て、まったく自然なばらつきです。ところが $\phi = 1.0$ ではレンジが 10.63 と、系列の標準偏差($\sqrt{20.7}=4.5$)より大きく振れています。窓平均が「ばらつき」ではなく「移動」として見えたら、単位根を疑うべきサインです。

分散も差が出ますが、平均ほど決定的ではありません。$\phi = 0.5$ でも最大/最小比が 1.94 と2倍近くあります。有限の窓幅で分散を推定すればこれくらいは普通に揺れるので、この程度で非定常と断じてはいけません。$\phi = 1.0$ の 3.04 も、単体では判断が微妙です。

lag-1 自己相関は最も安定した指標でした。$\phi = 0.5$ では 0.46〜0.64 で理論値 0.5 のまわりに集中し、$\phi = 1.0$ では 0.88〜0.94 とすべて1に近い値です。前に導出した $\mathrm{Corr}(X_t, X_{t+1}) = \sqrt{t/(t+1)}$ がほぼ1になるという理論と整合しています。窓を変えても自己相関がずっと0.9以上に張りついていたら、単位根の可能性が高いと考えてよいでしょう。

窓診断の結果を1枚のパネルにまとめて描いておくと、報告資料としても使えます。

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
titles = ["窓ごとの平均", "窓ごとの分散", "窓ごとの lag-1 自己相関"]
colors = {"$\\phi=0.5$(定常)": "#1f77b4", "$\\phi=1.0$(単位根)": "#d62728"}

for j in range(3):
    for name, (ms, vs, a1) in stats.items():
        y = [ms, vs, a1][j]
        axes[j].plot(centers, y, "o-", color=colors[name], label=name)
    axes[j].set_title(titles[j])
    axes[j].set_xlabel("窓の中心時刻")
    axes[j].grid(alpha=0.3)

axes[0].axhline(0.0, color="gray", ls=":", label="理論平均 0")
axes[1].axhline(4 / 3, color="gray", ls=":", label="理論分散 4/3")
axes[2].axhline(0.5, color="gray", ls=":", label="理論 $\\rho(1)=0.5$")
for j in range(3):
    axes[j].legend(fontsize=8)

plt.tight_layout()
plt.show()

窓ごとの平均・分散・lag-1自己相関を定常系列と単位根系列で比較した3パネル図。平均パネルでレンジ0.40対10.63の差が出ている

3つのパネルを並べると、定常系列(青)はどのパネルでも灰色の理論線のまわりで水平にとどまり、単位根系列(赤)は平均パネルで大きく上下し、自己相関パネルでは0.9付近の高い位置に張りつくという、一目でわかるパターンになります。定常性の診断は「値が理論線の周りで水平か」を見る作業だと思っておくと、実データでも迷いません。

ステップ4: 標本自己相関の減衰を比べる

ラグ方向に自己相関を並べたコレログラムは、単位根を見つける古典的な道具です。

import numpy as np
import matplotlib.pyplot as plt


def sample_acf(x, K):
    """ラグ 0..K の標本自己相関"""
    x = x - x.mean()
    n = len(x)
    c0 = np.sum(x * x) / n
    return np.array([np.sum(x[k:] * x[:n - k]) / n / c0 for k in range(K + 1)])


K = 30
acf_st, acf_rw = sample_acf(x_st, K), sample_acf(x_rw, K)
lags = np.arange(K + 1)

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
axes[0].bar(lags, acf_st, color="#1f77b4", width=0.6)
axes[0].plot(lags, 0.5 ** lags, "k--", lw=1.2, label="理論 $\\rho(k)=0.5^k$")
axes[0].set_title("定常 $\\phi=0.5$: 指数的に速く減衰")

axes[1].bar(lags, acf_rw, color="#d62728", width=0.6)
axes[1].set_title("単位根 $\\phi=1.0$: ほぼ直線でゆっくり減衰")

for ax in axes:
    ax.axhline(0, color="black", lw=0.8)
    ax.axhline(1.96 / np.sqrt(N), color="gray", ls=":", label="$\\pm 1.96/\\sqrt{N}$")
    ax.axhline(-1.96 / np.sqrt(N), color="gray", ls=":")
    ax.set_xlabel("ラグ $k$")
    ax.grid(alpha=0.3)
    ax.legend(fontsize=9)
axes[0].set_ylabel("標本自己相関 $\\hat{\\rho}(k)$")
plt.tight_layout()
plt.show()

print("phi=0.5 の rho(k), k=0..5:", np.round(acf_st[:6], 3))
print("     理論 0.5^k         :", np.round(0.5 ** np.arange(6), 3))
print("phi=1.0 の rho(k), k=0..5:", np.round(acf_rw[:6], 3))

出力は次のとおりです。

phi=0.5 の rho(k), k=0..5: [ 1.     0.537  0.224  0.06  -0.024 -0.033]
     理論 0.5^k         : [1.    0.5   0.25  0.125 0.062 0.031]
phi=1.0 の rho(k), k=0..5: [1.    0.972 0.941 0.912 0.882 0.857]

定常AR(1)と単位根過程のコレログラム比較。左は理論0.5^kに沿って急減衰、右はほぼ直線でゆっくり減衰しラグ30でも0.488

コレログラムにすると差は一層はっきりします。左は棒が黒い破線(理論 $0.5^k$)に沿ってラグ4までに信頼帯へ落ち込み、右はラグ30まで棒が高いまま並んで、破線の直線当てはめ(傾き $-0.0163$/ラグ)にきれいに乗っています。「指数的に落ちて信頼帯に入る」か「直線的に落ちて信頼帯の外に居座り続ける」かが、定常と単位根を分ける視覚的な署名です。

コレログラムの形の違いが決定的です。$\phi = 0.5$ の標本自己相関は 0.537, 0.224, 0.060 と急速に落ち、ラグ4以降は $\pm 1.96/\sqrt{N} = \pm 0.080$ の信頼帯の中に入って統計的にゼロと区別できなくなります。理論値 $0.5^k$ とも良く一致していて、定常AR(1)の $\rho(k) = \phi^k$ が確認できました。

対して $\phi = 1.0$ は 0.972, 0.941, 0.912, 0.882, 0.857 と、ほぼ一定の傾きで直線的に下がっていきます。1ラグあたりの下げ幅は0.016程度で、ラグ30まで行っても0.488と、まだ信頼帯のはるか外側です。この「ゆっくりした線形減衰」こそ単位根の典型的な署名です。理論的には $\mathrm{Corr}(X_t, X_{t+k}) = \sqrt{t/(t+k)}$ が、有限標本での平均化を通じて $1 – k/N$ に近い直線的減衰として現れるためです。実データのコレログラムがこの形をしていたら、まず差分を疑ってください。

目視での診断はここまでです。しかし「ゆっくり減衰している気がする」は主観であり、$\phi = 0.98$ のような境界近くのケースでは判断がつきません。ここで統計的検定の出番になります。

ADF検定で単位根を判定する

ADF検定(Augmented Dickey-Fuller test、拡張ディッキー・フラー検定)は、「単位根があるか」を統計的に判定する標準的な道具です。仕組みを理解しておきましょう。

出発点はAR(1)モデル $X_t = \phi X_{t-1} + \varepsilon_t$ です。単位根の有無は $\phi = 1$ かどうかで決まるので、そのままでは「係数がある特定の値と等しいか」の検定になります。ここで両辺から $X_{t-1}$ を引くという工夫をします。

$$ X_t – X_{t-1} = (\phi – 1) X_{t-1} + \varepsilon_t $$

$\gamma = \phi – 1$ とおけば

$$ \begin{equation} \Delta X_t = \gamma X_{t-1} + \varepsilon_t \end{equation} $$

となり、「$\phi = 1$ か」という問いが「$\gamma = 0$ か」という、ふつうの係数ゼロ検定の形に変換されました。これがディッキー・フラー検定の基本アイデアです。

ADF検定の正体を示す散布図。横軸に前期の値、縦軸に差分を取り、定常なら回帰直線が右下がり、単位根ならほぼ水平になる

この変換を散布図にすると、ADF検定が何を測っているのかが直感的に見えてきます。横軸に前期の値 $X_{t-1}$、縦軸にその後の差分 $\Delta X_t$ を取ると、定常な $\phi=0.5$ では傾き $\gamma = -0.462$(理論値 $\phi – 1 = -0.5$)の右下がりの関係が現れます。「値が高いところにいると次は下げられる」という平均回帰そのものです。一方 $\phi = 1.0$ では傾きが $-0.021$ とほぼ水平で、今どこにいるかが次の動きに影響しない、つまり引き戻す力が存在しないことを示しています。

実データはAR(1)とは限らず、より高次の自己相関を持つのがふつうです。そこで誤差項の系列相関を吸収するために、差分のラグ項を右辺に追加します。定数項とトレンド項も含めた一般形が

$$ \begin{equation} \Delta X_t = \alpha + \beta t + \gamma X_{t-1} + \sum_{i=1}^{p} \delta_i \Delta X_{t-i} + \varepsilon_t \end{equation} $$

です。この「ラグ項を足して拡張した(augmented)」ところがADF検定の名の由来です。

仮説は次のように設定されます。

  • 帰無仮説 $H_0$: $\gamma = 0$(単位根がある = 非定常
  • 対立仮説 $H_1$: $\gamma < 0$($|\phi| < 1$ = 定常

ここで帰無仮説が「非定常」である点を強く意識してください。p値が小さければ「非定常である」を棄却して定常と判断でき、p値が大きければ「非定常を棄却できなかった」だけです。「定常でないことが証明された」わけではありません。この非対称性が、ADF検定の解釈で最も間違えられるポイントです。

もうひとつ技術的な注意があります。$\gamma$ の t 統計量は、通常の t 分布に従いません。$H_0$ のもとで $X_{t-1}$ 自体が非定常なので、標準的な漸近理論が壊れ、ディッキー・フラー分布という別の分布に従います。この分布は左に大きく寄っていて、5%臨界値は $-2.86$ 程度と、通常のt分布の $-1.96$ よりずっと厳しい値になります。統計ソフトはこの専用の臨界値表を内部で使ってくれます。

実際に適用してみましょう。

import numpy as np
from statsmodels.tsa.stattools import adfuller

for name, x in [("phi=0.5(定常)", x_st), ("phi=1.0(単位根)", x_rw)]:
    res = adfuller(x, regression="c", autolag="AIC")
    print("%s : ADF統計量=%7.3f  p値=%.4f  採用ラグ=%d  5%%臨界値=%.3f"
          % (name, res[0], res[1], res[2], res[4]["5%"]))

# 単位根系列を1階差分すると定常になるか
d = np.diff(x_rw)
res = adfuller(d, regression="c", autolag="AIC")
print("phi=1.0 の1階差分 : ADF統計量=%.3f  p値=%.5f" % (res[0], res[1]))

出力は次のとおりです。

phi=0.5(定常) : ADF統計量=-12.833  p値=0.0000  採用ラグ=1  5%臨界値=-2.866
phi=1.0(単位根) : ADF統計量= -2.561  p値=0.1013  採用ラグ=1  5%臨界値=-2.866
phi=1.0 の1階差分 : ADF統計量=-22.355  p値=0.00000

理論どおりの結果です。$\phi = 0.5$ の系列ではADF統計量が $-12.833$ と臨界値 $-2.866$ を大きく下回り、p値は事実上ゼロ。単位根の帰無仮説を明確に棄却して「定常」と判断できます。$\phi = 1.0$ の系列では統計量が $-2.561$ で臨界値に届かず、p値0.1013は5%水準を超えているので帰無仮説を棄却できません。そして単位根系列を1階差分すると統計量は $-22.355$ に跳ね、完全に定常と判定されます。差分定常であることが数値で確認できました。

ここで、$\phi = 1.0$ の結果の言い回しに注意してください。正しくは「単位根がないとは言えなかった」であって、「単位根があると証明された」ではありません。実際、ADF検定の検出力(本当に定常なときに正しく定常と判定できる確率)は、$\phi$ が1に近づくにつれて急速に落ちます。定常なAR(1)を色々な $\phi$ で生成し、棄却率を測ってみましょう。

import numpy as np
import warnings
from statsmodels.tsa.stattools import adfuller
warnings.filterwarnings("ignore")

rng = np.random.default_rng(5)
for phi in [0.5, 0.9, 0.95, 1.0]:
    rej, R = 0, 200
    for _ in range(R):
        e = rng.normal(0, 1, 300)
        x = np.zeros(300)
        for t in range(1, 300):
            x[t] = phi * x[t - 1] + e[t]
        if adfuller(x[50:], regression="c", autolag="AIC")[1] < 0.05:
            rej += 1
    print("phi=%.2f : 「定常」と判定された割合 = %.1f%%" % (phi, 100 * rej / R))

出力は次のとおりです。

phi=0.50 : 「定常」と判定された割合 = 100.0%
phi=0.90 : 「定常」と判定された割合 = 95.0%
phi=0.95 : 「定常」と判定された割合 = 44.5%
phi=1.00 : 「定常」と判定された割合 = 6.5%

ADF検定の検出力を棒グラフで示した図。φ=0.5で100%、0.9で95%、0.95で44.5%、1.0で6.5%と、1に近づくほど検出力が崩れる

棒グラフにすると検出力の崖がよく見えます。$\phi = 0.9$ までは検出力80%の目安(緑の点線)を上回っているのに、$\phi$ を0.05だけ上げて0.95にしただけで44.5%まで落ち込み、目安の半分近くになります。$\phi = 1.0$ の6.5%だけは「失敗」ではなく、有意水準5%どおりの第一種の過誤(灰色の破線)である点も読み分けてください。ADF検定を使う人が必ず知っておくべき事実を示しています。$\phi = 0.5$ なら100%、$\phi = 0.9$ でも95%と高い検出力があるのに、$\phi = 0.95$ になると 44.5%、つまり本当は定常なのに半分以上のケースで「単位根を棄却できない」という結果になってしまいます。$\phi = 1.0$ での棄却率6.5%は、有意水準5%に対する第一種の過誤の実績値で、ほぼ設計どおりです。

$\phi = 0.95$ の系列は、$\rho(k) = 0.95^k$ なのでラグ50でも0.077の相関が残る、非常に長い記憶を持つ定常過程です。250サンプル程度では、単位根との区別がつかないのも無理はありません。「定常性の検定結果は、サンプル数と $\phi$ の組み合わせに強く依存する」ことを踏まえ、検定だけに判断を委ねず、目視診断や理論的な背景知識と組み合わせるのが賢明です。

トレンド定常と単位根をADF検定で見分ける

ADF検定には regression という重要な引数があります。回帰式にどの決定論的項を入れるかの指定で、"c" は定数項のみ、"ct" は定数項+線形トレンドです。この選択を誤ると、正反対の結論に至ります。実験してみましょう。

import numpy as np
from statsmodels.tsa.stattools import adfuller

rng = np.random.default_rng(7)
N = 600
t = np.arange(N)

x_trend = 0.02 * t + rng.normal(0, 1, N)             # トレンド定常
x_drift = np.cumsum(rng.normal(0.02, 1, N))          # ドリフト付き単位根

for name, x in [("トレンド定常", x_trend), ("ドリフト付き単位根", x_drift)]:
    rc = adfuller(x, regression="c", autolag="AIC")
    rt = adfuller(x, regression="ct", autolag="AIC")
    print("%-20s : c  -> 統計量=%7.3f p値=%.4f" % (name, rc[0], rc[1]))
    print("%-20s : ct -> 統計量=%7.3f p値=%.4f" % ("", rt[0], rt[1]))

出力は次のとおりです。

トレンド定常          : c  -> 統計量= -0.164 p値=0.9427
                      : ct -> 統計量=-24.546 p値=0.0000
ドリフト付き単位根    : c  -> 統計量= -2.543 p値=0.1054
                      : ct -> 統計量= -2.535 p値=0.3108

結果を読み解きましょう。トレンド定常な系列 $X_t = 0.02t + \varepsilon_t$ は、決定論的トレンドを除けば白色雑音そのものです。ところが regression="c"(定数項のみ)で検定すると、p値0.9427で「単位根を棄却できない」となってしまいます。回帰式にトレンド項がないため、右肩上がりの動きを「単位根による確率的トレンド」と誤認したのです。regression="ct" に変えてトレンド項を入れると、統計量は $-24.546$、p値は事実上ゼロで、明確に定常と判定されます。

一方、ドリフト付き単位根の系列は、"c" でも "ct" でも棄却できません(p値0.1054と0.3108)。決定論的トレンドを回帰式に入れても、確率的トレンドは取り除けないからです。

この非対称性が実務上の指針を与えてくれます。右肩上がりの系列に対しては "ct" で検定するのが安全です。"ct" で棄却できればトレンド定常なので回帰でトレンドを除去し、棄却できなければ単位根を疑って差分を取る。この使い分けを覚えておけば、前処理の選択で大きく外すことはなくなります。

なぜ定常性を気にするのか — 見せかけの回帰

「定常性を確認しないとどうなるのか」を体感できる、有名な現象があります。見せかけの回帰(spurious regression)です。まったく無関係な2つの非定常系列を回帰すると、統計的に強く有意な関係が「見つかって」しまう、という現象です。

import numpy as np
from scipy import stats

rng = np.random.default_rng(2024)
M, n = 500, 200

cnt_rw, cnt_iid = 0, 0
for _ in range(M):
    a = np.cumsum(rng.normal(0, 1, n))   # 独立なランダムウォークA
    b = np.cumsum(rng.normal(0, 1, n))   # 独立なランダムウォークB
    if stats.linregress(a, b).pvalue < 0.05:
        cnt_rw += 1
for _ in range(M):
    a = rng.normal(0, 1, n)              # 独立なi.i.d.系列A
    b = rng.normal(0, 1, n)              # 独立なi.i.d.系列B
    if stats.linregress(a, b).pvalue < 0.05:
        cnt_iid += 1

print("独立なランダムウォーク同士: p<0.05 となった割合 = %.1f%%" % (100 * cnt_rw / M))
print("独立な i.i.d. 系列同士     : p<0.05 となった割合 = %.1f%%" % (100 * cnt_iid / M))

出力は次のとおりです。

独立なランダムウォーク同士: p<0.05 となった割合 = 82.6%
独立な i.i.d. 系列同士     : p<0.05 となった割合 = 5.6%

見せかけの回帰の図。独立な2本のランダムウォーク、その散布図に引かれた有意な回帰直線、p値のヒストグラムが0付近に集中する様子

左は共通要因が何もない独立な2本のランダムウォーク、中央はその2本を散布図にしたもので、$r = 0.30$、p値 $1.3\times10^{-5}$ という「強く有意な関係」が引かれています。右のp値ヒストグラムが決定的で、本来は $[0,1]$ の一様分布になるはずのp値が、非定常のまま回帰すると0付近にごっそり積み上がる(赤)のに対し、i.i.d.系列同士(青)はきれいに一様分布しています。定常化を怠ると検定の土台そのものが壊れる、ということです。

この数字は衝撃的です。生成過程からして完全に独立な2本のランダムウォークを回帰したところ、82.6%のケースで有意水準5%の有意な関係が検出されました。本来は5%程度であるべき値が、16倍以上に膨れ上がっています。対照として置いたi.i.d.系列同士では5.6%と、理論どおり5%付近に収まりました。

なぜこうなるのでしょうか。ランダムウォークは長期にわたって同じ方向に漂う傾向があるため、2本のランダムウォークを並べると「たまたま両方とも上がっている」「片方が上がって片方が下がっている」といった見かけの連動が高い確率で生じます。回帰の理論はサンプルが独立(あるいは少なくとも定常)であることを前提に標準誤差を計算しているので、この強い持続性を無視した結果、標準誤差が過小評価され、t統計量が過大になるのです。

これは学術的な余興ではありません。「アイスクリームの売上と溺死者数」のような有名な疑似相関の多くは、単に両方が夏に増える(共通のトレンドを持つ)ことに起因します。時系列どうしの関係を論じるなら、まず定常化してから回帰する(あるいは共和分の枠組みを使う)のが鉄則です。定常性の確認は、統計学的なお作法ではなく、誤った結論を防ぐための実質的な防具だと理解してください。

定常性とエルゴード性は別物

最後に、混同されやすい2つの概念を整理しておきます。定常性は「アンサンブルの分布が時間シフトで変わらない」という性質でした。一方エルゴード性は「1本のサンプルパスの時間平均が、アンサンブル平均に一致する」という性質です。

$$ \lim_{T \to \infty} \frac{1}{T} \sum_{t=1}^{T} X_t = E[X_t] = \mu \quad \text{(確率1で)} $$

私たちが実データで平均を計算するとき、暗黙に使っているのはこのエルゴード性のほうです。定常であっても、エルゴード的とは限りません。

反例が明快です。確率変数 $A \sim \mathcal{N}(0,1)$ を一度だけ引き、すべての時刻で

$$ X_t = A \quad (\text{すべての } t) $$

とします。この過程は、どの時刻で切っても分布は $\mathcal{N}(0,1)$、任意の時刻の組の同時分布もシフト不変なので、強定常です。しかし1本のサンプルパスは、たとえば $A = 1.7$ を引いたなら $1.7, 1.7, 1.7, \dots$ という定数列で、何万個平均しても $1.7$ にしかなりません。アンサンブル平均 $0$ には決して収束しない、非エルゴード的な例です。

エルゴード性の十分条件として実用的なのが、自己共分散が減衰することです。広義定常過程について $R(\tau) \to 0$($\tau \to \infty$)、より正確には $\frac{1}{T}\sum_{\tau=0}^{T-1} R(\tau) \to 0$ であれば、標本平均は平均二乗の意味で $\mu$ に収束します(平均に関するエルゴード定理)。AR(1)では $R(\tau) = \phi^\tau \gamma_0$ が指数的に0に落ちるのでエルゴード的、$X_t = A$ では $R(\tau) = 1$ が減衰しないのでエルゴードではない、と綺麗に説明がつきます。

定常性とエルゴード性の違いを示す3パネル図。X_t=Aは強定常だが各パスが定数のまま、時間平均が引いた値に張りつき、自己共分散が減衰しない

左は $X_t = A$ の6本の実現値で、それぞれが引いた値の高さで水平に伸びています。中央が本質で、$X_t = A$ の時間平均(色つき)はサンプル数をいくら増やしても引いた値に張りついたままなのに、AR(1) $\phi=0.5$ の時間平均(黒)はアンサンブル平均0にきちんと収束していきます。右はその理由で、自己共分散が減衰するか否かがエルゴード性を分けていることを示しています。定常であることと、1本のパスから平均を推定できることは、別々の話なのです。

エルゴード性が成り立つとき、標本平均の分散も自己共分散から計算できます。定常過程の標本平均の分散は

$$ \mathrm{Var}\left( \frac{1}{n}\sum_{t=1}^{n} X_t \right) \approx \frac{1}{n} \sum_{\tau=-\infty}^{\infty} R(\tau) $$

と近似され、AR(1)($\mu=0$、$\sigma^2=1$)ではこの総和が $\sigma^2/(1-\phi)^2$ に等しくなります。$\phi = 0.5$、$n = 200$ なら $\frac{1}{(1-0.5)^2 \cdot 200} = 0.0200$ です。実際に2000回のシミュレーションで標本平均の分散を測ると 0.01982 となり、理論値とよく一致します。i.i.d.なら $\gamma_0/n = 1.333/200 = 0.00667$ だったはずですから、正の自己相関がある分だけ標本平均の精度が3倍ほど悪化していることになります。「サンプル数が200ある」と思っていても、正の自己相関があると実質的な情報量はもっと少ない — この考え方は有効サンプルサイズ(ESS)としてエルゴード性とは?時間平均とアンサンブル平均の一致を図解でわかりやすく解説で詳しく扱っています。

まとめると、定常性は「レシピが時間で変わらない」、エルゴード性は「1本のパスがレシピの全貌を語ってくれる」という別々の主張です。実データで統計量を推定するには、両方が必要になります。

まとめ

本記事では、広義定常過程の定義と検証について解説しました。

  • 強定常は「任意個の時刻を選んだ同時分布が、時間シフトで不変」という厳密な条件。分布全体を縛るので実データでの検証は困難
  • 広義定常(弱定常)は「2次モーメント有限・平均一定・自己共分散が時間差のみに依存」という、推定可能な量だけに絞った条件
  • 包含関係は非対称。2次モーメントが有限な強定常 $\Rightarrow$ 広義定常だが、逆は成り立たない。コーシー分布のi.i.d.列(強定常だが広義定常でない)と、偶数時刻が正規・奇数時刻が $\pm1$ の列(広義定常だが強定常でない)が反例
  • ガウス過程では両者が一致する。多変量正規分布が平均ベクトルと共分散行列だけで決まり、シフト量が共分散の引数の引き算で消えるため
  • 自己共分散関数は偶関数であり、$|R(\tau)| \le R(0)$(分散の非負性から導出)を満たし、正定値関数でなければならない。$\rho = (1, 0.9, 0)$ は固有値 $-0.273$ を持ち存在し得ない
  • AR(1)の定常条件は $|\phi| < 1$。分散 $\sigma^2/(1-\phi^2)$ が正であること、MA($\infty$)展開が収束することの両方から同じ条件が出る。$\phi = 1$ のランダムウォークは平均こそ一定だが $\mathrm{Var}(X_t) = t\sigma^2$ が増大するため非定常
  • 検証は3段構え。系列を眺める → 窓ごとの平均・分散・自己相関を比べる → ADF検定。窓平均のレンジ(0.40 対 10.63)と lag-1 自己相関(0.5前後 対 0.9以上)が実用的な判別指標
  • ADF検定は $\Delta X_t = \gamma X_{t-1} + \dots$ の $\gamma = 0$ を帰無仮説とする。帰無仮説が「非定常」なので、棄却できないことは非定常の証明ではない。$\phi = 0.95$ では検出力が44.5%まで落ちる
  • トレンド定常と単位根は regression="ct" で見分ける。トレンド定常な系列は "c" だとp値0.9427で誤判定される
  • 定常化を怠ると見せかけの回帰が起きる。独立なランダムウォーク同士でも82.6%の確率で「有意」が出る
  • 定常性とエルゴード性は別物。$X_t = A$ は強定常だが非エルゴード。自己共分散が減衰することが実用的な十分条件

定常性は、時系列を扱うあらゆる手法の入口に立つ番人です。ここを素通りすると、その先で得られた予測も検定も相関も、すべて足元が崩れます。逆に定常性の確認と適切な前処理さえ済ませてしまえば、AR・ARMA・状態空間モデル・スペクトル解析といった強力な道具がまとめて使えるようになります。

次のステップとして、以下の記事も参考にしてください。