ケィオスの時系列解析メモランダム

時系列解析,生体情報学,数学・物理などの解説です.

【時系列解析】ランダムウォーク解析の解析的考察(その1:白色ノイズの場合)

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

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)