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

本日のセミナーの課題を，私なりこなしました．パワースペクトルを推定するときは，時系列をまるごと1本としてパワースペクトルを計算するよりも，時系列を部分区間に分割して，複数のパワースペクトルの平均をとった方が，パワースペクトル推定のばらつきが減ります．研究室の学生以外，意味がよく分からないと思います．すいませんが，読み飛ばしてください．

平均をとったパワースペクトルの比較

　長さ $N$ の時系列全体からパワースペクトルを計算した結果と，長さ $L$ の部分区間に分割した時系列のパワースペクトルについて平均をとったものを比較します．

　以下のRスクリプトで，Lの値をいろいろと変えてみてください．

# データ長
N <- 3200
# 部分時系列の長さ
L <- 200
##################
# 分割数
n.sub <- floor(N/L)
##################
# 白色ノイズを生成
y <- rnorm(N)
####################
# 全体のパワースペクトル (分割なし)
m <- N-1
##################
# 周波数の定義
freq.N <- (0:m)/(2*m)
##################
a.cov <- acf(y,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
# fftでフーリエ・コサイン変換 [0:1/2]のみ代入
PSD.N <- Re(fft(c(a.cov,rev(a.cov[-1]))))[1:(m+1)]
####################
# 平均化
m <- L-1
for(i in 1:n.sub){
 y.sub <- y[(1+(i-1)*L):(i*L)]
 ######################
 # 自己共分散の推定 (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)
#################
# パワースペクトル（横軸のみ対数スケール）
par(mfrow=c(1,2),cex.axis=1.3,cex.lab=1.5)
########
plot(freq.N[-1],PSD.N[-1],"l",col=4,log="y",main=paste("Periodogram (N = ",N,")",sep=""),xlab="f",yaxt="n",xaxs="i",ylab="",las=1,xlim=c(0,0.5),ylim=c(0.001,10))
abline(h=0,lty=2)
# 対数目盛を描く
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:2,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:2),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########
plot(freq[-1],PSD.ave[-1],"l",col=4,log="y",main=paste("Averaged Periodogram (L = ",L,")",sep=""),xlab="f",yaxt="n",xaxs="i",ylab="",las=1,xlim=c(0,0.5),ylim=c(0.001,10))
abline(h=0,lty=2)
# 対数目盛を描く
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:2,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:2),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))

平均をとったピリオドグラムの分散

　省略します．

ラグを小さくとったピリオドグラムの比較と，分割平均したピリオドグラムの比較

　長さ $N$ の時系列1本に対して最大ラグ $L-1$ としてピリオドグラムを計算した場合と，長さ $L$ の部分時系列に分割して平均ピリオドグラムを計算した場合の比較です．

ピリオドグラムの計算結果．(左) 1本のサンプル時系列 (N=3200, L=200)の結果．(右) 100本のサンプル時系列の平均と標準偏差

　一致しません．

# データ長
N <- 3200
# 部分時系列の長さ
L <- 200
##################
# 分割数
n.sub <- floor(N/L)
##################
# 最大ラグ
m <- L-1
####################
# 周波数の定義
freq <- (0:m)/(2*m)
####################
N.samps <- 100
for(k in 1:N.samps){
 ##################
 # 白色ノイズを生成
 y <- rnorm(N)
 ##################
 a.cov <- acf(y,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
 # fftでフーリエ・コサイン変換 [0:1/2]のみ代入
 PSD.N <- Re(fft(c(a.cov,rev(a.cov[-1]))))[1:(m+1)]
 ####################
 # 平均化
 for(i in 1:n.sub){
  y.sub <- y[(1+(i-1)*L):(i*L)]
  ######################
  # 自己共分散の推定 (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))
  }
 }
 ##################
 PSD.ave <- rowMeans(PSD.samps)
 if(k == 1){
  PSD.L <- data.frame(PSD.N)
  PSD.a <- data.frame(PSD.ave)

  #################
  # パワースペクトル（横軸のみ対数スケール）
  par(mfrow=c(1,2),cex.axis=1.3,cex.lab=1.5)
  ########
  plot(freq,PSD.N,"l",col=4,log="y",xlab="f",yaxt="n",xaxs="i",ylab="",las=1,xlim=c(0,0.5),ylim=c(0.1,10),lwd=2)
  abline(h=0,lty=2)
  # 対数目盛を描く
  axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
  tmp<-paste(paste(sep="","expression(10^",-5:2,")"),collapse=",")
  v.label<-paste("axis(2,las=1,at=10^(-5:2),label=c(",tmp,"),tck=-0.03)",sep="")
  eval(parse(text=v.label))
  ########
  lines(freq,PSD.ave,col=2,lwd=2)
  legend("topright",legend=c("Single L-1 autocov","Averaged"),col=c(4,2),lwd=2)
 }else{
  PSD.L <- cbind(PSD.L,data.frame(PSD.N))
  PSD.a <- cbind(PSD.a,data.frame(PSD.ave))
 }
}
PSD.L.mean <- rowMeans(PSD.L)
PSD.a.mean <- rowMeans(PSD.a)
PSD.L.sd <- apply(PSD.L,1,sd)
PSD.a.sd <- apply(PSD.a,1,sd)
plot(freq,PSD.L.mean,"l",col=4,log="y",xlab="f",yaxt="n",xaxs="i",ylab="",las=1,xlim=c(0,0.5),ylim=c(0.1,10),lwd=2)
lines(freq,PSD.L.mean+PSD.L.sd,col=rgb(0, 0, 1, alpha=0.2),lwd=2)
lines(freq,PSD.L.mean-PSD.L.sd,col=rgb(0, 0, 1, alpha=0.2),lwd=2)
abline(h=0,lty=2)
legend("topright",legend=c("Single L-1 autocov","Averaged"),col=c(4,2),lwd=2)
# 対数目盛を描く
axis(2,at=10^(-5:5)%x%(1:9),label=FALSE,tck=-0.02)
tmp<-paste(paste(sep="","expression(10^",-5:2,")"),collapse=",")
v.label<-paste("axis(2,las=1,at=10^(-5:2),label=c(",tmp,"),tck=-0.03)",sep="")
eval(parse(text=v.label))
########
lines(freq,PSD.a.mean,col=2,lwd=2)
lines(freq,PSD.a.mean+PSD.a.sd,col=rgb(1, 0, 0, alpha=0.2),lwd=2)
lines(freq,PSD.a.mean-PSD.a.sd,col=rgb(1, 0, 0, alpha=0.2),lwd=2)