今回は，長時間相関，長期記憶，フラクタル性などを示す時系列のスケーリング解析をやってみようというお話です．
以下では，Detrending Moving Average Algorithm，略して，DMAと呼ばれる方法を紹介します．この業界では，Detrended Fluctuation Analysis，略して，DFAと呼ばれる方法が有名ですが，DFAは数理的に洗練されていない方法です．DFAの欠点としては，非線形フィルタを使っているためバイアスが生じる (解析結果がガタガタする)，自己相関関数やパワースペクトルとの関係が無駄に複雑ということがありますが，DMAではこれらの欠点はありません．とはいえ，DMAを使おうが，DFAを使おうが，多くのユーザにとっては，結果はほとんどかわらないので，DFAを使っても問題ありません．
今回はやってみるだけ
ここでは，長時間相関をもつ時系列を数値的に生成して，その時系列をDMAで解析してみます．Rで以下のスクリプトを実行してみてください．
######################
# 事前に以下の2つのパッケージをインストールすること
# install.packages("longmemo")
# install.packages("signal")
######################
require(longmemo)
require(signal)
######################
# DMAの次数
q <- 0
######################
# 時系列の長さ
n <- 10000
# Hurst指数の指定
H <- 0.8
# 時系列の分散
StDev <- 1
################
# scaleは平均0，標準偏差1にする関数
x <- scale(simFGN0(n,H))[,1]*StDev
#
############################
# スケールを何点とるか
n.s <- 30
# スケールの決定 (DMAで解析するスケールは奇数のみ)
s <- unique(round(exp(seq(log(q+3),log(n/4),length.out=n.s))/2)*2+1)
# スケール数の再計算 (重複があるので減る)
n.s <- length(s)
#########################
F2 <- c()
# 【STEP1】時系列の積分
y <- cumsum(x)
# q次DMA
for(i in 1:n.s){
# 【STEP2】Detrending
y.detrend <- y - sgolayfilt(y,p=q,n=s[i],m=0,ts=1)
# 【STEP3】2乗偏差の計算
F2[i] <- mean(y.detrend^2)
}
# スケーリング
log10F.s <- log10(F2)/2
log10F.min <- min(log10F.s)
log10F.max <- max(log10F.s)
#
par(mfcol=c(2,2),las=1,cex.axis=1.2,cex.lab=1.5,mar=c(5,5,3,1),cex.main=1.6)
par(pty="m")
plot(1:n-1,x,"l",col=4,xaxs="i",xlab="i",main="Sample time series",ylab="x[i]")
plot(1:n-1,y,"l",col=4,xaxs="i",xlab="i",main="Integrated series",ylab="x[i]")
# 【STEP4】両対数プロットで傾きを計算
plot(log10(s),log10F.s,col=4,pch=16,xlab="log10 s",ylab="log10 F(s)",ylim=c(log10F.min,log10F.max),main="DMA")
###
# 傾きの推定
log10s.min <- 1 # 区間の下限
log10s.max <- 3 # 区間の上限
lfit <- lm(log10F.s[log10(s) >= log10s.min & log10(s) <= log10s.max]~log10(s)[log10(s) >= log10s.min & log10(s) <= log10s.max])$coefficients
curve(lfit[[1]]+lfit[[2]]*x,xlim=c(log10s.min,log10s.max),add=TRUE,col=2,lty=2,lwd=2)
abline(v=c(log10s.min,log10s.max),col=gray(0.6),lty=2)
legend("bottomright",legend=c(sprintf("Estimated slope:%.2f",lfit[[2]])),lty=c(2),col=c(2),lwd=c(2))
###
plot(log10(s)[-1],diff(log10F.s)/diff(log10(s)),col=4,pch=16,xlab="log10 s",ylab="local slope",ylim=c(0,2),main="local slope estimation")
abline(h=H,lty=2,col=gray(0.6))
legend("topleft",legend=c("Theory","Slope between adjacent points"),lty=c(2,NA),pch=c(NA,16),col=c(gray(0.6),4))
結果は，
皆さんの手元にある時系列を分析したいときは，
x <- (時系列ベクトル)
の部分にデータを渡してください．
Rスクリプトの解説
サンプル時系列の生成では，
# Hurst指数の指定 H <- 0.8
の部分の数値を変えると，Hurst指数を指定できます．Hは，0より大きくて，1より小さい値にしてください．
DMAで解析するスケールは，奇数である必要があります．ですので，
# スケールの決定 (DMAで解析するスケールは奇数のみ) s <- unique(round(exp(seq(log(q+3),log(n/4),length.out=n.s))/2)*2+1)
のように，2の倍数に1を加える計算をしています．
DMAの本体は，
# 【STEP1】時系列の積分
y <- cumsum(x)
# q次DMA
for(i in 1:n.s){
# 【STEP2】Detrending
y.detrend <- y - sgolayfilt(y,p=q,n=s[i],m=0,ts=1)
# 【STEP3】2乗偏差の計算
F2[i] <- mean(y.detrend^2)
}
の部分です．DMAでは，トレンド除去にSavitzky-Golayフィルタを使います．DFAではこんなに簡単に書けません．
最後に両対数プロットを描き，傾きを計算しているのは，
# 【STEP4】両対数プロットで傾きを計算 plot(log10(s),log10F.s,col=4,pch=16,xlab="log10 s",ylab="log10 F(s)",ylim=c(log10F.min,log10F.max),main="DMA")
から下の部分です．
注意：DMAには高速アルゴリズムがある
今回はRでDMAをやってみました．最近はパソコンが早くなっているので，DMAの計算をRでも実行できます (計算時間が気にならないという意味です)．しかし，時系列が長くなったりデータ数が増えると，Rなんかではやってられません．その場合は，DMAの高速アルゴリズムを使う必要があります． 以下のリンクに高速アルゴリズムの論文があります．必要であれば，C言語のプログラムを提供しますので連絡してください．
