【時系列解析】２次自己回帰過程のパワースペクトル（解析的導出）

今回から，本格的に２次自己回帰過程のパワースペクトルに入ります．自己回帰過程の理解において，最も重要なポイントの一つは

「n次自己回帰過程のパワースペクトルは，1次自己回帰過程と2次自己回帰過程のパワースペクトルの和になっている」

ということです (ほとんどの教科書には書いてありませんが，というか，書いてある教科書を私は知りません)．ですので，1次と2次自己回帰過程について理解すれば，自己回帰過程を完全にマスターした気分に近づけると思います．

今回は，前回紹介した，以下のテクニックを使います．
chaos-kiyono.hatenablog.com

２次自己回帰過程のパワースペクトル

それでは，２次自己回帰過程

　 $y[n] = a_1 y[n-1] +a_2 y[n-2] + x[n]$

を考えます．ここで， $a_1$ ， $a_2$ は実数のパラメタ， $\{x[n]\}$ は平均0，分散 $\sigma^2$ の白色ノイズです（以前はwとしたのを，xにしました）．

$y [n]$ のフーリエ変換を， $Y(f_k)$ ，パワースペクトルの推定量を $S_y(f_k)$ とします． $x [n]$ についても，それぞれ， $X(f_k)$ ， $S_x(f_k)$ とします．そして， $S_x(f_k)=\sigma^2$ です (これは真のパワースペクトルですが)．Nを時系列の長さとすれば， $f_k=k/N$ です．上の式の両辺をフーリエ変換して，両辺の絶対値をとって２乗して，Nで割ると以下のようになります．
　 $\begin{eqnarray} Y(f_k) &=& a_1 Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-1} + a_2 Y(f_k) \left( e^{{\rm \bf i} 2 \pi f_k} \right)^{-2} + X(f_k) \\ \left(1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right)Y(f_k) &=& X(f_k) \\ \left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2 \frac{ \left|Y(f_k)\right|^2}{N} &=& \frac{\left|X(f_k)\right|^2}{N} \\ \left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2_y(f_k) &=& S_x(f_k) \\ S_y(f_k) &=& \frac{S_x(f_k)}{\left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2} \\ S_y(f_k) &=& \frac{\sigma^2}{\left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2} \end{eqnarray}$

複素数の絶対値の２乗が， $|z|^2 = z \bar{z}$ であることに注意して，オイラーの公式

　 $e^{{\rm \bf i} \theta} = \cos \theta + {\rm \bf i} \sin \theta$

を使えば，２次自己回帰過程のパワースペクトル

　 $\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{1 + a_1^2 + a_2^2 - 2 a_1 (1 - a_2) \cos(2 \pi f_k) - 2 a_2 \cos( 4 \pi f_k)} \end{eqnarray}$

がえられます．ここから，このパワースペクトルの性質を考えていきます．

パワースペクトルの例． $a_1 = 0.9 \sqrt{3}$ ， $a_2 = -0.81$ ．

パワースペクトルの例． $a_1 = 0.9 \sqrt{3}$ ， $a_2 = -0.81$ ．

パワースペクトルが1次自己回帰過程のパワースペクトルの和になる場合

２次自己回帰過程のパワースペクトルには，振動に対応するピークがあると思い込んでいませんか．違うときもあります．

　 $a_1^2+4 a_2 \ge 0$ のときは，振動しません．
（注意： $f=1/2$ のプラスマイナスを繰り返す振動はありえます）

ここでは，
　 $\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{\left|1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 e^{-2 {\rm \bf i} 2 \pi f_k} \right|^2} \\ &=& \frac{\sigma^2}{\left(1 - a_1 e^{-{\rm \bf i} 2 \pi f_k} - a_2 \left(e^{-{\rm \bf i} 2 \pi f_k}\right)^2 \right)\left(1 - a_1 e^{{\rm \bf i} 2 \pi f_k} - a_2 \left(e^{{\rm \bf i} 2 \pi f_k}\right)^2 \right) } \end{eqnarray}$

の分母に注目します． $e^{{\rm \bf i} 2 \pi f_k}$ の部分を $z$ で置き換えた式

$1 - a_1 z - a_2 z^2$

が実数の範囲で因数分解できて，

$1 - a_1 z - a_2 z^2 = \left(1- \frac{z}{\lambda_1} \right)\left(1- \frac{z}{\lambda_2} \right)$

とできるとします． $a_1^2+4 a_2 > 0$ が仮定されています (重解はやりません)． $\lambda_1$ と $\lambda_2$ は， $1 - a_1 z - a_2 z^2 = 0$ の解で

$\lambda_1 = \displaystyle \frac{-a_1 + \sqrt{a_1^2+4 a_2} }{2 a_2}$ ， $\lambda_2 = \displaystyle \frac{-a_1 - \sqrt{a_1^2+4 a_2} }{2 a_2}$

です．このとき，上の $S_y(f_k)$ の式は， $\lambda_1$ と $\lambda_2$ を使って，
　 $\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right) \left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right)} \\ &=& \frac{\sigma^2}{\left\{\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_1} \right)\right\} \left\{ \left(1 - \frac{ e^{{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right)\left(1 - \frac{ e^{-{\rm \bf i} 2 \pi f_k}}{\lambda_2} \right)\right\} } \\ &=& \frac{\sigma^2 \lambda_1^2 \lambda_2^2}{\left(1 + \lambda_1^2 - 2 \lambda_1 \cos \left(2 \pi f_k \right) \right) \left(1 + \lambda_2^2 - 2 \lambda_2 \cos \left(2 \pi f_k \right) \right)} \\ &=& \frac{\lambda_1^2 \lambda_2^3 \sigma^2}{\left(\lambda_1^2 \lambda_2-\lambda_1 \lambda_2^2-\lambda_1+\lambda_2\right) \left(1+\lambda_2^2-2 \lambda_2 \cos \left(\frac{2 \pi k}{N}\right)\right)}- \frac{\lambda_1^3 \lambda_2^2 \sigma^2}{\left(\lambda_1^2 \lambda_2-\lambda_1 \lambda_2^2-\lambda_1+\lambda_2\right) \left(1+\lambda_1^2-2 \lambda_1 \cos \left(\frac{2 \pi k}{N}\right)\right)} \end{eqnarray}$

となります (計算間違いや，打ち間違いがあればすいません)．つまり，以下のような形になり，２つの1次自己回帰過程のパワースペクトルの足し算になってます (下図参照)．

$\begin{eqnarray} S_y(f_k) &=& \frac{c_1}{1+\lambda_2^2-2 \lambda_2 \cos \left(\frac{2 \pi k}{N}\right)} + \frac{c_2}{1+\lambda_1^2-2 \lambda_1 \cos \left(\frac{2 \pi k}{N}\right)} \end{eqnarray}$

(そのうち，一般的な分解法について詳しく説明します)．

パワースペクトルの例． $a_1 = 0.9$ ， $a_2 = 0.7$ ．破線が１次自己回帰過程のパワースペクトル．青と緑の破線を足すと，赤実線に一致します．

パワースペクトルの例． $a_1 = 0.9$ ， $a_2 = 0.7$ ．破線が１次自己回帰過程のパワースペクトル．青と緑の破線を足すと，赤実線に一致します．

パワースペクトルが振動のピークを持つ場合

上の条件と違うとき，つまり，

　 $a_1^2+4 a_2 < 0$ のとき，振動しちゃいます．

このとき， $4 a_2 < - a_1^2 \le 0$ なので， $a_2$ はマイナスの値です．

ここでは，最初に導いた２次自己回帰過程のパワースペクトル，

　 $\begin{eqnarray} S_y(f_k) &=& \frac{\sigma^2}{1 + a_1^2 + a_2^2 - 2 a_1 (1 - a_2) \cos(2 \pi f_k) - 2 a_2 \cos( 4 \pi f_k)} \end{eqnarray}$

を使います． $0 < f < 1/2$ の領域で，パワースペクトルが振動に対応するピークを持つということは， $S_y(f_k)$ の分母が下に凸なグラフになり，最小値を一つ持ちます (小さい値で割ると，値が大きいということだぞ)．ということで， $a_2$ がマイナスであることに注意して，分母のcosを含む部分の最小値を評価すると ( $\cos( 2 \pi f)$ の２次関数なので平方完成)，
　 $\begin{eqnarray} - 2 a_1 (1 - a_2) \cos(2 \pi f_k) - 2 a_2 \cos( 4 \pi f_k) &=& - 2 a_2 \left\{2 \cos^2( 2 \pi f_k) - 1\right\}- 2 a_1 (1 - a_2) \cos(2 \pi f_k) \\ &=& - 4 a_2 \cos^2( 2 \pi f_k) - 2 a_1 (1 - a_2) \cos(2 \pi f_k) + 2 a_2 \\ &=& - 4 a_2 \left( \cos( 2 \pi f_k) - \frac{a_1 (a_2 - 1)}{4 a_2} \right)^2 + \frac{a_1^2 (a_2-1)^2}{4 a_2} + 2 a_2 \end{eqnarray}$

となり，最小値（頂点）になる $f_k$ は，

　 $\begin{eqnarray} \cos( 2 \pi f_k) &=& \frac{a_1 (a_2 - 1)}{4 a_2} \\ 2 \pi f_k &=& \cos^{-1} \left(\frac{a_1 (a_2 - 1)}{4 a_2}\right) \\ f_k &=& \frac{1}{2 \pi}\cos^{-1} \left(\frac{a_1 (a_2 - 1)}{4 a_2}\right) \\ \end{eqnarray}$

ということで，ピーク周波数は，

　 $\begin{eqnarray} f_k &=& \frac{1}{2 \pi}\cos^{-1} \left(\frac{a_1 (a_2 - 1)}{4 a_2}\right) \end{eqnarray}$

です．英語のWikipediaには，結果のみ示してあります．
en.wikipedia.org

学生が，記事の更新が早すぎると不平を言いますが，理解してほしい内容にたどり着くまでには，このペースでも，１ヶ月くらいはかかるかもしれません．

おまけ

グラフを書くのに使用したRスクリプト．

ピークがある場合．

# ２次自己回帰仮定のパワースペクトル
f.ar2 <- function(x,a,b){
 return(1/(1+a^2+b^2-2*a*(1-b)*cos(2*pi*x)-2*b*cos(4*pi*x)))
}
f.ar11 <- function(x,a,b){
 return(a^2*b^3/(a^2*b-a*b^2-a+b)/(1+b^2-2*b*cos(2*pi*x)))
}
f.ar12 <- function(x,a,b){
 return(-a^3*b^2/(a^2*b-a*b^2-a+b)/(1+a^2-2*a*cos(2*pi*x)))
}
a1 <- 0.9*sqrt(3)
a2 <- -0.81
lambda1 <- (-a1+sqrt(a1^2+4*a2))/(2*a2)
lambda2 <-  (-a1-sqrt(a1^2+4*a2))/(2*a2)
# グラフを描いてみる
par(mfrow=c(1,2),cex.axis=1.2,cex.lab=1.5)
curve(f.ar2(x,a1,a2),xlim=c(0,1/2),n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,las=1,xaxs="i")
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
# 両対数目盛
curve(f.ar2(x,a1,a2),xlim=c(0.01,1/2),log="xy",n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,xaxt="n",yaxt="n",xaxs="i")
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
####################

１次自己回帰過程の重ね合わせになる場合．

# ２次自己回帰仮定のパワースペクトル
f.ar2 <- function(x,a,b){
 return(1/(1+a^2+b^2-2*a*(1-b)*cos(2*pi*x)-2*b*cos(4*pi*x)))
}
f.ar11 <- function(x,a,b){
 return(a^2*b^3/(a^2*b-a*b^2-a+b)/(1+b^2-2*b*cos(2*pi*x)))
}
f.ar12 <- function(x,a,b){
 return(-a^3*b^2/(a^2*b-a*b^2-a+b)/(1+a^2-2*a*cos(2*pi*x)))
}
a1 <- 0.9
a2 <- 0.7
lambda1 <- (-a1+sqrt(a1^2+4*a2))/(2*a2)
lambda2 <-  (-a1-sqrt(a1^2+4*a2))/(2*a2)
# グラフを描いてみる
par(mfrow=c(1,2),cex.axis=1.2,cex.lab=1.5)
curve(f.ar2(x,a1,a2),xlim=c(0,1/2),n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,las=1,xaxs="i",ylim=c(0,3),yaxs="i")
curve(f.ar11(x,lambda1,lambda2),xlim=c(0,1/2),add=TRUE,n=1000,col=4,lwd=2,lty=2)
curve(f.ar12(x,lambda1,lambda2),xlim=c(0,1/2),add=TRUE,n=1000,col=3,lwd=2,lty=2)
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
# 両対数目盛
curve(f.ar2(x,a1,a2),xlim=c(0.01,1/2),log="xy",n=1000,ylab="S(f)",xlab="f",col=2,lwd=2,xaxt="n",yaxt="n",xaxs="i",ylim=c(0.01,3))
if(a1^2+4*a2 < 0) abline(v=acos(a1*(a2-1)/(4*a2))/(2*pi),col=8,lty=2)
curve(f.ar11(x,lambda1,lambda2),xlim=c(0.01,1/2),add=TRUE,n=1000,col=4,lwd=2,lty=2)
curve(f.ar12(x,lambda1,lambda2),xlim=c(0.01,1/2),add=TRUE,n=1000,col=3,lwd=2,lty=2)
# 対数目盛を描く
axis(1,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(1,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:5,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:5),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
####################