【Rで時系列解析】中点変位法 (midpoint displacement algorithm)で非整数ブラウン運動のサンプルパスを生成

１次元の非整数ブラウン運動 (fractional Brownian motion)のサンプルパス (時系列)の生成方法の一つ「中点変位法」(midpoint displacement algorithm)について説明します．

この話は，「フラクタル時系列って何？」，「非整数ブラウン運動のフラクタル次元は？」，「非整数ガウスノイズってフラクタル？」，「ピンクノイズ (1/fノイズ)はフラクタル？」という疑問につながっていきます．今日，Wikipedaで「ピンクノイズ」を調べてみると，
ja.wikipedia.org
そこに，

『ピンクノイズの波形は、フラクタル状になっていることが知られている。』

と書いてありました．『フラクタル状』ってことは，フラクタルだなと理解する人もいると思います．正直私も，「ピンクノイズの波形は、フラクタル的」と説明してしまうことがあります．本当はフラクタルとはいえないかもという気持ちで，「フラクタル的」と言っています．というのは，フラクタルっぽいルールでピンクノイズの波形を構成できるけど，そのフラクタル次元は２で，フラクタルの特徴である非整数次元を持ちません．Wikipedaで「ピンクノイズ」の説明を書いた方はどのような意味で書いたのか気になりました．前置きが長くなりましたが，今回は，フラクタル次元の２歩手前として，中点変位法について説明します．

非整数ブラウン運動のオリジナルの論文は，
　[Mandelbrot, Benoit B., and John W. Van Ness. "Fractional Brownian motions, fractional noises and applications." SIAM review 10, 422-437 (1968)]
中点変位法の論文は，
　[Fournier et al., Communications of the ACM, 25, 371 (1982)]
です．

非整数ブラウン運動の特徴

　中点変位法の根拠となる，非整数ブラウン運動の２つ特徴を説明します．
１．非整数ブラウン運動 $B_H(t)$ の増分の２乗偏差の平方根 (標準偏差のようなもの)は，増分をとる時間幅が $s$ のとき， $s^{H}$ に比例する．
　下の図のように，非整数ブラウン運動 $B_H(t)$ のサンプルパスをたくさん描くと，時間の経過と共に軌道が広がっていきます．この広がり具合が $s^{H}$ に比例するということです (標準偏差とはっきり言えないのは，この時系列は非定常 (幅が発散)なので，標準偏差は存在しないからです)．

つまり，

　 $\displaystyle \sqrt{{\rm E} \left( \left(B_H(t+s)-B_H(t) \right)^2 \right)} \propto s^{H}$

ということです．

２．非整数ブラウン運動の軌道 $B_H(t)$ 上の２点間には，見かけの直線トレンドがある．
　下の図のように，２点 $(0, B_H(0))$ ， $(1, B_H(1))$ をつなぐ，非整数ブラウン運動の軌道 $B_H(t)$ をたくさん考えると，２点を結んだ線分上を通る確率がもっとも高いということです．

中点変位法 (midpoint displacement algorithm)

　以下の手順で，サンプルパスを作ります．非整数ブラウン運動のハースト指数を $H$ とします．

１．平均0，分散1の正規乱数の値を２つ生成し (もちろん独立)，それを， $B_H(0)$ ， $B_H(1)$ とする (両方0にしたり，片方0にしても構いません)．

２．２点 $(0, B_H(0))$ ， $(1, B_H(1))$ の中点に，平均0，標準偏差 $2^{-H}$ の正規乱数を足す (足すのは縦軸 $B_H$ 方向です)．

３．隣り合う２点の間に中点をとる．この中点と最っとも近い隣の点までの距離が $s$ のとき，平均0，標準偏差 $s^{H}$ の正規乱数を足す．この作業を繰り返す．

詳しく説明するのが面倒になったので，サンプルのRスクリプトを下に載せておきます．

Rスクリプト

################
# ハースト指数
H <- 0.8
################
# 分割する回数
# 15~18くらいを超えると計算時間がかかります
n.step <- 14
################
# 最初の２点
x <- c(0,1)
y <- rnorm(2) 
# 0から出発するときは，y <- c(0,rnorm(1))
# 両端を0にするときは，y <- c(0,0)
################
i <- 2
for(j in 1:n.step){
 h <- (1/2)^j
 x.tmp <- x[1:i]
 y.tmp <- y[1:i]
 x[(1:(i-1))*2] <- (x.tmp[1:(i-1)]+x.tmp[2:i])/2
 y[(1:(i-1))*2] <- (y.tmp[1:(i-1)]+y.tmp[2:i])/2+h^(H)*rnorm(i-1)
 x[(1:i)*2-1] <- x.tmp
 y[(1:i)*2-1] <- y.tmp
 i <- i + (i-1)
}
# グラフの描画
plot(x,y,"l",col=4,las=1,lwd=2,main=paste("fBm with H =",H))