【Ｒで時系列解析】自分で指定した自己共分散関数に従うサンプル時系列の生成

今回は，自分で自己共分散関数を指定し，その自己共分散関数に従う確率過程のサンプル時系列を生成する方法の説明です．自己共分散とパワースペクトルの関係を理解した上で，Rのfftコマンドを使いこなせれば簡単です．とはいえ，私は理解力がないので，他の人の論文に書いてある方法をそのまま信じ，その記述のいい加減さ（間違い）に振り回されて，結構迷った経験があります．

ここで学ぶポイントは，「平均0のガウス過程は，自己共分散関数のみで完全に特徴づけられる」ということです．以下では，自己共分散関数を指定するだけで，いくらでもサンプル時系列を生成できます．他の情報は必要ありません．

計算の流れは以下です．

自己共分散関数 $C(k)$ を設定
自己共分散関数 $C(k)$ をフーリエ変換し，パワースペクトル $S(f)$ を求める
白色ノイズを生成し，そのフーリエ変換を計算
生成した白色ノイズのパワーが，2で求めた $S(f)$ に従うように変換してNで割る
逆フーリエ変換して時系列を計算

以下では例として，1次自己回帰過程の自己共分散関数
$\displaystyle C(k)=\frac{\sigma^2}{1-a^2} a^k$
を設定して，この $C(k)$ に従う，時系列を生成してみます．

Rスクリプトはこんな感じです．

# 時系列の長さ (時系列は奇数長になります)
N <- 1001
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# モデルの自己共分散関数を与える
# 例：AR(1)過程の自己共分散関数
AutoCov.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(AutoCov.model(0:M,a,sig2),AutoCov.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(3,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="y",col=4,pch=16,xlim=c(0,0.5),xaxs="i",xlab="f",ylab="PSD",main="Power Spectrum",lwd=2)
lines(f[f >= 0],PSD.model[f >= 0],col=2,lwd=2,lty=2)
legend("topright", legend=c("Target Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))
plot(0:M,est.acov.sim,"l",col=4,pch=16,xlim=c(0,M),xaxs="i",xlab="lag",ylab="Autocov",main="Autocovariance",lwd=2)
lines(lag[lag >= 0],acov.model[lag >= 0],col=2,lwd=2,lty=2)
legend("topright", legend=c("Target Autocovariance"), col=c(2), lty=c(2),lwd=c(2))

このスクリプトを実行すると，以下のような図が表示されます．