R Time Series Analysis 時系列解析(13) スペクトル分析

定常過程の母スペクトル

定常時系列を特徴を把握するための要素として、自己共分散 $\gamma_{k},k=0,1,2,\cdots,$ およびそれを基準化した自己相関係数 $\rho_{k},k=0,1,2,\cdots$ があった。定常過程の自己共分散に対して、次の無限級数を「自己共分散母関数」として定義する。

$g_{y}(z) \equiv \sum_{k=-\infty}^{\infty} \gamma_{k} z^k$

これを $2 \pi$ で割り、 $z$ に $e^{-i \omega}$ （ただし、 $\omega$ は周波数)を代入したものは系列 ${\gamma}$ のフーリエ変換となっている。(これを $y$ の母スペクトルと呼んでいる。)　※フーリエ変換はココがわかりやすいと思った。

$f(\omega) = \frac{1}{2\pi} \sum_{k=-\infty}^{\infty} \gamma_{k} e^{-i \omega k}$

$= \frac{1}{2\pi} \sum_{k=-\infty}^{\infty} \gamma_{k} \cos \omega k$

$= \frac{1}{2\pi} \gamma + \frac{1}{\pi} \sum_{k=1}^{\infty} \gamma_{k} \cos \omega k,-\pi \leq \omega \leq \pi$

つまり、自己共分散 $\gamma_{k}$ は、関数 $f(\omega)$ のフーリエ逆変換となっており、

$\gamma_{k} =\int_{-\infty}^{\infty} f(\omega) e^{-i \omega k} d \omega$

と表せる。

さまざまなモデルの母スペクトル

ここまで登場した様々なモデルに対して、母スペクトルが求められる。まず、線形定常過程をMA(∞)で表現すると、

$y_{t} = \sum_{j=0}^{\infty} \psi a_{t-j} = \psi (B) a_{t}, \psi_{0}=1$

自己共分散母関数は次のように表せることが知られている（要確認）。

$g_{y} (z) = \sigma_{a}^2 \psi (z) \psi (z^{-1})$

これを $2 \pi$ で割り、 $z$ に $e^{-i \omega}$ を代入したものは、 $\gamma_{k}$ のフーリエ変換 $f(\omega)$ となっている。

$f(\omega) = \frac{1}{2\pi} g_{y} (e^{-i \omega})$

$= \frac{1}{2\pi} \sigma_{a}^2 \psi(e^{-i \omega}) \psi(e^{i \omega})$

一方、ARMA(p,q)モデルは次のように表せる。

$\phi_{p} (B) y_{t} = \theta_{q} (B) a_{t}$

$\phi_{p} (B) = 1 - \phi_{1} B - \phi_{2} B^2 -\cdots - \phi_{p} B^p$

$\theta_{p} (B) = 1 + \theta_{1} B + \theta_{2} B^2 + \cdots + \theta_{p} B^p$

ただし、 $B y_{t} = y_{t-1}, B a_{t} = a_{t-1}$

これを変形して、

$y_{t} = \psi (B) a_{t} = \frac{\theta_{q} (B)}{\phi_{p} (B)} a_{t}$

よって、自己共分散系列 ${\gamma_{k}}$ のフーリエ変換 $f(\omega)$ は、

$f(\omega) = \frac{\sigma_{a}^{2}}{2 \pi} \frac{\theta_{q} (e^{-i \omega}) \theta_{q} (e^{i \omega})}{\theta_{q} (e^{-i \phi}) \theta_{q} (e^{i \phi})}$

$=\frac{\sigma_{a}^{2}}{2 \pi} |\frac{\theta (e^{-i \omega})}{\phi (e^{-i \omega})}|^2$

$=\frac{\sigma_{a}^{2}}{2 \pi} {(1+\theta_1 e^{-i \omega}+\cdots+\theta_{q} e^{-iq \omega})(1+\theta_{1} e^{i \omega}+\cdots+\theta_{q} e^{iq \omega})}{(1-\phi_1 e^{-i \omega}+\cdots+\phi_{q} e^{-iq \omega})(1+\phi_{1} e^{i \omega}+\cdots+\phi_{q} e^{iq \omega})}$

上記に倣えば、様々なモデルの母スペクトルを求められる。

ホワイトノイズの場合、 $y_{t} = a_{t}$ なので、

$f(\omega) = \frac{\sigma_{a}^2}{2\pi}$

AR(1)の場合、定常性の条件 $|\phi|$ <1を満たせば、 $y_{t} = (1 - \phi B)^{-1} a_{t}$ なので、

$f(\omega) = \frac{\sigma_{a}^2}{2\pi} \frac{1}{(1-\phi e^{-i \omega})(1 - \phi e^{i \omega})}$

$= \frac{\sigma_{a}^2}{2\pi} \frac{1}{(1+\phi^{2}-2\phi \cos \omega)}$

AR(2)の場合も、定常性の条件を満たせば(特性方程式の解の絶対値が1より小さければ)、

$y_{t} = (1 - \phi_{1} B - \phi_{2} B^2)^{-1} a_{t} = (1-\lambda_{1} B)^{-1}(1 - \lambda_{2} B)^{-1} a_{t}$ なので、

$f(\omega) = \frac{\sigma_{a}^2}{2\pi} \frac{1}{(1 - \lambda_{1}e^{-i \omega})(1 - \lambda_{2}e^{-i \omega})(1 - \lambda_{1}e^{i \omega})(1 - \lambda_{2}e^{i \omega})}$

$=\frac{\sigma_{a}^2}{2\pi} \frac{1}{(1 + \lambda_{1}^{2} - 2 \lambda_{1} \cos \omega)(1 + \lambda_{2}^{2} - 2 \lambda_{2} \cos \omega)}$

のように表せる。MA(1)やMA(2)、ARMA(1,1)でも同様。

母スペクトルでわかること

ここによると、「Aのフーリエ変換」とは「Aの周波数成分の分布を与える変換」という意味である(と理解している)。したがって、自己共分散系列 ${\gamma_{k}}$ のフーリエ変換はその周波数の分解を与えている。

$\omega$ の値が増加するにつれて、 $f(\omega)$ が増加するのであれば、周波数の大きい成分の方が自己共分散系列への影響が大きい
$\omega$ の値が減少するにつれて、 $f(\omega)$ が減少するのであれば、周波数の小さい成分の方が自己共分散系列への影響が大きい

と解釈が出来る。つまり、フーリエ変換の側から見ることで、ある系列の周波数成分の影響を確認することが出来る。

【疑問点】

読んだ本では、時系列そのものではなく、自己共分散系列のフーリエ変換を考察していたが、そもそも興味の対象は時系列そのものであるので、そちらのフーリエ変換を考えるべきだと感じた。わざわざ自己共分散系列について考える必要があるのであろうか。

Rで母スペクトルを描く

Rで母スペクトルの挙動を図示してみる。

#AR(1)モデルの母スペクトル
ar1 <- function(sig2,phi) {
  w <- seq(0,pi,0.005*pi)
  y <- (sig2/(2*pi))*(1/(1+phi^2-2*phi*cos(w)))
  ttl <- paste("AR(1) : φ = ",phi)
  yl <- c(0,max(y)*1.05)
  plot(w,y,type="l",ylim=yl,main=ttl,axes=FALSE)
  lbl <- c("0","0.2pi","0.4pi","0.6pi","0.8pi","pi")
  axis(1,at=seq(0,pi,0.2*pi),label=lbl)
  axis(2)
}
#AR(2)モデルの母スペクトル
ar2 <- function(sgi2,phi1,phi2) {
  w <- seq(0,pi,0.005*pi)
  y1 <- (1+phi1^2+phi2^2)-2*phi1*(1-phi1)*cos(w)
  y2 <- -2*phi2*cos(2*w)
  y <- (sig2/(2*pi))*(1/(y1+y2))
  ttl <- paste("AR(2) : φ1 = ",phi1," φ2 = ",phi2)
  yl <- c(0,max(y)*1.05)
  plot(w,y,type="l",ylim=yl,main=ttl,axes=FALSE)
  lbl <- c("0","0.2pi","0.4pi","0.6pi","0.8pi","pi")
  axis(1,at=seq(0,pi,0.2*pi),label=lbl)
  axis(2)
}

#図示
par(mfrow=c(3,2))
ar1(1.2,0.8)
ar1(1.2,0.3)
ar1(1.2,-0.5)
ar2(1.2,0.5,0.2)
ar2(1.2,0.8,0.2)
ar2(1.2,0.5,-0.3)
par(mfrow=c(1,1))

結果がこれ。

f:id:nakhirot:20130714151131p:plain

AR(1)では係数φの正負で周波数成分の影響の大きさの傾向が異なることがわかる。また、AR(2)では係数φ1, φ2の絶対値が周波数成分の影響の大きさに影響しているように見える。

Spectrum関数で標本スペクトルを求める

Rには標本時系列[TeX:y]の標本のスペクトルを計算し、グラフに表示するspectrum関数というものがあるとのこと。引数methodの設定がポイント。

method="pgram"とすると、ペリオドグラムから標本のスペクトルを計算する。
method="ar"とすると、データにARモデルを当てはめて標本のスペクトルを計算する。(つまりは上記で行った計算プロセス)

ペリオドグラムとは、時系列[TeX:y_k]をフーリエ級数展開

[TeX:y_t = \sum_{k=0}^{[T/2] (a_k \cos \omega_k t + b_k \sin \omega_k t,(t=1,2,\cdots,T),\omega_k = \frac{2 \pi k}{T}]

したときに、各周波数成分[TeX:\omega_{k}]の強さを表すもので、次の通り。

[TeX: \begin{eqnarray}I(\omega_k)=\left\{ \begin{array}{ll}T a_{0}^2 & (k=0) \\ \frac{T}{2} (a_{k}^{2} + b_{k}^{2}) & (k=1,2,\cdots,[(T-1)/2]) \\ T a_{T/2}^{2} & (k = T/2 (When T is even) \\\end{array} \right.\end{eqnarray}]

実際に、サンプルでAR(1), AR(2)の系列を発生させて、そのスペクトルを求めてみる。なお、method="pgram"とした場合、ペリオドグラムは変動が大きいことから、平均移動の一種であるmodified Daniell smoothersで平滑化するためのオプションspansを指定する。（あとで調べよう）

#AR(1) φ=0.8 系列を作る
yar108 <- arima.sim(n=100,list(ar=0.8),innov=rnorm(100,sd=1.4))

#AR(2) φ1=0.5,φ2=-0.3系列を作る
yar2 <- arima.sim(n=100,list(ar=c(0.5,-0.3)),innov=rnorm(100,sd=1.4))

#spectrum()関数によるグラフを描く。
par(mfrow=c(2,2))
spectrum(yar108,method="pgram")
spectrum(yar108,method="pgram",spans=c(5,5))
spectrum(yar108,method="ar")
spectrum(yar2,method="pgram",spans=c(5,5))
par(mfrow=c(1,1))

f:id:nakhirot:20130715153728p:plain

右側の青い線が気になるが、、、そのうち調べよう。