【Rで時系列解析】1/fゆらぎを積分すると傾き2だけ変化（数値実験）

前回は，時系列の積分がパワースペクトルにもたらす変化

「 $1/f^{\beta}$ ゆらぎを積分すると， $1/f^{\beta+2}$ ゆらぎになる」

を解析的に示しました．

chaos-kiyono.hatenablog.com

今回は，数値実験です．

非整数ガウスノイズ（左）とその積分（右），下は対応するパワースペクトルのlog-logプロット．

使ったRスクリプトは以下です．パワースペクトルは，N.sampsで指定した数のサンプル時系列の平均で求めています．時系列は，1例を示しているだけです．

# longmemoがないときは以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
####################################
# ハースト指数(0 < H < 1)
H <- 0.6
# 時系列の長さ
n <- 2^10
####################################
# サンプル数
N.samps <- 100
####
S.x.mean <- numeric(n)
S.y.mean <- numeric(n)
for(i in 1:N.samps){
 ####################################
 # 非整数ガウスノイズの生成 
 x <- simFGN0(n,H)
 ####################################
 # 積分 (xの平均をゼロにする必要がある．なーんでか？)
 x <- x-mean(x)
 y <- cumsum(x)
 ####################################
 # フーリエ変換
 X <- fft(x)
 Y <- fft(y)
 # パワースペクトルの推定
 S.x <- abs(X)^2/n
 S.y <- abs(Y)^2/n
 #
 S.x.mean <- S.x.mean+S.x
 S.y.mean <- S.y.mean+S.y
}
S.x.mean <- S.x.mean/N.samps
S.y.mean <- S.y.mean/N.samps
################################################
f <- (0:(n-1))/n
log.S.x <- log10(S.x.mean[f > 0 & f <= 1/2])
log.S.y <- log10(S.y.mean[f > 0 & f <= 1/2])
y1 <- min(log.S.x,log.S.y)
y2 <- max(log.S.y)
log.f <- log10(f[f > 0 & f <= 1/2])
# 直線フィット
fit.x <- lm(log.S.x[log.f < -1]~log.f[log.f < -1])
fit.y <- lm(log.S.y[log.f < -1]~log.f[log.f < -1])
# プロット
par(mfrow=c(2,2))
par(pty="m")
plot(0:(n-1),x,"l",col=2,xlab="i",ylab="x[i]",xaxs="i")
plot(0:(n-1),y,"l",col=4,xlab="i",ylab="y[i]",xaxs="i")
par(pty="s")
plot(log.f,log10(S.x.mean[f > 0 & f <= 1/2]),col=2,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sx(f)",xaxs="i",main=paste("slope =",format(fit.x$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.x,lty=2,col=1,lwd=2)
plot(log.f,log10(S.y.mean[f > 0 & f <= 1/2]),col=4,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sy(f)",xaxs="i",main=paste("slope =",format(fit.y$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.y,lty=2,col=1,lwd=2)
########################################

ちなみに，上のRスクリプトで

 ####################################
 # 積分 (xの平均をゼロにする必要がある．なーんでか？)
 x <- x-mean(x)

の部分にある
x <- x-mean(x)
を消して，H>0.5のときの数値実験をしてみてください．積分時系列の傾きが-2になってしまうと思います．なぜでしょうか？

上の図を作成したRスクリプトです．

# longmemoがないときは以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
####################################
# ハースト指数(0 < H < 1)
H <- 0.9
# 時系列の長さ
n <- 2^10
####################################
# サンプル数
N.samps <- 100
####
S.x.mean <- numeric(n)
S.y.mean <- numeric(n)
for(i in 1:N.samps){
 ####################################
 # 非整数ガウスノイズの生成 
 x <- simFGN0(n,H)
 ####################################
 # 積分 (xの平均をゼロにする必要がある．なーんでか？)
 y <- cumsum(x)
 ####################################
 # フーリエ変換
 X <- fft(x)
 Y <- fft(y)
 # パワースペクトルの推定
 S.x <- abs(X)^2/n
 S.y <- abs(Y)^2/n
 #
 S.x.mean <- S.x.mean+S.x
 S.y.mean <- S.y.mean+S.y
}
S.x.mean <- S.x.mean/N.samps
S.y.mean <- S.y.mean/N.samps
################################################
f <- (0:(n-1))/n
log.S.x <- log10(S.x.mean[f > 0 & f <= 1/2])
log.S.y <- log10(S.y.mean[f > 0 & f <= 1/2])
y1 <- min(log.S.x,log.S.y)
y2 <- max(log.S.y)
log.f <- log10(f[f > 0 & f <= 1/2])
# 直線フィット
fit.x <- lm(log.S.x[log.f < -1]~log.f[log.f < -1])
fit.y <- lm(log.S.y[log.f < -1]~log.f[log.f < -1])
# プロット
par(mfrow=c(2,2))
par(pty="m")
plot(0:(n-1),x,"l",col=2,xlab="i",ylab="x[i]",xaxs="i")
plot(0:(n-1),y,"l",col=4,xlab="i",ylab="y[i]",xaxs="i")
par(pty="s")
plot(log.f,log10(S.x.mean[f > 0 & f <= 1/2]),col=2,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sx(f)",xaxs="i",main=paste("slope =",format(fit.x$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.x,lty=2,col=1,lwd=2)
plot(log.f,log10(S.y.mean[f > 0 & f <= 1/2]),col=4,pch=16,"l",lwd=2,ylim=c(y1,y2),las=1,
  xlab="log f",ylab="log Sy(f)",xaxs="i",main=paste("slope =",format(fit.y$coefficients[[2]],digit=3)),xaxs="i")
abline(fit.y,lty=2,col=1,lwd=2)
########################################