【時系列解析】ランダムウォーク解析の解析的考察（その３：パワースペクトルとの関係）

今回も，ランダムウォーク解析についてです．ランダムウォーク解析 (fluctuation analysisと呼ばれることも)は，Natureに出版された論文でも使われています (以下のリンク)．
www.nature.com
とはいえ，ランダムウォーク解析は，今では実用的な時系列解析法とはいえません．それでも，教育的な題材としては役に立ちます．

以下では，ゆらぎ関数 $F(s)$ とパワースペクトル $S(f)$ の関係を導いていきます．前回の記事は，これです．
chaos-kiyono.hatenablog.com

基本事項の確認

1. ラグオペレータ $L$ を，時系列 $\{x[i]\}$ に作用させると，

$\displaystyle L x [i] = x[i-1]$ ，　 $\displaystyle L^m x [i] = x[i-m]$

となります．つまり， $L$ が1つかかるごとに，時間が１戻ります．

2. 時系列 $\{x[i]\}$ のフーリエ変換を $X(f)$ とします．時間の世界でのラグオペレータ $L$ をフーリエ変換すると， $e^{-{\bf i } 2 \pi f}$ になります．つまり，

$\displaystyle L x [i]$ をフーリエ変換すると， $\displaystyle e^{{\bf i} 2 \pi f} X(f)$

となります．フーリエ変換すると， $L$ が $e^{{\bf i} 2 \pi f}$ に変わると覚えてください．

3. 時系列 $\{x[i]\}$ のパワースペクトルの推定量 $S_x(f)$ は，長さ $N$ の時系列 $\{x[i]\}$ のフーリエ変換を $X(f_k)$ (ここで， $f_k = k/N$ )として，

$\displaystyle S_x(f_k) = \frac{|X(f_k)|^2}{N}$

です．

4. $S_x(f_k)$ (ここで， $f_k = k/N$ )の全面積は，時系列 $\{x[i]\}$ の2乗平均 (平均0なら分散)になります．

$\displaystyle \sum_{k=0}^{N-1} S_x(f_k) \Delta f_k=\sum_{k=0}^{N-1} \frac{S_x(f_k)}{N} = \left\langle x[i]^2 \right\rangle$

これは，パーセヴァルの等式 (Parseval's identity)と呼ばれる関係から来ています．

ランダムウォーク解析を周波数の世界で考える

　ランダムウォーク解析では，平均0の観測時系列 $\left\{x[i] \right\}$ を積分し，

　 $\begin{equation} \displaystyle y[i] = \sum_{j=1}^{i} x[j] \end{equation}$

その増分

　 $\begin{eqnarray} \Delta_{s} y[i] &=& y[i+s] - y[i] = \sum_{j=1}^{s} x[i+j] \end{eqnarray}$

の2乗平均の期待値 (ゆらぎ関数: fluctuation function)

　 $\begin{equation} F^2(s) = \left\langle \Delta_{s} y^2[i] \right\rangle \end{equation}$

を求めるのでした．

　増分の計算の部分の計算から考えます．計算が楽になるように和をとる区間を１シフトします．つまり，「1からs」を，「0から (s-1)」にしちゃいました．そうすると，

　 $\begin{eqnarray} \Delta_{s} y[i] &=& \sum_{j=0}^{s-1} x[i+j] = x[i]+x[i+1] + \cdots + x[i+s-1] \end{eqnarray}$

です．これについて考えます．

ラグオペレータ $L$ を使えば，

　 $\begin{eqnarray} \Delta_{s} y[i] &=& x[i]+L^{-1} x[i] + \cdots + L^{-(s-1)} x[i] \\ &=& (1+L^{-1}+L^{-2}+\cdots + L^{-(s-1)} ) x[i] \\ &=& \frac{1-L^{s}}{1-L^{-1}} x[i] \end{eqnarray}$

です．最後の行になるのは，等比数列の和の公式を使ったからです．ちなみに等比数列の和の公式は，(初項－末項×公比)/(1－公比)でしたよね．

次に，両辺をフーリエ変換します． $\Delta_{s} y[i]$ のフーリエ変換を $Y_{\Delta s} (f_k)$ ， $x[i]$ のフーリエ変換を $X(f_k)$ と表せば，

　 $\begin{eqnarray} Y_{\Delta s} (f_k) &=& \frac{1-\left(e^{-{\bf i} 2 \pi f} \right)^{-s}}{1-\left(e^{-{\bf i} 2 \pi f} \right)^{-1}} X(f_k) \end{eqnarray}$

となります．これを使えば，パワースペクトルの関係式を求めることができます．つまり，「両辺の絶対値をとり，２乗して， $N$ でわる」の３ステップです．

　 $\begin{eqnarray} \left|Y_{\Delta s} (f_k)\right| &=& \left|\frac{1-\left(-e^{{\bf i} 2 \pi f_k} \right)^{-s}}{1-\left(-e^{{\bf i} 2 \pi f_k} \right)^{-1}} X(f_k) \right| \\ \left|Y_{\Delta s} (f_k)\right|^2 &=& \left|\frac{1-e^{-{\bf i} 2 \pi f_k s}}{1-e^{-{\bf i} 2 \pi f_k}} X(f_k) \right|^2 \\ \left|Y_{\Delta s} (f_k)\right|^2 &=& \left|\frac{1-e^{-{\bf i} 2 \pi f_k s}}{1-e^{-{\bf i} 2 \pi f_k}}\right|^2 \left|X(f_k) \right|^2 \\ \frac{\left|Y_{\Delta s} (f_k)\right|^2}{N} &=& \frac{\left|1-e^{-{\bf i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} \frac{\left|X(f_k) \right|^2}{N} \end{eqnarray}$

ここで， $\Delta_{s} y[i]$ のパワースペクトルを $S_{\Delta s} (f_k)$ ， $x[i]$ のパワースペクトルを $S(f_k)$ と表せば，

　 $\begin{eqnarray} \frac{\left|Y_{\Delta s} (f_k)\right|^2}{N} &=& \frac{\left|1-e^{-{\bf_k i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} \frac{\left|X(f_k) \right|^2}{N} S_{\Delta s} (f_k) &=& \frac{\left|1-e^{-{\bf i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} S(f_k) \\ \end{eqnarray}$

さらに，

　 $\begin{eqnarray} \frac{\left|1-e^{-{\bf i} 2 \pi f_k s}\right|^2}{\left|1-e^{-{\bf i} 2 \pi f_k}\right|^2} &=& \frac{(1-e^{-{\bf i} 2 \pi f_k s})(1-e^{{\bf i} 2 \pi f_k s})}{(1-e^{-{\bf i} 2 \pi f})(1-e^{{\bf i} 2 \pi f_k})} \\ &=& \frac{2-2 \cos(2 \pi f_k s)}{2 - 2 \cos (2 \pi f_k)} \\ &=& \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} \\ \end{eqnarray}$

となるので，結局，

　 $\begin{eqnarray} S_{\Delta s} (f_k) &=& \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \\ \end{eqnarray}$

　ここで，ゆらぎ関数

　 $\begin{equation} F^2(s) = \left\langle \Delta_{s} y^2[i] \right\rangle \end{equation}$

に登場してもらい，上の基本事項の4を使えば，

　 $\begin{equation} \displaystyle F^2(s) = \frac{1}{N}\sum_{k=0}^{N-1} \frac{\sin^2(\pi f_k s)}{\sin^2(\pi f_k)} S(f_k) \end{equation}$

という関係式が得られます．これで，時系列のパワースペクトルとゆらぎ関数を解析的に結びつけることができました．正しいか心配なので，数値実験してみました．下の図は，同じ時系列 (H = 0.8の非整数ガウスノイズ)のランダムウォーク解析 (左)とパワースペクトル解析 (中)をした結果です．右の図は，ランダムウォーク解析の結果 (青丸)と，上の式を使ってパワースペクトル解析の結果を $F(s)$ に変換したもの (赤バツ)の比較です．わずかに違いましたが，丸め誤差の影響だと思います．ポイントは，「ランダムウォーク解析とパワースペクトル解析は同じ情報を持っている」ということです．とはいえ，ランダムウォーク解析の方が，ギザギザしないので，直線の当てはめはやりやすいです．一方で，パワースペクトル解析の結果は，平滑化してやる必要があるということです．

今日は，明日締め切りの原稿があるので，ここまでにしておきます．ミスを見つけたら教えてください．

おまけ

　作図に使用したRスクリプトはこれです．

# パッケージのインストールが必要な場合は以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
###################################
# パラメタの設定：ハースト指数(0 < H < 1)
H <- 0.8
###################################
# 時系列の長さ
N <- 10000
######################
# 解析するスケールs
s <- unique(round(exp(seq(log(1),log(N/10),length.out=20))))
n.s <- length(s)
######################
# 時系列の生成
x <- simFGN0(N,H)
x <- x - mean(x) # 平均0に
#####################
# ランダムウォーク解析
###
# 積分
y <- cumsum(x)
###
F.s <- c()
for(i in 1:n.s){
 D.y <- y[1:(N-s[i])]-y[(1+s[i]):N]
 F.s[i] <- sqrt(mean(D.y^2))
}
freq <- (0:(N-1))/N
S.f <- abs(fft(x))^2/N
par(mfrow=c(1,3))
###
# log-logプロット
par(pty="s")
# log-logをとって直線をあてはめ
logs <- log10(s)
logFs <- log10(F.s)
fit <- lm(logFs~logs)
#
plot(logs,logFs,col=4,xlab="log10 s",ylab="log10 F(s)",las=1,
 main="Random Walk Analysis")
abline(fit,col=1,lty=2,lwd=2)
# log-logをとって直線をあてはめ
fit <- lm(log10(S.f[freq>0 & freq<1/10])~log10(freq[freq>0 & freq<1/10]))
#
plot(log10(freq[freq>0 & freq<1/2]),log10(S.f[freq>0 & freq<1/2]),"l",col=2,xlab="log10 f",ylab="log10 S(f)",las=1,
 main="Spectral Analysis")
abline(fit,col=1,lty=2,lwd=2)

F2.Sf <- c()
for(i in 1:n.s){
 F2.Sf[i] <- 0
 for(k in 1:N){
  fk <- k/N
  F2.Sf[i] <- F2.Sf[i]+sin(pi*fk*s[i])^2/sin(pi*fk)^2*S.f[k]
 }
}
plot(logs,logFs,col=4,xlab="log10 s",ylab="log10 F(s)",las=1,main="Comparison")
points(log10(s),log10(F2.Sf/N)/2,pch=4,col=2,xlab="log10 s",ylab="log10 F(s)",las=1)