【Rで時系列解析】白色ノイズのm階積分 (和分)過程のパワースペクトル (数値実験)

前に， $m$ 階積分 (和分)過程のパワースペクトルの記事を書いたときに導いた結果，

分散 $\sigma^2$ の白色ノイズ $\{x [n]\}_{n=1}^N$ の $m$ 階差分過程の時系列 $\{y [n]\}_{n=1}^N$ のパワースペクトル $S_y(f_k)$ は，

　 $\displaystyle S_y(f_k) = \sigma^2 \left(2- 2 \cos 2 \pi f_k \right)^{m} = \sigma^2 \left( 2 \sin \pi f \right)^{2m}$

$m$ 階積分 (和分)過程の時系列 $\{y [n]\}_{n=1}^N$ のパワースペクトル $S_y(f_k)$ は，

　 $\displaystyle S_y(f_k) = \frac{\sigma^2}{\left(2- 2 \cos 2 \pi f_k \right)^{m}} = \frac{\sigma^2}{\left( 2 \sin \pi f \right)^{2m}}$

について，数値実験で確認していなかったので，やっておきます．前の記事はこれです．
chaos-kiyono.hatenablog.com

数値実験結果

　以下に掲載したRスクリプトを実行すると，サンプル時系列１本と，100本のサンプル時系列から計算した平均パワースペクトルが描かれます．

　Rスクリプトでは，

m <- 1

の部分で， $m$ 階積分 (和分)を設定します．mをマイナスの値にすると，差分になります．

平均をとるサンプル時系列の数は，

N.samps <- 100

で設定しています．

時系列の長さは，

N <- 1001

で設定しています．

　下に，結果の図をいくつか掲載しておきますが，上のパラメタをいろいろと変えて，自分で数値実験してみてください．数値実験結果と解析的結果が一致するので，

　 $m$ 階積分 (和分)過程のパワースペクトルの解析的結果は正しそう

ということです．ただし，差分過程は，低周波数側で解析的結果とのずれがみられます．絶対これだという原因は，私にはわからないので，原因が分かる人は，教えてください．

1階積分 (和分)過程 ( $m=1$ )．青実線が数値実験結果，赤破線が解析的結果．

1階積分 (和分)過程 ( $m=1$ )．青実線が数値実験結果，赤破線が解析的結果．

2階積分 (和分)過程 ( $m=2$ )．青実線が数値実験結果，赤破線が解析的結果．

2階積分 (和分)過程 ( $m=2$ )．青実線が数値実験結果，赤破線が解析的結果．

3階積分 (和分)過程 ([tex:m=13)．青実線が数値実験結果，赤破線が解析的結果．

　注意してほしいことは，この記事の内容は，長期記憶過程とか，弱定常の世界から，はみ出しているということです．枠にとらわれないことは重要ですので，以下の記事も参考にしてください．

chaos-kiyono.hatenablog.com

Rスクリプト

#########################
# m階積分(和分)
# m < 0のとき，m階差分
m <- 1
###########################
# 平均をとるサンプル時系列数
N.samps <- 100
##########################
# 時系列の長さ (時系列は奇数長になります)
N <- 1001
# 奇数に変換
N <- round(N/2)*2+1
#########################
# パラメタの設定
sig2 <- 1
sig <- sqrt(sig2)
#########################
# 数値実験
if(m < 0){
 N.sim <- N+2*m
}else{
 N.sim <- N
}
for(i in 1:N.samps){
　#########################
　# 白色ノイズの生成
　x.sim <- rnorm(N,0,sig)
　#########################

 if(m > 0){
  for(j in 1:m){
   # 【重要】 平均をきっちり0にしてから足さないと期待通りにならない
   x.sim <- x.sim-mean(x.sim)
   x.sim <- cumsum(x.sim)
  }
 }else if(m < 0){
  for(j in 1:(-m)){
   x.sim <- diff(x.sim)[1:(N-2*j)]
  }
 }
# x.sim <- x.sim-((x.sim[N]-x.sim[1])/(N-1)*(1:N)+x.sim[1])
#　plot(1:N,x.sim,"l")

　#########################
　# PSDの推定
　tmp <- fft(x.sim)
　PSD <- abs(tmp)^2/N
 if(i == 1){
  PSD.samps <- data.frame(PSD)
 }else{
  PSD.samps <- cbind(PSD.samps,data.frame(PSD))
 }
}
# サンプルの平均
PSD <- rowMeans(PSD.samps)
#########################
# 結果の描画
# モデルのパワースペクトルを与える
PSD.model <- function(f,m,sig2){
 return((f!=0)*sig2/(2*abs(sin(pi*f)))^(2*m))
}
#########################
# f = 0を避けて計算
M <- (N.sim-1)/2
f <- c((1:M)/N.sim,-(M:1)/N.sim)
# f = 0のとき，パワースペクトルを0に設定
PSD.model <- c(0,PSD.model(f,m,sig2))
f <- c(0,f)
##################################
par(mfrow=c(1,2),cex.lab=1.4)
par(pty="m")
plot(0:(N.sim-1),x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main=paste("Sample Time Series: m =",m),lwd=2,las=1)
#
par(pty="s")
plot(f[f > 0],PSD[f > 0],"l",log="xy",col=4,pch=16,xlim=c(f[2],0.5),xaxs="i",xlab="f",ylab="",main="Averaged Power Spectrum",lwd=4,las=1)
lines(f[f > 0],PSD.model[f > 0],col=2,lwd=2,lty=2)
legend("topright", legend=c(paste("Analytical: beta =",2*m)), col=c(2), lty=c(2),lwd=c(2))
axis(1,at=10^(-15:15)%x%(1:9),label=FALSE,tck=-0.015)
axis(2,at=10^(-15:15)%x%(1:9),label=FALSE,tck=-0.015)