【Ｒで時系列解析】1次自己回帰過程の自己共分散関数とパワースペクトルの推定

今日のセミナーの課題は，
問題　数値的に生成した1次自己回帰過程のサンプル時系列の，自己共分散関数とパワースペクトル（ピリオドグラム）を推定し，理論値と比較せよ．推定は，複数のサンプル時系列を使って平均として求めよ．

解答としては，以下のような図を想定しています．ただし，この図では，ラグと周波数の領域を欲張りすぎで，実用上は，結果として参照するラグの上限はN/10（最大でもN/4）程度，周波数の上限は10/N（最小でも4/N）程度にしてください．

実行例（時系列の長さN=2^13=8192，サンプル数N.samps=2^6=64）

上の図のような数値実験の結果を作るＲスクリプトは，以下です（今回は，フーリエ・コサイン変換にfftを使いました）．

# 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.99
# w[n]の標準偏差sig^2
sig2 <- 1
# 時系列の長さ
N <- 2^13
####
y <- c()
####################
# 平均をとるサンプル時系列数
N.samps <- 2^6
####################
# 最大ラグの決定
m <- N-1
####################
# 周波数の定義
freq <- (0:m)/(2*m)
####################
PSD0 <- c()
####################
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]
 }
 ######################
 # 自己共分散の推定 (sample autocovariance)
 a.cov <- acf(y,lag.max=m,type="covariance",plot=FALSE)$acf[,,1]
# fftでフーリエ・コサイン変換
 PSD0 <- Re(fft(c(a.cov,rev(a.cov[-1]))))
 ####################
 # 平滑化（これがないと，スペクトルがマイナスになることがあります）
 PSD <- numeric(m+1)
 PSD[2:m] <- PSD0[2:m]*0.5+0.25*PSD0[1:(m-1)]+0.25*PSD0[3:(m+1)]
 PSD[1] <- (PSD0[1]+PSD0[2])/2
 PSD[m+1] <- (PSD0[m]+PSD0[m+1])/2
 ######################
 if(ii == 1){
  PSD.samps <- data.frame(PSD)
  AutoCov.samps <- data.frame(a.cov)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
  AutoCov.samps <- cbind(AutoCov.samps,data.frame(a.cov))
 }
}
f <- freq #[freq>=0 & freq<=0.5]
Spec <- rowMeans(PSD.samps)
AutoCov <- rowMeans(AutoCov.samps)
#################
# 図を描く
par(mfrow=c(1,3))
#################
# パワースペクトル（横軸のみ対数スケール）
plot(f[-1],Spec[-1],"l",col=4,log="x",main="Power spectrum density",xlab="f",xaxt="n",xaxs="i",ylab="",las=1)
abline(h=0,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))
########
# 解析的に計算したパワースペクトル
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))
#######
# パワースペクトル（両対数スケール）
plot(f[-1],Spec[-1],"l",col=4,log="xy",main="Power spectrum density",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))
########
# 解析的に計算したパワースペクトル
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))
###########################################
# 自己共分散関数（横軸のみ対数スケール）
plot(1:m,AutoCov[-1],"l",log="x",col=4,main="Autocovariance",xlab="f",xaxt="n",xaxs="i",ylab="",las=1)
abline(h=0,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))
########
# 解析的に計算した自己共分散関数
curve(sig2/(1-a^2)*a^x,add=TRUE,col=2,lwd=2,lty=2,n=1000)

実行するごとに結果は異なり，自己共分散の短いラグ側でずれが大きくなることがある．