統計学 学習ノート

はじめに

気温は今日から明日を予想できるが、サイコロはできない

「明日の売上は、今日の売上からどこまで予想できるの？」

今日が暑ければ、明日もたぶん暑いですよね。一方で、今日のサイコロの目から、明日の目は予想できません。2級では、この違いを系列相関とコレログラムで読み取りました。でも、読み取るだけでは、予測はできません。

ここまでの記事の手法は、データの順番を入れ替えても結果が変わりませんでした。ところが、時系列では順番そのものが情報です。しかも、「標本」は1本の系列しかありません。歴史は1回しか起きないので、同じ条件で繰り返し観測できないのです。

実は、この難しさを乗り越える鍵が、「性格が時間によって変わらない」という仮定です。性格が同じなら、1本の長い系列の中の別々の時点を、繰り返しの観測の代わりに使えます。

そこで登場するのが、今回の主役 定常性 と ARIMAモデル です。定常性は時系列の性格が一定であることを表し、ARIMAモデルはその性格を数式で書き表します。

前回の記事では、形の違うデータを地図にする手法を学びました。今回の記事では、自己相関、ARIMAモデル、スペクトラム、状態空間モデル、そして回帰の誤差の系列相関を調べるダービンワトソン比について解説します。

定常性と自己相関：時系列の性格を測る

定常なら、平均も波の大きさも時間で変わらない

はじめに、すべての土台である定常性と自己相関について解説します。

定常性（stationarity） とは、時間がたっても統計的な性格が変わらないことです。時系列 $X_t$ が次の3つを満たすとき、弱定常（2次定常）と言います。

【弱定常の3条件】

1. 平均が一定： $E[X_t] = \mu$
2. 分散が一定で有限： $V[X_t] = \gamma_0$
3. 共分散が時間差だけで決まる： $\text{Cov}(X_t, X_{t+k}) = \gamma_k$

波打ち際に例えると、定常な時系列は、毎日同じように寄せては返す波です。波は上下しますが、平均の水位も、波の大きさも変わりません。一方で、満ち潮のように水位が上がり続けるなら、定常ではありません。

例えば、確率過程の記事で学んだランダムウォーク $X_t = X_{t-1} + \varepsilon_t$ は、定常ではありません。$X_0 = 0$ なら $V[X_t] = t\sigma^2$ で、分散が時間とともに増えるからです。トレンドがある系列も、平均が動くので定常ではありません。

平均0、分散 $\sigma^2$ で、時点が違えば無相関な $\varepsilon_t$ を 白色雑音（ホワイトノイズ） と呼びます。予測の手がかりが何もない、最も単純な定常過程です。

定常なら、ラグ $k$ の 自己相関 を次のように定義できます。

$$\rho_k = \frac{\gamma_k}{\gamma_0}$$

2級で学んだ系列相関 $r_k$ は、この $\rho_k$ の推定値です。白色雑音なら、$r_k$ は近似的に平均0、分散 $1/n$ の正規分布に従います。したがって、$|r_k| > 1.96/\sqrt{n}$ なら、有意水準5%で自己相関があると判断します。例えば $n = 100$ なら、境界は ±0.196 です。

次に、偏自己相関 です。ラグ2の自己相関を考えましょう。今日と一昨日の値が似ているのは、「今日が昨日に似ていて、昨日が一昨日に似ている」からだけかもしれません。伝言ゲームで、1人目と3人目の話が似ていても、2人目が伝えただけで、直接は話していないのと同じです。

偏自己相関 $\phi_{kk}$ は、間にある $X_{t-1}, \dots, X_{t-k+1}$ の影響を取り除いた、$X_t$ と $X_{t-k}$ の相関です。偏相関係数の公式に $\rho_1, \rho_2$ を入れると、ラグ2は次の式です。

$$\phi_{22} = \frac{\rho_2 - \rho_1^2}{1 - \rho_1^2}$$

例えば $\rho_1 = 0.8$、$\rho_2 = 0.64$ なら、$\phi_{22} = (0.64 - 0.64) \div 0.36 = 0$ です。ラグ2の相関0.64は、すべて伝言で説明できたことになります。

以上が、定常性と自己相関についての解説です。

ARIMAモデル：過去の値と過去のショックで書く

ARは偏自己相関が、MAは自己相関がラグの後で切れる

次に、時系列を数式で書くARIMAモデルについて説明します。

ARIMAモデルは、3つの部品からできています。1つずつ見ていきましょう。

1. AR：過去の自分を引きずる

AR（自己回帰）モデル は、今の値を過去の値で回帰するモデルです。AR(1)は次の形です。

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

お風呂の湯の温度に例えると、今の温度は「少し前の温度の一部 + 新しいゆらぎ」です。$|\phi| < 1$ なら過去の影響は薄れて定常になり、$\phi = 1$ ならランダムウォークです。

定常のとき、両辺の分散をとると $\gamma_0 = \phi^2\gamma_0 + \sigma^2$ なので、分散が求まります。また、両辺に $X_{t-k}$ をかけて期待値をとると、$\gamma_k = \phi\gamma_{k-1}$ です。

$$\gamma_0 = \frac{\sigma^2}{1 - \phi^2}, \qquad \rho_k = \phi^k$$

例えば $\phi = 0.8$ なら、自己相関は 0.8, 0.64, 0.512 と、等比数列で減っていきます。

AR(2) $X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \varepsilon_t$ では、同じ操作で ユール・ウォーカー方程式 が得られます。

$$\rho_k = \phi_1\rho_{k-1} + \phi_2\rho_{k-2} \quad (k \ge 1)$$

一般のAR($p$)が定常である条件は、特性方程式 $1 - \phi_1 z - \dots - \phi_p z^p = 0$ の解の絶対値がすべて1より大きいことです。

2. MA：過去のショックの余韻

MA（移動平均）モデル は、今の値を過去のショックで書くモデルです。MA(1)は次の形です。

$$X_t = \varepsilon_t + \theta\varepsilon_{t-1}$$

鐘の音に例えると、今日打った音と、昨日の音の余韻だけが聞こえる状態です。一昨日の音は完全に消えています。分散は $(1 + \theta^2)\sigma^2$、ラグ1の共分散は $\theta\sigma^2$ で、ラグ2以降の共分散は0です。

$$\rho_1 = \frac{\theta}{1 + \theta^2}, \qquad \rho_k = 0 \quad (k \ge 2)$$

例えば $\theta = 0.5$ なら $\rho_1 = 0.5 \div 1.25 = 0.4$ です。$\theta = 2$ でも同じ $\rho_1 = 0.4$ になるので、自己相関だけでは $\theta$ を決められません。そこで、$|\theta| < 1$ の側を選びます。これを 反転可能性 と言います。なお、$|\rho_1|$ は0.5を超えません。

ARとMAを合わせたものが、ARMA($p, q$)モデル です。モデルの見分け方は、次の表のとおりです。

モデル	自己相関	偏自己相関
AR($p$)	徐々に減衰する	ラグ $p$ の後で切れる
MA($q$)	ラグ $q$ の後で切れる	徐々に減衰する
ARMA($p, q$)	徐々に減衰する	徐々に減衰する

ここがポイントです。「切れる」と「減衰する」が、ARとMAで入れ替わります。AR(1)の偏自己相関はラグ1の後で0で、MA(1)の自己相関もラグ1の後で0です。

3. I：階差で定常に直す

現実の系列の多くは、トレンドがあって定常ではありません。そこで、隣との差をとります。これを 階差（差分）と言います。

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

例えば 10, 12, 15, 19, 24 の階差は 2, 3, 4, 5 で、さらに階差をとると 1, 1, 1 です。1回の階差で直線のトレンドが消え、2回の階差で2次曲線のトレンドが消えます。

$d$ 回の階差をとった系列が ARMA($p, q$) に従うとき、元の系列は ARIMA($p, d, q$)モデル に従うと言います。I は integrated（和分）の略です。例えば、ランダムウォークは ARIMA(0, 1, 0) です。

月次データの季節変動には、12か月前との差 $X_t - X_{t-12}$ をとる 季節階差 を使います。対数をとってから階差をとると、2級で学んだ成長率に近い量になります。

【ARIMAモデルの当てはめ方（ボックス・ジェンキンス法）】

1. 同定： グラフと自己相関を見て、階差の回数 $d$ と次数 $p, q$ を選ぶ
2. 推定： 最尤法などで係数を推定する
3. 診断： 残差が白色雑音に見えるか確かめる
4. 予測： 残差に問題がなければ予測する

例題：AR(2)の自己相関とAR(1)の予測

(1) $X_t = 0.5 X_{t-1} + 0.3 X_{t-2} + \varepsilon_t$ の $\rho_1$、$\rho_2$ を求めましょう。

【解き方】

1. $k = 1$ を入れる： $\rho_1 = 0.5 + 0.3\rho_1$（$\rho_0 = 1$、$\rho_{-1} = \rho_1$ を使う）
2. 解く： $\rho_1 = 0.5 \div 0.7 \approx 0.714$
3. $k = 2$ を入れる： $\rho_2 = 0.5 \times 0.714 + 0.3 \approx 0.657$

(2) 平均50のAR(1) $X_t - 50 = 0.8(X_{t-1} - 50) + \varepsilon_t$ で、最新の値が60でした。1期先と2期先の予測値を求めましょう。

【解き方】

1. 1期先： 平均からのズレ10が0.8倍になるので、50 + 8 = 58
2. 2期先： ズレがさらに0.8倍になるので、50 + 6.4 = 56.4

答えは58と56.4です。遠い将来の予測は、平均50に近づきます。予測の誤差の分散は、1期先が $\sigma^2$、2期先が $(1 + 0.8^2)\sigma^2 = 1.64\sigma^2$ で、先へ行くほど広がります。

試験のポイント：自己相関と偏自己相関の形から、モデルの種類を答える問題がよく出ます。「切れるのはどちらか」を確認しましょう。

以上が、ARIMAモデルについての説明です。

スペクトラムとペリオドグラム：周期に分解する

波を周波数に分けると、周期4のピークが現れる

続いて、時系列を周期の成分に分ける見方について解説します。

ここまでは、時間の順に時系列を見てきました。今度は見方を変えます。太陽の光をプリズムに通すと、色ごとに分かれますよね。同じように、時系列を「速い波」や「ゆっくりした波」に分けると、どの周期の波がどれだけ強いかが分かります。

スペクトラム（スペクトル密度） は、時系列の分散を、周波数ごとに分けて表した関数です。自己共分散から、次のように定義します。

$$f(\omega) = \frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\gamma_k\, e^{-ik\omega} \qquad (-\pi \le \omega \le \pi)$$

$\omega$ は角周波数で、周期とは「周期 = $2\pi/\omega$」の関係にあります。自己共分散とスペクトラムは、フーリエ変換でつながった表と裏です。積分すると $\int_{-\pi}^{\pi} f(\omega)\,d\omega = \gamma_0$ で、分散に戻ります。つまり、スペクトラムは分散の内訳です。

いくつかの例を見てみましょう。

【スペクトラムの例】

• 白色雑音： $f(\omega) = \sigma^2/(2\pi)$。すべての周波数が同じ強さで、平らになる
• AR(1)： $f(\omega) = \dfrac{\sigma^2}{2\pi(1 - 2\phi\cos\omega + \phi^2)}$

白色雑音の名前は、すべての色を同じ強さで含む光が白いことから来ています。AR(1)で $\phi = 0.8$ なら、$\omega = 0$ の強さは $\omega = \pi$ の81倍です。正の自己相関をもつ系列は、ゆっくりした波が主役なのです。

実際のデータから計算する量が、ペリオドグラム です。周波数 $\omega_j = 2\pi j/n$ ごとに、次の量を求めます。

$$I(\omega_j) = \frac{1}{n}\left|\sum_{t=1}^{n} x_t\, e^{-i\omega_j t}\right|^2$$

データ $x_t$ は平均を引いておきます。言葉で言うと、データとその周波数の波を重ね合わせて、どれだけ揃っているかを測っています。周期が合う周波数で、値が大きくなります。なお、定数の $1/n$ を $1/(2\pi n)$ とする流儀もあります。

例題：ペリオドグラムのピークを読む

データ 1, 0, −1, 0, 1, 0, −1, 0（$n = 8$）のペリオドグラムを、$j = 1, 2$ について求めましょう。

【解き方】

1. $j = 2$ のとき $\omega_2 = \pi/2$ で、$\sum_t x_t e^{-i\pi t/2}$ を計算する。0でない4項がどれも $-i$ になり、和は $-4i$
2. 絶対値の2乗を $n$ で割る： $I(\omega_2) = 16 \div 8 = 2$
3. $j = 1$ のとき $\omega_1 = \pi/4$ で、4項の和は打ち消し合って0。したがって $I(\omega_1) = 0$

答えは $I(\omega_1) = 0$、$I(\omega_2) = 2$ です。ピークは $j = 2$ にあり、周期は $n/j = 8 \div 2 = 4$ です。データを見ると、確かに4期ごとに同じ値を繰り返しています。

読み方の例も見てみましょう。120か月の月次データで $j = 10$ にピークがあれば、周期は 120 ÷ 10 = 12 か月で、季節変動があると分かります。

ただし、ペリオドグラムは、データを増やしてもギザギザが収まりません。各周波数の値の分散が、$n$ が増えても小さくならないのです。つまり、スペクトラムの一致推定量ではありません。そこで、隣り合う周波数の値を平均して、滑らかにしてから使います。

以上が、スペクトラムとペリオドグラムについての解説です。

状態空間モデル：見えない状態を追いかける

上下する読みの裏の、見えない状態を追いかける

次に、観測の裏にある状態を追いかける状態空間モデルについて説明します。

毎朝、体重計に乗るとしましょう。読みは 60.4kg、59.8kg、60.6kg と上下します。でも、1日で本当の体重がそんなに変わるはずはありません。読みには、水分や計り方の誤差が入っています。本当に知りたいのは、見えない「本当の体重」です。

状態空間モデル は、見えない 状態 と、それを誤差つきで観測した値を、2本の式で表すモデルです。最も単純な ローカルレベルモデル は、次の形です。

$$\text{観測方程式}\quad y_t = \mu_t + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \sigma_\varepsilon^2)$$ $$\text{状態方程式}\quad \mu_t = \mu_{t-1} + \eta_t, \qquad \eta_t \sim N(0, \sigma_\eta^2)$$

$\mu_t$ が見えない本当の体重で、ランダムウォークのようにゆっくり動きます。$y_t$ が体重計の読みです。

状態の推定に使うのが、カルマンフィルタ です。「予測」と「更新」を、1時点ごとに繰り返します。$a_t$ を状態の推定値、$P_t$ をその分散とします。

【カルマンフィルタの手順】

1. 予測： 状態の予測は $a_{t-1}$ のままで、分散は $P_{t-1} + \sigma_\eta^2$ に増える
2. カルマンゲイン： $K_t = \dfrac{P_{t-1} + \sigma_\eta^2}{P_{t-1} + \sigma_\eta^2 + \sigma_\varepsilon^2}$
3. 更新： $a_t = a_{t-1} + K_t(y_t - a_{t-1})$、$P_t = (1 - K_t)(P_{t-1} + \sigma_\eta^2)$

言葉で言うと、「新しい推定値 = 予測 + ゲイン × 予測の外れ」です。ゲイン $K_t$ は0から1の間の値で、自分の予測と新しい観測のどちらを信じるかの配分です。体重計の誤差が大きいと、$K_t$ は小さくなります。一言で言うと、信用度に応じた加重平均です。

例題：カルマンフィルタを手で回す

$\sigma_\eta^2 = 1$、$\sigma_\varepsilon^2 = 4$ のローカルレベルモデルで、$a_0 = 20$、$P_0 = 3$ とします。$y_1 = 26$、$y_2 = 22$ を観測したときの $a_1, a_2$ を求めましょう。

【解き方】

1. 予測の分散： $3 + 1 = 4$
2. ゲイン： $K_1 = 4 \div (4 + 4) = 0.5$
3. 更新： $a_1 = 20 + 0.5 \times (26 - 20) = 23$、$P_1 = 0.5 \times 4 = 2$
4. 次の予測の分散： $2 + 1 = 3$
5. ゲイン： $K_2 = 3 \div (3 + 4) \approx 0.429$
6. 更新： $a_2 = 23 + 0.429 \times (22 - 23) \approx 22.57$

答えは $a_1 = 23$、$a_2 \approx 22.57$ です。観測が26まで跳ねても、推定は23までしか動きません。

状態空間モデルの強みは、次のとおりです。

【状態空間モデルの強み】

1. 欠測に強い： 観測がない時点は、更新を飛ばして予測だけを行えばよい
2. 組み立てやすい： トレンド、季節、回帰の項を、状態の成分として足していける
3. 幅が広い： ARIMAモデルも状態空間モデルの形に書き直せる

補足（ここは読み飛ばしてOK）： ローカルレベルモデルでは、時間がたつと $K_t$ は一定の値に近づきます。すると、更新の式は指数平滑法と同じ形になります。

以上が、状態空間モデルについての説明です。

系列相関とダービンワトソン比：回帰の誤差を点検する

残差がゆっくり動くと、DW比は0に近づく

最後に、時系列データの回帰で必ず確かめる点を見ていきましょう。

重回帰分析の記事で、最小二乗法の標準誤差は、誤差が互いに無相関であることを前提にしていました。ところが、時系列データの回帰では、残差がプラスの時期とマイナスの時期が続くことがよくあります。この状態が、誤差の 系列相関 です。

系列相関があると、どうなるのでしょうか？ 説明変数が誤差と無相関なら、係数の推定値に偏りは生じません。でも、標準誤差が正しくなくなります。正の系列相関があると、標準誤差は小さく出やすく、$t$ 値が大きくなり、効果がないのに有意に見えることがあります。答案の採点で、同じ子のテストを別人のものとして数えてしまうようなものです。情報の量を水増ししているのです。

調べる道具が、ダービンワトソン比（DW比） です。回帰の残差 $e_t$ から、次の量を計算します。

$$DW = \frac{\sum_{t=2}^{n}(e_t - e_{t-1})^2}{\sum_{t=1}^{n} e_t^2}$$

分子を展開すると、$n$ が大きいとき、残差のラグ1の自己相関 $r_1$ を使って次のように近似できます。

$$DW \approx 2(1 - r_1)$$

$r_1$ は −1 から1の間なので、DW比は0から4の間の値です。

【DW比の読み方】

• 2に近い： 系列相関なし
• 0に近い： 正の系列相関（残差がゆっくり動く）
• 4に近い： 負の系列相関（残差が毎回符号を変える）

厳密な検定では、標本の大きさと説明変数の数で決まる2つの値 $d_L < d_U$ を、表から読みます。正の系列相関に対する判定は、次のとおりです。

【DW検定の判定】

1. $DW < d_L$： 正の系列相関がある
2. $d_L \le DW \le d_U$： 判定を保留する
3. $DW > d_U$： 正の系列相関があるとは言えない

判定保留の区間があるのが、この検定の特徴です。

例題：DW比を計算する

ある回帰の残差が、1, 2, 1, −1, −2, −2, −1, 0, 1, 1 でした。DW比を求めましょう。表の値は $d_L = 0.88$、$d_U = 1.32$ だとします。

【解き方】

1. 隣との差： 1, −1, −2, −1, 0, 1, 1, 1, 0
2. 分子（差の2乗の和）： 1 + 1 + 4 + 1 + 0 + 1 + 1 + 1 + 0 = 10
3. 分母（残差の2乗の和）： 1 + 4 + 1 + 1 + 4 + 4 + 1 + 0 + 1 + 1 = 18
4. 割り算： $DW = 10 \div 18 \approx 0.56$

答えは約0.56で、$d_L$ より小さいので、正の系列相関があると判断します。実際、$r_1 = 12 \div 18 \approx 0.67$ で、近似式の値は 2 × 0.33 ≈ 0.67 です。

系列相関が見つかったら、誤差のAR(1)構造を入れた一般化最小二乗推定を使うか、変数の不足を見直します。

補足（ここは読み飛ばしてOK）： 無関係な2本のランダムウォークを回帰すると、決定係数が高く、$t$ 値も有意に見えることがあります。これを見せかけの回帰と言います。このときDW比は0に近い値になるので、DW比は警報の役割も果たします。

試験のポイント：$DW \approx 2(1 - r_1)$ を覚えておくと、DW比から $r_1$ の値を読めます。「2なら無相関、小さいと正、大きいと負」です。

以上が、系列相関とダービンワトソン比についての解説です。

細かい式を忘れてしまっても問題ありません。ここでは、「定常性を確かめ、自己相関の形でモデルを選び、残差が白色雑音かを点検する」という点だけ押さえておきましょう。

今回のまとめ

定常か見て、モデルを選び、残差が白色雑音かを点検する

それでは、今回の内容を振り返ります。

今回は、時間の順番そのものを情報として使う時系列解析について紹介しました。時間の領域で見れば自己相関とARIMAモデル、周波数の領域で見ればスペクトラム、見えない状態に注目すれば状態空間モデルです。入口は違っても、「過去の情報から、どこまでを説明できるか」を問う点は同じです。

今回のポイントを振り返りましょう。

振り返り

• 定常性は性格が変わらないこと。平均・分散が一定で、共分散はラグだけで決まる。
• AR(1)の自己相関は $\phi^k$ で減衰し、MA(1)の自己相関はラグ1の後で切れる。偏自己相関では、これが入れ替わる。
• 定常でない系列は、階差で定常に直す。ARIMA($p, d, q$) の $d$ は階差の回数。
• スペクトラムは分散の周波数ごとの内訳。ペリオドグラムのピークから周期を読む。
• 状態空間モデルは体重計と本当の体重。カルマンフィルタは予測と観測の加重平均を繰り返す。
• ダービンワトソン比は約 $2(1 - r_1)$ で、回帰の残差の系列相関を調べる。

キーワードのおさらい

• 自己相関: 同じ系列の、ラグ $k$ 離れた値どうしの相関
• 偏自己相関: 間の時点の影響を取り除いた自己相関
• スペクトラム: 時系列の分散を周波数ごとに表した関数
• ペリオドグラム: データから計算する、周波数ごとの強さ
• ARIMAモデル: 階差をとった系列に、ARとMAのモデルを当てはめるモデル
• 定常性: 平均・分散・自己共分散の構造が時間によって変わらないこと
• 階差: 隣り合う時点の値の差をとること
• 状態空間モデル: 見えない状態と観測を、2本の式で表すモデル
• 系列相関: 時点の違う値どうしの相関。回帰では誤差の相関が問題になる
• ダービンワトソン比: 回帰の残差のラグ1の系列相関を調べる統計量

さて、ここまでは量的なデータの分析でした。話をクロス表に戻しましょう。2級で学んだカイ二乗検定は、度数が少ないと使えませんでした。関連の「強さ」を測る数も、まだ手元にありません。

次回は、オッズ比や正確検定など「分割表の解析」について紹介します。

それでは、次の記事へ進んでいきましょう。