前回に続き，ARFIMA(0, d, 0)を，無限次の自己回帰過程 (AR過程: autoregressive process)で表現する話です．今回は数値実験です．
前回の記事は，これです．
chaos-kiyono.hatenablog.com
時系列を生成するのが目的であれば，ARFIMAの自己共分散から生成するとか，より洗練された方法があると思います．しかしここでは，あえて，自己回帰過程の差分方程式を使って，時系列を生成します．
1. 基本事項：long AR過程表現
論文を探していると，long AR近似という表現がありました．今回は，有限だけどかなり長い自己回帰過程で，ARFIMA(0, d, 0)を近似します．使う，基本事項は以下です．
を平均0と分散
の白色ノイズとすれば，ARFIMA(0, d, 0)は以下のように表すことができます．
2. 数値計算のポイント
ガンマ関数を使った数値実験では，式の通りにガンマ関数の値を計算してもうまくいきません．なぜなら，ガンマ関数の値が，無茶苦茶でかくなるので，コンピュータで計算すると表現できる最大値を超えてしまします．例えば，Rで，
gamma(x)
とすれば，の値を計算してくれます．
を計算してみると，
> gamma(171) [1] 7.257416e+306
となり，になってます．
で限界を迎え，
> gamma(172) [1] Inf
のようにRの世界では，大きすぎて限界を超えた値 (Inf)になります．
ですので，数値計算では，の代わりに
の値を教えてくれる関数を使うことが多いです．Rでは，
lgamma(x)というコマンドです．
今回はポイントは，ガンマ関数の部分を
として計算することです．
3. 数値計算結果
Rスクリプトは，一番下に載せておきます．Rスクリプトでは，求めたいサンプル時系列の長さNに対して，ARFIMA(0,d,0)をN次自己回帰過程で近似しています．計算では，合計3N回計算を繰り返して，前から2N分を過渡状態として捨てています．
結果の図をいくつか載せておきます．時系列の長さをとし，200本のサンプル時系列のパワースペクトルの平均を計算しました．図中の赤破線が解析的な結果です．もちろん，ARFIMA(0,d,0)では許されない，
の結果も載せておきます．
Rスクリプト
###########################
# 平均をとるサンプル時系列数
N.samps <- 100
##########################
# 時系列の長さ
N <- 1000
#########################
# パラメタの設定
d <- 0.4
sig2 <- 1
sig <- sqrt(sig2)
#########################
# 数値実験
j <- 1:N
coef <- sign(d)*exp(lgamma(j-d)-lgamma(j+1)-lgamma(-d))
for(i in 1:N.samps){
#########################
if(d == 0){
x.sim <- rnorm(N,0,sig)
}else{
x0 <- numeric(N)
eps <- rnorm(2*N,0,sig)
for(k in 1:(2*N)){
x0[N+k] <- x0[(N+k-1):k] %*% coef + eps[k]
}
x.sim <- x0[(2*N+1):(3*N)]
}
#########################
# 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.theol <- function(f,m,sig2){
return((f!=0)*sig2/(2*abs(sin(pi*f)))^(2*m))
}
#########################
# 周波数の設定とPSDの解析的結果
f <- c((1:floor(N/2))/N,-(floor(N/2):1)/N)
PSD.ana <- c(0,PSD.theol(f,d,sig2))
f <- c(0,f)
#########################
par(mfrow=c(1,3),cex.lab=1.4)
plot(0:N,c(1,coef),"l",col=4,lwd=2,xlab="j",ylab="AR coefficients",xaxs="i",las=1)
abline(h=0,lty=2,col=8)
plot(0:(N-1),x.sim,"l",xaxs="i",col=4,xlab="i",ylab="x[i]",main=paste("Synthetic Time Series: d =",d),lwd=2,las=1)
#
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="Power Spectrum",lwd=4,las=1)
lines(f[f > 0],PSD.ana[f > 0],col=2,lwd=2,lty=2)
legend("topright", legend=c(paste("Analytical PSD with slope of -",2*d)), 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)
