【時系列解析】なぜ，自己共分散（自己相関）が急激に減衰しないとパワースペクトルは定義できないのか？

北川源四郎先生の本「Rによる時系列モデリング入門」(33 ページ)では，

　「ラグ $k$ が大きくなるとき自己共分散関数 $C_k$ が急激に減衰し」

　 $\displaystyle p(f) = \sum_{k=-\infty}^{\infty} C_k e^{- 2 \pi i k f}$

の定義を与える部分に書いてあります (以下では， $\displaystyle p(f)$ を $\displaystyle S(f)$ と書きます)．何で？この条件「自己共分散関数 $C_k$ が急激に減衰」の部分を読み飛ばしてはいけないというのが今回のお話です．

　Wikipediaのウィーナー＝ヒンチンの定理 (Wiener–Khinchin theorem)の説明
ja.wikipedia.org

にはそんなこと書いてないけど，と思うかもしれません．

1次自己回帰過程 ( $a=1$ ， $\sigma^2=1$ )の自己共分散とパワースペクトル

1次自己回帰過程 ( $a=1$ ， $\sigma^2=1$ )の自己共分散とパワースペクトル

自己共分散とパワースペクトルの関係

　時系列 $\{y[i]\}$ の自己共分散関数を，

　 $\begin{eqnarray}C(\, j\, ) = E \left[ y [ i ]\, y[i - j ]\right] \end{eqnarray}$

とします．ここでは， $E$ で期待値を表しています．このとき，次のような有限の長さ $N$ の時系列のパワースペクトルの推定量の期待値 $\bar{S}(f_k)$ を計算してみます．ここでは，時系列のフーリエ変換を使ってパワースペクトルを計算すると (パワースペクトルの定義は複数あります (関連記事1))，

　 $\begin{eqnarray} \bar{S}(f_k) &=& {\rm E} \left[ \frac{\left|Y_k \right|^2}{N} \right] \nonumber \\ &=& {\rm E} \left[ \frac{1}{N} \left(\sum_{n=0}^{N-1} y[n] \, e^{- i 2 \pi f_k n} \right)\left(\sum_{n=0}^{N-1} y[n] \, e^{+ i 2 \pi f_k n} \right) \right] \nonumber \\ &=& \frac{1}{N} \sum_{j=-(N-1)}^{N-1} \left(N-|\, j \,|\right) C(\, j \,)\, e^{- i 2 \pi f_k j} \nonumber \\ &=& \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)\, e^{- i 2 \pi f_k j} \end{eqnarray}$

となります (ちょっと自信がありませんので自分で確かめてください)．途中で，

　 $\displaystyle \frac{1}{N}\sum_{k=0}^{N-1} e^{i 2 \pi n k/N } = \left\{ \begin{array}{ll} \displaystyle 1 & \displaystyle (n = 0) \\ \displaystyle 0 & \displaystyle (n \neq 0) \end{array}\right.$

の公式を使いました．これは，離散の単位インパルス (ディラックじゃありません)

　 $\displaystyle \delta [n] = \left\{ \begin{array}{ll} \displaystyle 1 & \displaystyle (n = 0) \\ \displaystyle 0 & \displaystyle (n \neq 0) \end{array}\right.$

のフーリエ変換が

　 $\displaystyle \sum_{n=0}^{N-1} \delta [n] e^{- i 2 \pi n k/N } = 1$

になることの逆変換です．フーリエ変換と逆フーリエ変換には1対1の関係がありますので片方が分かれば，逆も言えます．

　 $\begin{eqnarray} \bar{S}(f_k) &=& \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)\, e^{- i 2 \pi f_k j} \end{eqnarray}$
となったので，さらに， $N \to \infty$ とすれば，

　 $\begin{eqnarray} S(f_k) &=& \sum_{j=-\infty}^{\infty} C(\, j \,)\, e^{- i 2 \pi f_k j} \end{eqnarray}$

になりそうですが，本当にそうでしょうか．問題は，

　 $\displaystyle \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)$

の部分です． $|\, j \,|$ がすごーく大きいときは， $N$ がそこそこ大きくても，

　 $\displaystyle \left(1- \frac{|\, j \,|}{N}\right) \approx 1$

と近似するには無理があります．

　このような無理をしないために， $|\, j \,|$ がすごーく大きくなる前に， $|\, j \,|$ の増加にともない急速に

　 $\displaystyle C(\, j \,) \approx 0$

となってくれる必要があるのです．つまり，

　「ラグ $k$ が大きくなるとき自己共分散関数 $C_k$ が急激に減衰し」

の条件が必要なのです．とはいえ，なぜ，自己共分散関数 $C_k$ のフーリエ変換でパワースペクトルを定義しているのかについては，私は十分に理解できていません．弱定常過程の前提ありきなのかなと想像しています．どなたか，詳しい方に教えていただきたいです．

練習問題 (すごく重要だと思います．お気に入りの問題です)

　次の確率過程ではない例

　 $\begin{eqnarray} y [i]&=& A \sin \left(2 \pi f_k i \right) \quad (f_k = k/N) \end{eqnarray}$

の自己共分散関数 (減衰せずに振動し続けとる (下図参照))
　 $\begin{eqnarray} C(\, j \,) &=& \frac{A^2}{2} \cos \left(2 \pi f_k j \right) \end{eqnarray}$
を使って，パワースペクトルを計算してみてください．この自己共分散関数は，

　 $\begin{eqnarray} C(\, j \,) &=& \lim_{N \to \infty}\frac{1}{N} \cdot A \sin \left(2 \pi f_k i + \theta \right) \cdot A \sin \left(2 \pi f_k (i+j) + \theta \right) \\ &=& \lim_{N \to \infty} A^2 \left( \frac{1}{2} \cos\left(2 \pi f_k j \right) + \frac{1}{N} \frac{\sin(\cdots) - \sin (\cdots)}{\sin (2 \pi f_k)} \right) \\ &=& \frac{A^2}{2} \cos \left(2 \pi f_k j \right) \end{eqnarray}$
として，無理矢理計算しました． $\theta$ は任意の値です．

正弦波の自己共分散関数とパワースペクトル ( $f_k=50/N$ ， $N=1000$ )

正弦波の自己共分散関数とパワースペクトル ( $f_k=50/N$ ， $N=1000$ )

　長さ $N$ の時系列 $\{ y [i] \}$ のフーリエ変換を使ったパワースペクトルの計算結果 (これは， $N \to \infty$ で収束しないぞ！)は，

　 $\begin{eqnarray} \hat{S}(f_k) = \left\{ \begin{array}{ll} \displaystyle \frac{A^2}{4} N & (k=j, k=N-j) \\ \\ \ 0 & ({\rm otherwise}) \end{array} \right. \end{eqnarray}$

となります (上の図右側)．「えー，なんでNがかかってんの？」と感じること，疑問に思うことが大切です．パワースペクトルというのは，正式にはパワースペクトル密度です．パワースペクトル密度の面積が，時系列の分散に一致するのでした．この場合，長方形で近似することにすれば (三角形でもいいですが)，周波数の幅 $\Delta f=1/N$ を使って，全面積は
　
　 $\displaystyle \frac{A^2}{4} N \cdot \frac{1}{N} + \frac{A^2}{4} N \cdot \frac{1}{N} = \frac{A^2}{2}$

となります．これは，時系列の分散 $\displaystyle C(0)=\frac{A^2}{2}$ と一致します (正しそうです)．

　ということで，上に示したパワースペクトル $\hat{S}(f_k)$ と，自己共分散関数 $C(\, j \,)$ のフーリエ変換 (ウィーナー＝ヒンチンの定理)を使ったパワースペクトルの計算を比較してください (発散しちゃうぞ！)． $C(\, j \,)$ のフーリエ変換が存在しないことに気づいてください．

　ちなみに，上で導いた関係式

　 $\begin{eqnarray} \bar{S}(f_k) &=& \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) C(\, j \,)\, e^{- i 2 \pi f_k j} \end{eqnarray}$

で計算すると，
　 $\begin{eqnarray} \sum_{j=-(N-1)}^{N-1} \left(1- \frac{|\, j \,|}{N}\right) \left\{ \frac{A^2}{2} \cos \left(2 \pi f_k j \right) \right\} e^{- i 2 \pi \frac{k}{N} j} &=& \left\{ \begin{array}{ll} \displaystyle \frac{A^2}{4} N & (k=j, k=N-j) \\ \\ \ 0 & ({\rm otherwise}) \nonumber \\ \end{array} \right. \end{eqnarray}$

となり，上に書いたパワースペクトル $\hat{S}(f_k)$ と一致します．

1. パワースペクトルの定義
chaos-kiyono.hatenablog.com

2. 1次自己回帰過程
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com
chaos-kiyono.hatenablog.com

おまけ

　作図に使ったRスクリプト

sig2 <- 1
a <- 0.9
par(mfrow=c(1,2))
########
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^abs(x),col=4,lwd=2,n=1000,xlim=c(-4*pi/log(a),4*pi/log(a)),xlab="k",ylab="C(k)",main="Autocovariance",xaxs="i")
# 解析的に計算したパワースペクトル
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),col=4,lwd=2,n=1000,xlim=c(-1/2,1/2),xlab="f",ylab="S(f)",main="Power spectral density",xaxs="i")

k <- 50
N <- 1000
A <- 1
par(mfrow=c(1,2))
########
# 解析的に計算した自己共分散関数
curve(1/2*cos(2*pi*k/N*x),xlim=c(-N/k*4,N/k*4),col=4,lwd=2,xlab="k",ylab="C(k)",main="Autocovariance",xaxs="i",n=1000)
# 解析的に計算したパワースペクトル
curve(0*x,col=4,lwd=2,n=100,xlim=c(-1/2,1/2),ylim=c(0,N/3),xlab="f",ylab="S(f)",main="Power spectral density",xaxs="i")
points(c(k/N,(N-k)/N-1),c(N/4*A^2,N/4*A^2),pch=16,col=4)
points(c(k/N,(N-k)/N-1),c(0,0),pch=16,col="white")
points(c(k/N,(N-k)/N-1),c(0,0),pch=1,col=4)