時系列解析の手法のひとつに、「マルチフラクタル解析」という、どこか響きのかっこいいものがあります。ただしこの解析は、マルチフラクタル性という特殊な性質を持つ時系列にしか適用できません。ですから、私はあまり積極的に使おうとは思わないのです。しかも時系列が短いと、推定精度がとんでもなく悪くなる。これもまた、私がマルチフラクタル解析をいまひとつ信用しきれない理由です。
研究の世界では、ときどき難しそうな言葉を並べて相手を煙に巻こうとする研究者がいます。もしあなたがフラクタル時系列の論文を投稿したら、査読者から「マルチがどうのこうの」といったコメントを受けるかもしれません。ところが残念なことに、そうした査読者の多くがマルチフラクタル性を本当に理解しているかというと、正直ちょっと怪しいのです。
少し嫌味っぽく（むしろ偉そうに）「マルチフラクタル」という言葉を出しましたが、もちろんこれは重要で面白い概念でもあります。せっかくなので、今後は、私なりにわかりやすく説明してみようと思います。今回はその前段として、マルチフラクタル性を持つ時系列をどうやって生成できるのか、実際の方法を紹介してみましょう。
マルチフラクタル時系列は分散のゆらぎが1/f
1/fゆらぎとは、時系列のパワースペクトルが周波数 f に反比例する「1/f 型」になる性質のことです。数理的にはかなり特異な確率過程なのですが、なぜか自然界や生体の信号には、この1/fゆらぎがやたらと顔を出します。ただし、すべてのゆらぎが1/fになるわけではありません。ですので、1/fゆらぎは“特別な”性質です。
マルチフラクタルではない、いわゆるモノフラクタルな時系列では、値の分布は正規分布（ガウス分布）に従います。例外としてレビー安定分布などもありますが、実世界で頻繁に登場するわけではありません。
ところが、マルチフラクタルな時系列は少し様子が違います。局所的に分散がゆらぐため、値の分布は非ガウス的になります。そして何より特異なのは、その局所分散の変動が1/fゆらぎの性質をもっている、という点です。ここで多くの人が勘違いするのですが、単に局所分散がゆらぐだけではマルチフラクタルではありません。それが、1/fゆらぎになってはじめて「マルチフラクタル」なのです。
この基本的な性質を踏まえて、マルチフラクタル時系列を生成するモデルのひとつが マルチフラクタル・ランダムウォーク（Multifractal Random Walk, MRW） です。今回紹介するのは、その増分過程をシミュレーションするためのRスクリプトです。スクリプト中の x が増分過程に対応し、それを積分して得られる X が、MRW のサンプル軌道になります。――少々風変わりな数理的構築ですが、これをいじってみると、マルチフラクタルの“実体”が少し身近に感じられるかもしれません。ただし、今回は簡単化するために、増分が白色過程になるマルチフラクタル・ランダムウォークを生成します。
Multifractal random walk | Phys. Rev. E
Rスクリプト
## ------------------------------------------------------------
## Multifractal Random Walk (MRW) sample path generator
## based on Bacry–Delour–Muzy (PRE 2001) and related notes.
## This implementation is restricted to the H = 1/2 case,
## i.e., the driving process ε_t is Gaussian white noise.
## ------------------------------------------------------------
set.seed(100) # 再現性のための乱数種
## -----------------------
## Model parameters
## -----------------------
L <- 2048 # integral scale T に相当（相関が消えるスケール）
lambda2 <- 0.06 # intermittency係数 λ^2
dt <- 1 # 離散サンプリング間隔 l に相当（以下の共分散に使う）
N <- 100001 # 生成する系列長
N <- 2 * round(N/2) + 1
M <- (N - 1) / 2 # 半長（周回埋込みのため）
## ------------------------------------------------------------
## (A) MRW の「対数ボラティリティ」ω(t) の理論自己共分散（Bacry ら）
## 離散化版：τ=0.. での ρ(τ) を以下で定義（PRE 2001 の離散近似）。
## 連続時間の定義（例：E[ω] と ρ(τ) の形）は文献の式(5)(6)を参照。
## - E[ω] = -Var(ω) = -λ^2{ log(T/l) + 1 }
## - ρ(τ) = Cov(ω_t, ω_{t+τ}) は
## 0 ≤ τ < l : λ^2{ log(T/l) + 1 - τ/l }
## l ≤ τ < T : λ^2 log(T/τ)
## それ以外 : 0
## ここでは離散化して l = dt, T = L として実装。
## ------------------------------------------------------------
mrw_covariance <- function(tau, L, lambda2, l = 1) {
a <- abs(tau)
out <- numeric(length(a))
# τ = 0
out[a == 0] <- lambda2 * (log(L / l) + 1)
# 0 < τ < l
idx1 <- (a > 0) & (a < l)
out[idx1] <- lambda2 * (log(L / l) + 1 - a[idx1] / l)
# l ≤ τ < L
idx2 <- (a >= l) & (a < L)
out[idx2] <- lambda2 * log(L / a[idx2])
# τ ≥ L は 0 のまま
out
}
## ------------------------------------------------------------
## (B) 周回埋込み（circulant embedding）用の共分散ベクトルを構築
## 長さ N の巡回共分散 [0,1,...,M, M,...,1] を組み立てる
## ------------------------------------------------------------
lags_two_sided <- c(0:M, -(M:1)) # 0,1,...,M, -M,...,-1
cov_omega_circ <- mrw_covariance(lags_two_sided, L, lambda2, l = dt)
## パワースペクトル（固有値）— 非負性を数値的に確保
PSD <- Re(fft(cov_omega_circ))
PSD <- pmax(PSD, 0) # 丸め誤差に伴う微小負値を切り上げ
## ------------------------------------------------------------
## (C) 角周波数領域で ω を生成（スペクトル合成）
## 実数時系列を得るため、FFT(ホワイトノイズ)を振幅整形
## 注: N は奇数なのでナイキスト点はありません
## ------------------------------------------------------------
WN <- rnorm(N) # 実白色雑音（時間領域）
WN_fft <- fft(WN) # そのフーリエ変換（共役対称を持つ）
omega_fft <- sqrt(PSD) * WN_fft
omega_raw <- Re(fft(omega_fft, inverse = TRUE)) / N # 時間領域へ
## ------------------------------------------------------------
## (D) 理論平均の設定（E[exp(2ω)] = 1 となる条件）
## PRE 2001 / Bacry ノートの定義に合わせて：
## E[ω] = -Var(ω) = -λ^2{ log(T/l) + 1 }
## 生成した ω をゼロ平均化したうえで、理論平均を付与。
## ------------------------------------------------------------
mean_theory <- -lambda2 * (log(L / dt) + 1)
omega_sim <- omega_raw - mean(omega_raw) + mean_theory
## ------------------------------------------------------------
## (E) MRW の「インクリメント」x_t = ε_t * exp(ω_t) を生成（式(7)）
## ε_t ~ N(0,1) を独立（ω と無相関）にとる。
## 上の平均条件により E[exp(2ω)] = 1 なので Var(x_t) = Var(ε_t) = 1
## ------------------------------------------------------------
eps <- rnorm(N) # ガウス白色雑音（イノベーション）
x <- eps * exp(omega_sim)
## （必要なら）MRW の「累積和」X_t = Σ x_t を作る：
X <- cumsum(x) # 価格や位置の「標本軌道」を見たい場合
## ------------------------------------------------------------
## (F) 幾何級数ラグで ω の自己共分散を推定し、理論と比較
## ------------------------------------------------------------
autocov_at_lag <- function(z, k) {
n <- length(z)
mu <- mean(z)
sum((z[1:(n - k)] - mu) * (z[(1 + k):n] - mu)) / n
}
# 幾何ラグ（1,2,4,...）と 0 を含む等比系列
lags_geo <- c(0, unique(round(exp(seq(log(1), log(L - 1), length.out = 50)))))
cov_emp <- sapply(lags_geo, function(k) autocov_at_lag(omega_sim, k))
cov_th <- mrw_covariance(lags_geo, L, lambda2, l = dt)
## ------------------------------------------------------------
## (G) 可視化
## ------------------------------------------------------------
op <- par(mfrow = c(3, 1), xaxs = "i", mar = c(5, 5, 3, 2))
# (1) MRW インクリメント（あるいは X を見たければ x -> X）
plot(seq_len(N) - 1, x, type = "l", col = 4, main = "MRW increments x_t = ε_t exp(ω_t)",
xlab = "t", ylab = "x_t")
# (2) ω の共分散（通常軸）
plot(lags_geo, cov_emp, type = "l", col = 4, lwd = 2,
ylim = c(0, lambda2 * (log(L / dt) + 1)),
main = "Cov(ω_t, ω_{t+τ}) vs. τ (linear scale)",
xlab = "lag τ", ylab = "covariance")
abline(h = 0, lty = 2, col = gray(0.7))
lines(lags_geo, cov_th, col = 2, lty = 2, lwd = 2)
legend("topright", c("empirical", "theory"), col = c(4, 2), lwd = 2, lty = c(1, 2), bty = "n")
# (3) ω の共分散（横軸のみ対数）
idx <- lags_geo > 0
plot(lags_geo[idx], cov_emp[idx], type = "l", col = 4, lwd = 2, log = "x",
ylim = c(0, lambda2 * (log(L / dt) + 1)),
main = "Cov(ω_t, ω_{t+τ}) vs. τ (log-x scale)",
xlab = "lag τ (log scale)", ylab = "covariance")
abline(h = 0, lty = 2, col = gray(0.7))
lines(lags_geo[idx], cov_th[idx], col = 2, lty = 2, lwd = 2)
par(op)
