【Rで時系列解析】弱定常過程のスペクトル表現定理 - ケィオスの時系列解析メモランダム

弱定常な確率過程のスペクトル表現定理が意味することを理解するために数値実験をしてみました．今回は，前回の記事
chaos-kiyono.hatenablog.com
の続きです．

　弱定常過程のパワースペクトルについて理解してほしい部分は，「パワースペクトルは，確率過程を周波数成分に分解したときの振幅を，どのような意味で指定しているのか」ということです．パワースペクトルは，周波数成分の振幅の情報をもっていますが，これは，振幅の２乗の期待値の意味であり，振幅を決定論的に決めているわけではありません．

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

　例として１次自己回帰過程を考えます．

　 $y_{n}=a y_{n-1} + w_n$

ここで， $a$ は， $-1 < a < 1$ の定数， $w_n$ は，平均0，分散 $\sigma^2$ の白色ノイズです．このとき，パワースペクトル $S(f)$ は，

　 $\displaystyle S(f)=\frac{\sigma^2}{1+a^2-2 a \cos (2 \pi f)}=\frac{\sigma^2}{(1-a)^2-2 a( \cos (2 \pi f)-1)}$

です．

　自己相関至上主義 (英: autocorrelationism)では，パワースペクトル $S(f)$ (と平均値)が分かれば，この過程のことはすべて理解できるということになります．では，このパワースペクトルは，弱定常な確率過程の何を指定しているのでしょうか？このことを前回の記事で説明したつもりですが，伝わった気がしないし，説明が不足している部分があると思いますので，数値実験して，説明を補足します．

１次自己回帰過程を周波数成分に分解

　前回の記事で説明したスペクトル表現定理からわかることは，平均0の弱定常な離散確率過程の時系列

　 $\left\{x[0], x[1], x[2], \cdots, x\left[\left(N-1\right) \right] \right\}$

を，

　 $\displaystyle x[n] = \sum_{j = 1}^{M} \left[A_{j} \cos \left(2 \pi f_j n \right) + B_{j} \sin \left(2 \pi f_j n \right) \right]$

と分解したとき ( $f_j = j/N$ ， $M = \lfloor N/2 \rfloor$ )，

　 $\begin{eqnarray} {\rm E} \left(A_{j} \right) &=& 0\\ {\rm E} \left(B_{j} \right) &=& 0\end{eqnarray}$

です．さらに， $S \left( f_j \right)$ をこの過程の片側パワースペクトルとすると，

　 $\begin{eqnarray} {\rm E} \left(A_{j} \, A_{k} \right) &=& \left\{\begin{array}{lc} \displaystyle \frac{S \left( f_j \right)}{N} & \left( j = k \right) \\ \displaystyle 0 & \left( j \neq k \right) \end{array}\right. \\ {\rm E} \left(B_{j} \, B_{k} \right) &=& \left\{\begin{array}{lc} \displaystyle \frac{S \left( f_j \right)}{N} & \left( j = k \right) \\ \displaystyle 0 & \left( j \neq k \right) \end{array}\right. \\ {\rm E} \left(A_{j} \, B_{k} \right) &=& 0 \quad \left({\rm all} \ j \ {\rm and} \ k \right)\end{eqnarray}$

です． $S \left( f_j \right)$ は，時系列から推定したパワースペクトルではなく，真のパワースペクトルということに注意してください．

　つまり，弱定常過程を

　 $\displaystyle x[n] = \sum_{j = 1}^{M} \left[A_{j} \cos \left(2 \pi f_j n \right) + B_{j} \sin \left(2 \pi f_j n \right) \right]$

と分解したとき， $A_{j}$ と $B_{j}$ は平均0，分散 $\displaystyle \frac{S \left( f_j \right)}{N}$ の正規乱数ということです．前回， $A_{j}$ と $B_{j}$ が正規分布に従う確率変数とは説明していませんでした．

　このことを数値実験で確かめてみました．以下に掲載したRスクリプトで，N.samps <- 10000とした結果が下の図です．

　上の図の一番左は，真のパワースペクトル (赤破線)と，サンプルの平均として推定したパワースペクトル (青実線)です．各段は，パワースペクトル中に縦線で示した周波数成分の結果を示してます．左から２番目の図は，縦軸に係数 $A_{j}$ ，横軸に係数 $B_{j}$ をプロットしたものです．たくさん点があるのは，10000個のサンプルの結果をプロットしたからです．図中の赤破線の円は，パワースペクトルが表す振幅を描いたものです．右側の２つの図は， $A_{j}$ と $B_{j}$ の分布を描いたものです．赤破線は，理論的に予言される正規分布です．

　上の図の左から２番目の図から分かるように， $A_{j}$ と $B_{j}$ で与えられる振幅は，赤破線の円のまわりにあるのではなく，原点を中心としてばらついているということです．右側の２つの図から，0が中心であり，正規分布にしたがってばらつくということが確認できます．

まとめ

　パワースペクトルは，パソコンを使って簡単に計算できるので，意味を理解せずに，とりあえずパワースペクトルを計算してみる人が多いように感じます．私自身の理解は，まだまだと感じていますので，これからも基礎を大事にして勉強したいと思います．

Rスクリプト

　この記事で使用したRスクリプトは以下です．

# 1次自己回帰過程 (first-order autoregressive process) 
# y[n] = a y[n-1] + w[n]
# w[n]: white noise with zero mean and variance sig^2
####################
# パラメタ (parameter)
a <- 0.9
# w[n]の標準偏差sig^2
sig2 <- 1
# 時系列の長さ
N <- 256
####
y <- c()
####################
# 平均をとるサンプル時系列数
N.samps <- 100
####################
# 最大ラグの決定
m <- N-1
####################
# 周波数の定義
freq <- (0:(N-1))/N
####################
for(ii in 1:N.samps){
 ####################
 # 白色ノイズの生成 
 w <- rnorm(N,0,sqrt(sig2))
 ####################
 # 空の変数を事前に準備
 # 初期値をy[1]=w[1]とする
 y[1] <- w[1]
 for(n in 2:N){
  y[n] <- a*y[n-1]+w[n]
 }
 ######################
 PSD <- abs(fft(y))^2/N
 Re.Y <- Re(fft(y))/sqrt(N)
 Im.Y <- Im(fft(y))/sqrt(N)
 ######################
 if(ii == 1){
  PSD.samps <- data.frame(PSD)
  Re.Y.samps <- Re.Y
  Im.Y.samps <- Im.Y
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
  Re.Y.samps <- cbind(Re.Y.samps,Re.Y)
  Im.Y.samps <- cbind(Im.Y.samps,Im.Y)
 }
}
f <- freq[freq>=0 & freq<=0.5]
Spec <- rowMeans(PSD.samps)[freq>=0 & freq<=0.5]
#################
N.plot <- 4
n.f <- length(f)
f.samp <- round(exp(seq(log(2),log(n.f-1),length=N.plot)))
#################
# 図を描く
par(mfrow=c(N.plot,4),cex.lab=1.5,las=1,cex.axis=1.3)
###
for(j in 1:N.plot){
par(pty="m")
#################
plot(f[-1],Spec[-1],"l",col=4,log="xy",main="Power spectrum density",xlim=c(f[2]/1.1,tail(f,1)*1.1),lwd=4,xlab="f",xaxt="n",yaxt="n",xaxs="i",ylab="")
# 対数目盛を描く
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))
abline(v=f[f.samp[j]],lty=2,col="lightpink",lwd=2)
########
# 解析的に計算したパワースペクトル
curve(sig2/(1-2*a*cos(2*pi*x)+a^2),add=TRUE,col=2,lwd=2,lty=2,n=1000,xlim=c(1/N,1/2))
#
points(f[f.samp[j]],sig2/(1-2*a*cos(2*pi*f[f.samp[j]])+a^2),pch=16,cex=2,col=2)
#####
par(pty="s")
r <- sqrt(sig2/(1-2*a*cos(2*pi*f[f.samp[j]])+a^2))/sqrt(2)
#colm <- densCols(Re.Y.samps[f.samp[j],],Im.Y.samps[f.samp[j],], colramp = colorRampPalette(cm.colors(16)[8:1]))
plot(Re.Y.samps[f.samp[j],],Im.Y.samps[f.samp[j],],pch=16,col=rgb(43/255,137/255,206/255,alpha=0.2),xlim=c(-r*3.5,r*3.5),ylim=c(-r*3.5,r*3.5),xlab="Re(Y) / sqrt(N)",ylab="Im(Y) / sqrt(N)",main=paste(f[f.samp[j]],"Hz"))
curve(sqrt(r^2-x^2),xlim=c(-r,r),add=TRUE,col=2,lty=2,lwd=2)
curve(-sqrt(r^2-x^2),xlim=c(-r,r),add=TRUE,col=2,lty=2,lwd=2)
par(pty="m")
r <- sqrt(sig2/(1-2*a*cos(2*pi*f[f.samp[j]])+a^2))/sqrt(2)
hist(Re.Y.samps[f.samp[j],],xlab="Re(Y) / sqrt(N)",breaks="FD",col=4,border=4,freq=FALSE,ylab="",main=paste(f[f.samp[j]],"Hz"))
curve(dnorm(x,0,r),add=TRUE,col=2,lty=2,lwd=2)
hist(Im.Y.samps[f.samp[j],],xlab="Re(Y) / sqrt(N)",breaks="FD",col=4,border=4,freq=FALSE,ylab="",main=paste(f[f.samp[j]],"Hz"))
curve(dnorm(x,0,r),add=TRUE,col=2,lty=2,lwd=2)
}