【Ｒを使った時系列解析】非整数ガウスノイズの生成 - ケィオスの時系列解析メモランダム

今回は，Ｒを使って非整数ガウスノイズ（fractional Gaussian noise）を生成します．しかも，一切，パッケージに頼らずにです．

そもそも，非整数ガウスノイズとは何かというと，ポイントは以下です：

ガウスノイズなので，分布は正規分布に従います．
自己共分散関数が，パラメタH（ハースト指数）を使って

　　 $C[k] = \frac{1}{2} \left( \left|k+1 \right|^{2H} - 2 \left|k \right|^{2H} + \left|k - 1 \right|^{2H} \right)$
　で与えられます．

この情報だけで，平均0の非整数ガウスノイズのサンプル時系列は生成できます．方法の概略は以前の記事
chaos-kiyono.hatenablog.com
を参照してください．

ハースト指数H，長さN (奇数)の，サンプル時系列x.simを作るＲスクリプトは以下です．

# 時系列の長さ
N <- 1001
# 偶数はいやなので，念のため奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# モデルの自己共分散関数を与える
# 例：非整数ガウスノイズの自己共分関数
AutoCov.model <- function(k,H){
 return((abs(k+1)^(2*H)-2*abs(k)^(2*H)+abs(k-1)^(2*H))/2)
}
# パラメタの設定
H <- 0.8
# 自己共分散の値を計算
acov.model <- c(AutoCov.model(0:M,H),AutoCov.model(M:1,H))
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

これで，正しくできているのかを確かめるためには，サンプル時系列を1本作ってもだめです．たくさんのサンプル時系列から求めた平均のパワースペクトルや自己共分散関数を，理論と比べる必要があります．

ということで，常に数値的に確認することが私のモットーです．使用するＲスクリプトはこれです．
自分で設定を変更してほしいパラメタは，ハースト指数H，時系列の長さN，サンプル数N.ensです．

# 時系列の長さ
N <- 1001
# 奇数に変換
N <- round(N/2)*2+1
# 半分
M <- (N-1)/2
#########################
# モデルの自己共分散関数を与える
# 例：非整数ガウスノイズの自己共分関数
AutoCov.model <- function(k,H){
 return((abs(k+1)^(2*H)-2*abs(k)^(2*H)+abs(k-1)^(2*H))/2)
}
# パラメタの設定
H <- 0.9
sig2 <- 1
# 【注意】[-M,M]区間ではなく，[0,2M-1]にしている
acov.model <- c(AutoCov.model(0:M,H),AutoCov.model(M:1,H))
lag <- c(0:M,-(M:1))
#########################
# 自己共分散関数のフーリエ変換
fft.model <- fft(acov.model)
PSD.model <- Re(fft.model)
f <- c((0:M)/N,-(M:1)/(N))
#########################
# グラフ用パラメタ設定
par(mfrow=c(2,2),pty="m",cex=1.1)
##############################################
# サンプル時系列の総数
N.ens <- 500
PSD <- c()
ACOR <- c() 
for(i in 1:N.ens){
#########################
# 白色ノイズの生成
WN <- rnorm(N)
#　ホワイトノイズのフーリエ変換
fft.WN <- fft(WN)
#########################
fft.sim <-  sqrt(PSD.model)*fft.WN
#########################
x.sim <- Re(fft(fft.sim,inverse=TRUE))/N
###############################################
TMP <- acf(x.sim,plot=FALSE,type="covariance",lag.max=round(N/10))
Ac.y <- TMP$acf[,,1]
X <- fft(x.sim)
S.y <- abs(X)^2/N
 if(i == 1){
  PSD <- S.y
  ACOR <- Ac.y
 }else{
  PSD <- PSD+S.y
  ACOR <- ACOR + Ac.y
 }
}
PSD <- PSD/N.ens
ACOR <- ACOR/N.ens
################################################
plot(f[order(f)],PSD[order(f)],log="y",col=4,pch=16,"l",lwd=2,xlim=c(-0.5,0.5),xaxs="i",xlab="f",ylab="PSD")
lines(f[order(f)],PSD.model[order(f)],col=2,lwd=2,lty=2)
legend("topright", legend=c("Target Power Spectrum"), col=c(2), lty=c(2),lwd=c(2))
plot(f[f>0],PSD[f>0],log="xy",col=4,pch=16,"l",lwd=2,xlab="f",ylab="PSD")
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))
########################################
k.lag <- 1:length(ACOR)-1
plot(k.lag,ACOR,pch=16,col=4,"o",lwd=2,xlim=c(0,N/10),xaxs="i",xlab="k",ylab="Autocovariance")
curve(H*(2*H-1)*abs(x)^{2*(H-1)},add=TRUE,col=8,lwd=2)
lines(lag[order(lag)],acov.model[order(lag)],col=2,lty=2,lwd=2)
legend("topright", legend=c("Target Autocovariance"), col=c(2), lty=c(2),lwd=c(2))
plot(k.lag[k.lag>0],abs(ACOR[k.lag>0]),log="xy",pch=16,col=4,"o",lwd=2,xlab="k",ylab="Autocovariance")
curve(H*(2*H-1)*abs(x)^{2*(H-1)},add=TRUE,col=8,lwd=2)
lines(lag[lag>0],abs(acov.model[lag>0]),col=2,lty=2,lwd=2)
legend("topright", legend=c("Target Autocovariance"), col=c(2), lty=c(2),lwd=c(2))

結果は以下に並べてあります．ここで，わかることは，