【時系列解析】ランダムウォーク解析の解析的考察（その１：白色ノイズの場合）

ランダムウォーク解析のつづきです．今回は，この解析法について解析的に考えてみます．

chaos-kiyono.hatenablog.com

白色ノイズを分析すると

　白色ノイズは，無相関な時系列です．平均0，分散 $\sigma^2$ の白色ノイズの時系列 $\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}$

は，どんな感じでしょうか．

計算してみます．今回は，確率変数について大文字，小文字の区別はしません．全部，小文字です．弱定常であること，つまり，時間シフトしても，平均と自己共分散が変らないことを使えば，

　 $\begin{eqnarray} F^2(s) &=& \left\langle \Delta_{s} y^2[i] \right\rangle = \displaystyle \left\langle \left(\sum_{j=1}^{s} x[i+j]\right)^2 \right\rangle \\ &=& \left\langle x[i+1]^2 + x[i+2]^2 + \cdots + 2 \left( x[i+1] x[i+2]+ \cdots + \right) \right\rangle \\ &=& s \left\langle x[i]^2 \right\rangle + 2 \sum_{k=1}^{s-1} (s-k) \left\langle x[i]x[i+k] \right\rangle \end{eqnarray}$

となります．白色ノイズの場合， $\left\langle x[i]^2 \right\rangle = \sigma^2$ ， $\left\langle x[i]x[i+k] \right\rangle = 0$ ( $k \neq 0$ )なので，

　 $\begin{eqnarray} F^2(s) &=& \sigma^2 s \\ F(s) &=& \sigma s^{1/2} \\ \end{eqnarray}$

になります．ですので， $\log s$ に対して $\log F(s)$ をプロットして，傾きを推定すると $1/2=0.5$ になるはずです．

今回の計算には，弱定常の性質と，和の期待値の性質を使いました．

chaos-kiyono.hatenablog.com

次回の考察では，ラグオペレータに活躍してもらいますので，復習しておいてください．

chaos-kiyono.hatenablog.com

おまけ

図の作成に使ったRスクリプトです．

# 時系列の長さ
N <- 1000
# サンプル数
N.samp <- 100
######################
# 解析するスケールs
s <- unique(round(exp(seq(log(1),log(N/10),length.out=20))))
n.s <- length(s)
###
# ゆらぎ関数の準備
F.s <- c()
logFs.table <- data.frame()

par(mfrow=c(1,3))
for(j in 1:N.samp){
 #####################
 # サンプルの時系列 （白色ノイズ）
 # S(f) ~ f^0 (beta=0)の例
 x <- rnorm(N)
 x <- x - mean(x) # 平均0に
 #####################
 # ランダムウォーク解析
 ###
 # 積分
 y <- cumsum(x)
 ###
 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))
 }
 if(j == 1){
  par(pty="m")
  plot(1:N,x,"l",xaxs="i",xlim=c(0,N),main="サンプル時系列 (解析対象)",col=4,xlab="i")
  plot(1:N,y,"l",xaxs="i",xlim=c(0,N),main="積分時系列",col=4,xlab="i")
  logFs.table <- data.frame(log10(F.s))
 }else{
  logFs.table <- cbind(logFs.table,data.frame(log10(F.s)))
 }
}
logFs <- rowMeans(logFs.table)
# log-logをとって直線をあてはめ
logs <- log10(s)
fit <- lm(logFs~logs)
###
# log-logプロット
par(pty="s")
plot(logs,logFs,col=4,xlab="log10 s",ylab="log10 F(s)",las=1,
 main=paste("解析結果： 推定スロープ =",format(fit$coefficients[[2]],digits=2)))
# 直線フィット
abline(fit,col=2,lty=2)