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

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

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

0次DMCAのC言語プログラム"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

(変数の定義は論文[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"の内容は

です．