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

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

【Rで時系列解析】Detrending moving average algorithm (DMA)を使った長時間相関スケーリング解析

今回はDetrending moving average algorithm (DMA)をRでやってみます.DMAの参考文献はこれです(ダウンロードできない人が多いと思います.オープンアクセスの良い論文を見つけられませんでした):
Detrending Moving Average algorithm: a brief review | IEEE Conference Publication | IEEE Xplore

ここでは,signalパッケージのSavitzky-Golayフィルタを使って,積分時系列の移動平均を計算し,この移動平均をトレンドとして除きます.

まずは,解析する時系列の準備です.これまで何度かお世話になっているlongmemoパッケージを使います.分析したい時系列がある場合は,このRスクリプトは無視して,xに時系列をベクトルとして格納してください.

# パッケージのインストールが必要な場合は以下を実行
# install.packages("longmemo")
# パッケージのロード
require(longmemo)
###################################
# パラメタの設定:ハースト指数(0 < H < 1)
H <- 0.6
# 時系列の長さ
n <- 2^10
####################################
# 非整数ガウスノイズの生成 
x <- simFGN0(n,H)

以下は,DMAの説明用のグラフ作成も一緒にやるRスクリプトです.

#
# xを時系列ベクトルとする
#
N <- length(x)
###############################
# Savitzky-Golay フィルタにsignalパッケージを使用 
# install.packages("signal")
require(signal)
###############################
# 平均0化
x <- x - mean(x)
# 積分
y <- cumsum(x)
###############################
# 作図用のパラメタ
par(mfrow=c(3,2))
par(pty="m",cex=1.1,las=1,cex.lab=1.2)
par(mar=c(5,4,2,2))
# 解析対象の時系列
plot(1:N,x,"l",col=4,xaxs="i",xlab="i",ylab="x[i]",lwd=2,main="Original time series")
plot(1:N,y,"l",col=4,xaxs="i",xlab="i",ylab="y[i]",lwd=2,main="Integrated series")
###############################
# scalesの設定
scale.min <- 4
scale.max <- N/4
scale <- unique(round(10^seq(log10(scale.min),log10(scale.max),0.05)/2)*2+1)
n.scale <- length(scale)
# プロットする3つのスケールの指定
p.scale <- c(3,round(3+(n.scale-3)/2),n.scale-2)
###############################
F2 <- c()
for(i in 1:n.scale){
 F2[i] <- 0
 i.edge <- (scale[i]-1)/2
 y.sg <- sgolayfilt(y,p=0,n=scale[i],m=0,ts=1)
 if(is.element(i, p.scale)){
  plot(1:N,y,"l",col=4,xaxs="i",lwd=2,xlab="i",ylab="y[i]",main=paste("scale =",scale[i]))
  lines((i.edge+1):(N-i.edge),y.sg[(i.edge+1):(N-i.edge)],col=2,lty=2,lwd=2)
 }
 F2[i] <- mean((y[(i.edge+1):(N-i.edge)]-y.sg[(i.edge+1):(N-i.edge)])^2)
}
###############################
# 結果のプロット
par(pty="s")
fit <- lm(log10(F2)/2~log10(scale))
plot(log10(scale),log10(F2)/2,xlab="log10 n",ylab="log10 F(n)",main=paste("slope =",format(fit$coefficients[[2]],digit=2)))
points(log10(scale[p.scale]),log10(F2[p.scale])/2,pch=16,col=2)
abline(fit,lty=2,col=2,lwd=2)

このスクリプトを実行すると,下のような図が出力されます.

DMAの結果の例

上の図では,左上に解析する時系列x,その右が積分時系列,下の時系列は,各スケールで区分的に直線 (赤実線)を当てはめた結果,右下がスケーリング解析(\log_{10} nに対する\log_{10} F(n))の結果です.

実際の時系列解析では,こんなに図は必要ないので,スクリプトはこんな感じになります.

#
# xを時系列ベクトルとする
#
N <- length(x)
###############################
# Savitzky-Golay フィルタにsignalパッケージを使用 
# install.packages("signal")
require(signal)
###############################
# 平均0化
x <- x - mean(x)
# 積分
y <- cumsum(x)
###############################
# scalesの設定
scale.min <- 4
scale.max <- N/4
scale <- unique(round(10^seq(log10(scale.min),log10(scale.max),0.05)/2)*2+1)
n.scale <- length(scale)
###############################
F2 <- c()
for(i in 1:n.scale){
 F2[i] <- 0
 i.edge <- (scale[i]-1)/2
 y.sg <- sgolayfilt(y,p=0,n=scale[i],m=0,ts=1)
 F2[i] <- mean((y[(i.edge+1):(N-i.edge)]-y.sg[(i.edge+1):(N-i.edge)])^2)
}
###############################
# 結果のプロット
par(pty="s")
fit <- lm(log10(F2)/2~log10(scale))
plot(log10(scale),log10(F2)/2,xlab="log10 n",ylab="log10 F(n)",main=paste("slope =",format(fit$coefficients[[2]],digit=2)))
points(log10(scale[p.scale]),log10(F2[p.scale])/2,pch=16,col=2)
abline(fit,lty=2,col=2,lwd=2)