ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

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

北川 源四郎先生の本「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)の自己共分散とパワースペクトル

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

 時系列\{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/NN=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)と一致します.

おまけ

 作図に使った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)