【Rで時系列解析】パワースペクトルの平均化処理 - ケィオスの時系列解析メモランダム

　パワースペクトルを推定する数値実験を，少しずつやっていきたいと思います．

　平均化処理ってやったほうがいいかも，というのが今回の感想です．

　Rでパワースペクトルを推定するとき，"spectrum"というコマンドを使えばいいんでしょ，と考えている人も多いと思います．Rのコマンド"spectrum"では，デフォルトの前処理として自動で直線トレンドが除かれ，窓関数がかけられます．さらに，"span"のオプションを指定すれば平滑化をしてくれます．私は正直言って，何も考えずに，Rのコマンド"spectrum"を使ってきました．最近，それじゃあ，研究者としてダメなんじゃないかと思いはじめたので，自分で納得できるパワースペクトルの推定法を追求してみたいと思います．

　ということで，今回は１本の時系列を部分時系列に分割し，部分時系列のパワースペクトルを周波数ごとに平均することで，パワースペクトルを推定してみます．トレンド除去とか，窓関数をかけるとか，平滑化するとか，そういうテクニックもありますが，今回はそれは使いません．単に平均化処理をするだけです．トレンドあり，トレンドなしを比較してみます．

検討する確率過程

　今回は，以下の確率過程のサンプル時系列を生成し，パワースペクトルを推定します．

1．１次自己回帰過程
　差分方程式は　

　　 $y[n] = a_1 y[n-1] + w[n]$

です．ここで， $a_1$ は実数のパラメタ， $\{w[n]\}$ は平均0，分散 $\sigma^2$ の白色ノイズです．このとき，自己共分散 (自己相関)関数 $C(k)$ は

$\begin{eqnarray} C(k) &=& \frac{ \sigma^2 }{1-a_1^2} a_1^{k} \end{eqnarray}$

です．

2．２次自己回帰過程
　差分方程式は　

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

です．ここで， $a_1$ ， $a_2$ は実数のパラメタ， $\{w[n]\}$ は平均0，分散 $\sigma^2$ の白色ノイズです．このとき，自己共分散 (自己相関)関数 $C(k)$ は

$\begin{eqnarray} C(0) &=& \frac{\sigma^2 (1-a_2)}{(1+a_2)(1-a_1-a_2)(1+a_1-a_2)} \\ C(1) &=& \frac{a_1}{1-a_2} C(0) \\ C(k) &=& a_1 C(k-1) + a_2 C(k-2) \qquad (C \ge 2) \end{eqnarray}$

で計算できます．

3．非整数ガウスノイズ (fractional Gaussian noise)
　自己共分散 (自己相関)関数 $C(k)$ が

$\displaystyle C(k) = \frac{\sigma^2}{2} \left( \left|k+1 \right|^{2H} - 2 \left|k \right|^{2H} + \left|k - 1 \right|^{2H} \right)$

で与えられる過程です．ここで， $\sigma^2$ は，分散です．

実験その１：トレンドがない理想的なサンプル時系列の分析

　上の各過程について長さNのサンプル時系列を生成します．時系列の生成方法はこの記事を参考にしてください．追加のパッケージは，必要ありません．

chaos-kiyono.hatenablog.com

　パワースペクトルの推定では，「１本まるごと使ってパワースペクトル (ピリオドグラム)を計算するパターン」と，「時系列を分割して求めたパワースペクトルの平均処理をするパターン」を比較します．

　パワースペクトル (ピリオドグラム)は，まず，自己共分散 (自己相関)関数を推定し，それをフーリエコサイン変換して求めます．時系列を分割して求めたパワースペクトルは，周波数ごとに平均値を求め，その結果をパワースペクトルの推定とします．

　ここでたくさん図を載せてもしょうがないので，各自で数値実験してください．Rスクリプトは，下の方にまとめて掲載しておきます．Rスクリプトでは，Nが時系列の長さ，n.subが分割数になってます．以下が，Rスクリプトの実行例です．

　図は，左がサンプル時系列，中が１本まるごとのパワースペクトルの推定結果，右が１６分割で平均処理したパワースペクトルの推定結果です．赤い破線が，解析的な結果で，この線に一致しているほど良いということです．

トレンドなしの１次自己回帰過程 ( $a_1=0.9$ )

トレンドなしの2次自己回帰過程 ( $a_1=0.9 \sqrt{3}$ ， $a_2=-0.81$ )

トレンドなしの非整数ガウスノイズ ( $H=0.9$ )

トレンドなしの非整数ガウスノイズ ( $H=0.9$ )．N= $2^{14}$ ，n.sub=16．

実験その2：直線トレンドがあるサンプル時系列の分析

　直線のトレンドを時系列に足しました (ちょっと極端な例です)．Rスクリプトは掲載してません．下のRスクリプトを自分で改造して，実験してみてください．図は，左がサンプル時系列，中が１本まるごとのパワースペクトルの推定結果，右が１６分割で平均処理したパワースペクトルの推定結果です．赤い破線が，解析的な結果で，この線に一致しているほど良いということです．

　ここでは，図の真ん中にある１本まるごとのパワースペクトルの推定結果に変な構造が現れる点に注目してください．これは，直線トレンドの影響です．図の右側のように分割するとトレンドの影響がずいぶん軽減されます．

直線トレンドありの１次自己回帰過程 ( $a_1=0.9$ )

直線トレンドありの2次自己回帰過程 ( $a_1=0.9 \sqrt{3}$ ， $a_2=-0.81$ )

直線トレンドありの非整数ガウスノイズ ( $H=0.9$ )

直線トレンドありの非整数ガウスノイズ ( $H=0.9$ )．N= $2^{14}$ ，n.sub=16．左がサンプル時系列．

まとめ

　直線トレンドありのとき，上図中央の両対数プロットに直線的変化が現れてます．これは，直線トレンドが生み出す構造です．時系列を分割するだけで，直線トレンドの影響がずいぶん軽減できていることが分かると思います．

Rスクリプト

1．トレンドなしの1次自己回帰過程の数値実験

# 部分区間の数
n.sub <- 16
# 時系列の長さ (時系列は奇数長になります)
N <- 2^14
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# モデルの自己共分散関数を与える
# 例：AR(1)過程の自己共分散関数
AR1.model <- function(t,a,sig2){
 return(sig2*a^t/(1-a^2))
}
# パラメタの設定
a <- 0.9
sig2 <- 1
# 【注意】[-M,M]区間ではなく，[0,2M-1]にしている
acov.model <- c(AR1.model(0:M,a,sig2),AR1.model(M:1,a,sig2))
lag <- c(0:M,-(M:1))
#########################
# 自己共分散関数のフーリエ変換
fft.model <- fft(acov.model)
PSD.model <- Re(fft.model)
f <- c((0:M)/N,-(M:1)/(N))
#########################
# 白色ノイズの生成
WN <- rnorm(N)
#　ホワイトノイズのフーリエ変換
fft.WN <- fft(WN)
#########################
# サンプル時系列の生成
fft.sim <-  sqrt(PSD.model)*fft.WN
x.sim <- Re(fft(fft.sim,inverse=TRUE))/N
#########################
# 結果の描画
par(mfrow=c(1,3),cex.main=1.5,cex.axis=1.4,cex.lab=1.5,las=1)
plot(0:(N-1),x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main="Synthetic Time Series",lwd=2)
#　自己共分散の推定
est.acov.sim <- acf(x.sim,plot=FALSE,type="covariance",lag.max=M)$acf[,,1]
# PSDの推定
est.fft.sim <- fft(x.sim)
est.psd.sim <- abs(est.fft.sim)^2/N
est.psd.sim <- est.psd.sim[f >= 0]
#
plot(f[f > 0],est.psd.sim[f > 0],"l",log="xy",col=4,pch=16,xlim=c(f[2],0.5),xaxs="i",xlab="f",ylab="Periodogram",main=paste("Number of subintervals =",1,"(No averaging)"),lwd=2,xaxt="n",yaxt="n")
lines(f[f > 0],PSD.model[f > 0],col=2,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))
########
legend("bottomleft", legend=c("Analytical Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))
##########################
# 部分区間
i.sub <- seq(1,N,length.out=n.sub+1)
m <- i.sub[2]-2
##########################
for(i in 1:n.sub){
 y.sub <- x.sim[i.sub[i]:(i.sub[i+1]-1)]
 ######################
 # 自己共分散の推定 (sample autocovariance)
 a.cov <- acf(y.sub,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
 # fftでフーリエ・コサイン変換 [0:1/2]のみ代入
 PSD <- Re(fft(c(a.cov,rev(a.cov[-1]))))[1:(m+1)]
 ######################
 if(i == 1){
  PSD.samps <- data.frame(PSD)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
 }
}
##################
# 周波数の定義
freq <- (0:m)/(2*m)
PSD.ave <- rowMeans(PSD.samps)
#################
plot(freq[freq > 0],PSD.ave[freq > 0],"l",log="xy",col=4,pch=16,xlim=c(freq[2],0.5),xaxs="i",xlab="f",ylab="Periodogram",main=paste("Number of subintervals =",n.sub),lwd=2,xaxt="n",yaxt="n")
lines(f[f > 0],PSD.model[f > 0],col=2,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))
########
legend("bottomleft", legend=c("Analytical Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))

2．トレンドなしの2次自己回帰過程の数値実験

# 部分区間の数
n.sub <- 16
# 時系列の長さ (時系列は奇数長になります)
N <- 2^14
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# モデルの自己共分散関数を与える
# 例：AR(2)過程の自己共分散関数
AR2.model <- function(n,a1,a2,sig2){
 Cov <- c()
 Cov[1] <- sig2*(1-a2)/(1-a1^2-a2-a1^2*a2-a2^2+a2^3)
 Cov[2] <- Cov[1]*a1/(1-a2)
 for(k in 3:(n+1)){
  Cov[k] <- a1*Cov[k-1] + a2*Cov[k-2]
 }
 return(Cov)
}
# パラメタの設定
a1 <- 0.9*sqrt(3)
a2 <- -0.81
sig2 <- 1
# 【注意】[-M,M]区間ではなく，[0,2M-1]にしている
acov.model <- c(AR2.model(M,a1,a2,sig2),rev(AR2.model(M,a1,a2,sig2)[-1]))
lag <- c(0:M,-(M:1))
#########################
# 自己共分散関数のフーリエ変換
fft.model <- fft(acov.model)
PSD.model <- Re(fft.model)
f <- c((0:M)/N,-(M:1)/(N))
#########################
# 白色ノイズの生成
WN <- rnorm(N)
#　ホワイトノイズのフーリエ変換
fft.WN <- fft(WN)
#########################
# サンプル時系列の生成
fft.sim <-  sqrt(PSD.model)*fft.WN
x.sim <- Re(fft(fft.sim,inverse=TRUE))/N
#########################
# 結果の描画
par(mfrow=c(1,3),cex.main=1.5,cex.axis=1.4,cex.lab=1.7,las=1,mar=c(5,5,2,2))
plot(0:(N-1),x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main="Synthetic Time Series",lwd=2)
#　自己共分散の推定
est.acov.sim <- acf(x.sim,plot=FALSE,type="covariance",lag.max=M)$acf[,,1]
# PSDの推定
est.fft.sim <- fft(x.sim)
est.psd.sim <- abs(est.fft.sim)^2/N
est.psd.sim <- est.psd.sim[f >= 0]
#
plot(f[f > 0],est.psd.sim[f > 0],"l",log="xy",col=4,pch=16,xlim=c(f[2],0.5),xaxs="i",xlab="f",ylab="Periodogram",main=paste("Number of subintervals =",1,"(No averaging)"),lwd=2,xaxt="n",yaxt="n")
lines(f[f > 0],PSD.model[f > 0],col=2,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))
########
legend("bottomleft", legend=c("Analytical Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))
##########################
# 部分区間
i.sub <- seq(1,N,length.out=n.sub+1)
m <- i.sub[2]-2
##########################
for(i in 1:n.sub){
 y.sub <- x.sim[i.sub[i]:(i.sub[i+1]-1)]
 ######################
 # 自己共分散の推定 (sample autocovariance)
 a.cov <- acf(y.sub,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
 # fftでフーリエ・コサイン変換 [0:1/2]のみ代入
 PSD <- Re(fft(c(a.cov,rev(a.cov[-1]))))[1:(m+1)]
 ######################
 if(i == 1){
  PSD.samps <- data.frame(PSD)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
 }
}
##################
# 周波数の定義
freq <- (0:m)/(2*m)
PSD.ave <- rowMeans(PSD.samps)
#################
plot(freq[freq > 0],PSD.ave[freq > 0],"l",log="xy",col=4,pch=16,xlim=c(freq[2],0.5),xaxs="i",xlab="f",ylab="Periodogram",main=paste("Number of subintervals =",n.sub),lwd=2,xaxt="n",yaxt="n")
lines(f[f > 0],PSD.model[f > 0],col=2,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))
########
legend("bottomleft", legend=c("Analytical Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))

3．トレンドなしの非整数ガウスノイズの数値実験

# 部分区間の数
n.sub <- 16
# 時系列の長さ (時系列は奇数長になります)
N <- 2^14
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# モデルの自己共分散関数を与える
# 非整数ガウスノイズ過程の自己共分散関数
fGn.model <- function(k,H,sig2){
 return(sig2*(abs(k+1)^(2*H)-2*abs(k)^(2*H)+abs(k-1)^(2*H))/2)
}
# パラメタの設定
H <- 0.9
sig2 <- 3
# 【注意】[-M,M]区間ではなく，[0,2M-1]にしている
acov.model <- c(fGn.model(0:M,H,sig2),fGn.model(M:1,H,sig2))
lag <- c(0:M,-(M:1))
#########################
# 自己共分散関数のフーリエ変換
fft.model <- fft(acov.model)
PSD.model <- Re(fft.model)
f <- c((0:M)/N,-(M:1)/(N))
#########################
# 白色ノイズの生成
WN <- rnorm(N)
#　ホワイトノイズのフーリエ変換
fft.WN <- fft(WN)
#########################
# サンプル時系列の生成
fft.sim <-  sqrt(PSD.model)*fft.WN
x.sim <- Re(fft(fft.sim,inverse=TRUE))/N
#########################
# 結果の描画
par(mfrow=c(1,3),cex.main=1.5,cex.axis=1.4,cex.lab=1.5,las=1)
plot(0:(N-1),x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main="Synthetic Time Series",lwd=2)
#　自己共分散の推定
est.acov.sim <- acf(x.sim,plot=FALSE,type="covariance",lag.max=M)$acf[,,1]
# PSDの推定
est.fft.sim <- fft(x.sim)
est.psd.sim <- abs(est.fft.sim)^2/N
est.psd.sim <- est.psd.sim[f >= 0]
#
plot(f[f > 0],est.psd.sim[f > 0],"l",log="xy",col=4,pch=16,xlim=c(f[2],0.5),xaxs="i",xlab="f",ylab="Periodogram",main=paste("Number of subintervals =",1,"(No averaging)"),lwd=2,xaxt="n",yaxt="n")
lines(f[f > 0],PSD.model[f > 0],col=2,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))
########
legend("bottomleft", legend=c("Analytical Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))
##########################
# 部分区間
i.sub <- seq(1,N,length.out=n.sub+1)
m <- i.sub[2]-2
##########################
for(i in 1:n.sub){
 y.sub <- x.sim[i.sub[i]:(i.sub[i+1]-1)]
 ######################
 # 自己共分散の推定 (sample autocovariance)
 a.cov <- acf(y.sub,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
 # fftでフーリエ・コサイン変換 [0:1/2]のみ代入
 PSD <- Re(fft(c(a.cov,rev(a.cov[-1]))))[1:(m+1)]
 ######################
 if(i == 1){
  PSD.samps <- data.frame(PSD)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
 }
}
##################
# 周波数の定義
freq <- (0:m)/(2*m)
PSD.ave <- rowMeans(PSD.samps)
#################
plot(freq[freq > 0],PSD.ave[freq > 0],"l",log="xy",col=4,pch=16,xlim=c(freq[2],0.5),xaxs="i",xlab="f",ylab="Periodogram",main=paste("Number of subintervals =",n.sub),lwd=2,xaxt="n",yaxt="n")
lines(f[f > 0],PSD.model[f > 0],col=2,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))
########
legend("bottomleft", legend=c("Analytical Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))