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

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

【時系列のフラクタル解析】DMCAの実際(最適ラグの推定)

今回は,Detrending Moving-Average Cross-Correlation Analysis (DMCA)でラグを推定します.ここでは,DMCAC言語プログラムをコンパイルしたもの(各自コンパイルしてください),windowsコマンドプロンプト,Rを使います.

下記の論文"Generalized theory for detrending moving-average cross-correlation analysis: A practical guide"の解説シリーズです.
doi.org
DMCAC言語ファイルは以下のリンクからダウンロードできます:
https://www.dropbox.com/s/ymugqqan5g6trhl/DMCA.zip?dl=0

0次DMCAC言語プログラム"xDMA0Lags.c"をコンパイラしたものが,"xDMA0Lags.exe"だとして説明します.同様に,2次DMCAは"xDMA2Lags.exe",4次DMCAは"xDMA4Lags.exe"とします.

ここでは,自分でラグを設定し,生成した時系列を対象として,正しくラグが推定できるかどうかを検証していきます.分析の流れは以下です:

  1. 長時間相互相関時系列の生成(過去のブログを参考にしてください)
  2. コマンドプロンプトDMCAの実行
  3. RでDMCA結果を読み込んで最適ラグを推定
長時間相互相関時系列の生成

サンプル時系列の生成については,過去のブログを参考にしてください:
chaos-kiyono.hatenablog.com

コマンドプロンプトDMCAの実行

解析したい時系列が含まれるフォルダのアイコン上で「Shift」+「右クリック」->「PoweShellウィンドウをここで開く」を選択してください.

PoweShellが起動しますので,そこで

cmd

を実行してください.

ここでは,見出し行(ヘッダ)を含まず,2変量データが2列で書かれているファイル"samp_series.csv"を分析するとします.2変量間のラグは20とします.

0次DMCAについては,

xDMA0Lags.exe -c 1 2 -L 30 samp_series.csv > result_xdma0.csv

を実行すると,"result_xdma0.csv"にDMCAの結果が保存されます."-c 1 2"は,1列と2列を読み込む指定,"-L 30"はラグを-30から30の範囲で調べる指定です.

出力ファイル"result_xdma0.csv"の内容は以下です:

ラグ スケールs 相互F関数 補正s log相互F関数 相互相関係数 log自己F関数1 log自己F関数2
lag s F^2 (s) \tilde{s} \log_{10} F_{12}(\tilde{s}) \rho (\tilde{s}) \log_{10} F_{1} (\tilde{s}) \log_{10} F_{2} (\tilde{s})

(変数の定義は論文[doi.org/10.1016/j.csfx.2020.100022]参照)

ここでまでの,PowerShellでの操作は以下のようになります:

PowerShell上ではこんな感じです
RでDMCA結果を読み込んで最適ラグを推定

DMCA結果のファイルを読み込むために,"result_xdma0.csv"を保存したフォルダへのパスを

DIR <- "...... ここに指定 ......"

に指定してください.例えば,

D:\Document\時系列解析

の場合は,バックスラッシュ("\",日本語環境では¥マーク)を,"/"に書き直す必要がありますので,下記のRスクリプト2行目を

DIR <- "D:/Document/時系列解析"

に変更してください.さらに,読み込むファイル名"result_xdma0.csv"を,下記のRスクリプトの4行目に指定します.

FN <- "result_xdma0.csv"

以上の点に注意して下記のスクリプトを実行してください.

#  分析結果があるフォルダのパス
DIR <- "...... ここに指定 ......"
# xDMA*LagSweeperの分析結果ファイル名
FN <- "result_xdma0.csv"
#####################################
# ラグ推定後の結果出力ファイル名
FN.OUT <- "lag_optimized_result.csv"
#####################################
# 評価するスケール範囲
log_s1 <- 0.8
log_s2 <- 3
#####################################
setwd(DIR)
DAT <- read.csv(FN,header=FALSE)
#####################################
par(mfrow=c(2,2))
par(cex.lab=1.4,cex.axis=1.2)
#####################################
LAGs <- unique(DAT$V1)
N.LAGs <- length(LAGs) 
#########################################
# 相関係数の面積の計算
Area <- c()
for(i in 1:N.LAGs){
 lag <- LAGs[i]
 TMP <- DAT[DAT$V1==lag,]
 Area[i] <- abs(sum(TMP$V6[TMP$V4>=log_s1 & TMP$V4 <= log_s2][-1]*diff(TMP$V4[TMP$V4>=log_s1 & TMP$V4 <= log_s2])))
}
#########################################
# ラグの推定値
Lag.est <- LAGs[which.max(Area)]
#########################################
for(i in 1:N.LAGs){
 lag <- LAGs[i]
 TMP <- DAT[DAT$V1==lag,]
 if(i==1){
  plot(TMP$V4,TMP$V6,"l",col="lightpink",lwd=1,ylim=c(-1,1),xlab="log s",ylab="Cross Correlation",las=1,main=paste("Estimated Lag = ",Lag.est))
  abline(h=c(-1,0,1),lty=2,col=8)
  abline(v=c(log_s1,log_s2),lty=2,col=8)
 }else{
  lines(TMP$V4,TMP$V6,col="lightpink",lwd=1)
 }
 TMP <- DAT[DAT$V1==Lag.est,]
 lines(TMP$V4,TMP$V6,col=2,lwd=2)
}
plot(LAGs,Area,"l",col="lightpink",lwd=2,ylim=c(0,max(Area)*1.1),yaxs="i",xlab="Lag",ylab="Abs Corr Area",main=paste("Estimated Lag = ",Lag.est),las=1)
abline(v=Lag.est,col=2,lty=2,lwd=2)
#########################################
y.min <- min(TMP$V5)
y.max <- max(TMP$V5)
for(i in 1:N.LAGs){
 lag <- LAGs[i]
 TMP <- DAT[DAT$V1==lag,]
 if(i==1){
  plot(TMP$V4,TMP$V5,"l",col="lightpink",lwd=1,ylim=c(y.min,y.max),xlab="log s",ylab="log F12(s)",las=1,main=paste("Estimated Lag = ",Lag.est))
 }else{
  lines(TMP$V4,TMP$V5,col="lightpink",lwd=1)
 }
 TMP <- DAT[DAT$V1==Lag.est,]
 lines(TMP$V4,TMP$V5,col=2,lwd=2)
}
#########################################
TMP <- DAT[DAT$V1==Lag.est,]
y.min <- min(c(TMP$V5,TMP$V7,TMP$V8))
y.max <- max(c(TMP$V5,TMP$V7,TMP$V8))
plot(TMP$V4,TMP$V7,"l",col=3,lwd=2,ylim=c(y.min,y.max),xlab="log s",ylab="log F12(s)",las=1,main=paste("Estimated Lag = ",Lag.est))
lines(TMP$V4,TMP$V8,col=4,lwd=2)
lines(TMP$V4,TMP$V5,col=2,lwd=2)
#########################################
# 最適Lagの結果
DAT.OUT <- data.frame(log_s=TMP$V4,cor=TMP$V6,log_F12=TMP$V5,log_F1=TMP$V7,log_F2=TMP$V8)
write.table(DAT.OUT, FN.OUT, sep = ",", append=FALSE, quote=FALSE, col.names=TRUE,row.names=FALSE) 

きちんと分析出ていれば,下のような図が出力され,最適ラグの結果が"lag_optimized_result.csv"として保存されます(FN.OUTでファイル名を指定).

DMCAの結果の例

"lag_optimized_result.csv"の内容は

\tilde{s} \rho (\tilde{s}) \log_{10} F_{12}(\tilde{s}) \log_{10} F_{1} (\tilde{s}) \log_{10} F_{2} (\tilde{s})

です.