Detrending Moving-Average Algorithm、あるいは Detrending Moving-Average Analysis はDMAと略され、1 本の時系列における長時間自己相関や自己アフィン性を評価するための解析手法です。それに対して、Detrending Moving-Average Cross-Correlation Analysisは、2 本の時系列に対して長時間相互相関を評価するための解析手法であり、DMCA と略されます。本記事では、この DMCA について解説していきます。とくに、ここでは我々の研究グループが提案している、Savitzky–Golay フィルタを用いてトレンド除去を行うDMCA に焦点を当てます。
オリジナルの DMA および DMCA は、生体信号時系列などに含まれる非定常トレンドを移動平均フィルタによって除去することを基本的な考え方としています。移動平均フィルタは、部分時系列に定数関数をフィットし、その値を用いて平滑化する方法と解釈することができます。
それに対してSavitzky–Golayフィルタは、下図のように部分時系列に多項式関数をフィットし、その区間中央の値を用いて平滑化する手法です。その意味で、Savitzky–Golay フィルタは移動平均フィルタの自然な一般化とみなすことができます。そのような背景から、我々は移動平均の代わりに Savitzky–Golay フィルタを用いる DMA および DMCA の数理的基礎を体系的に構築してきました。
本稿ではまず、DMCA によって、時系列のどのような特性を評価できるのかという点に焦点を当てて解説します。
を示し、右端のパネルは長時間相関特性の解析結果を示している。
参考文献 doi.org
1.２本の時系列に隠された共通成分を検出
DMCA の代表的な適用例の一つは、2 本の時系列に共通して含まれる長時間相関成分の存在を検出することです。たとえば、下図のような状況を想定します。
は、共通の入力
の影響を受けている。
上図は、2 本の観測時系列が、それぞれ固有のノイズ成分に加えて、共通の長時間相関成分を含んでいる状況を模式的に示しています。左側の は、非整数ガウスノイズ (fGn)に代表される長時間相関をもつ共通信号成分であり、これが 2 本の経路に分岐して、それぞれの時系列に同時に加わっていす。
一方、 および
は、各時系列に固有の独立なノイズ成分を表しており、それぞれが加算された結果として、観測時系列
が得られます。このような場合、個々の時系列を単独で解析しても、共通成分 の存在はノイズに埋もれてしまい、明確には捉えにくくなります。
DMCA を用いると、2 本の時系列を同時に扱うことで、両者に共通して含まれる長時間相関成分のみが強調され、図に示した に由来する相互相関のスケーリング特性を検出することが可能になります。この仕組みが、DMCA が「2 本の時系列に隠れた共通ダイナミクス」を明らかにできる理由です。
2. RでDMCAの数値実験
上記の過程 を数値的に生成し、実際にDMCAを適用した結果をRで確認するスクリプトを以下に示しました。
■ この数値実験で何を確認しているのか
このスクリプトの目的は、2 本の時系列に「共通の長時間相関成分」が弱く混入しているときに、DMCA によってその共通成分の存在とスケーリング特性が可視化できることを、数値実験で確かめることです。
2.1 想定しているモデル
数値データは、以下のモデルに基づいて生成されます。
- 共通成分：長時間相関をもつ確率過程（本研究では非整数ガウスノイズ成分）
- 個別成分：各時系列に固有のノイズ（こちらも長時間相関をもつことを許します）
これらが加算されたものとして観測されるモデルです。
Rスクリプトでは、次のように 2 本の系列を作っています：
- 共通成分：
x_common - 系列 1：
x1 = (1-r_c)*noise1 + r_c*x_common - 系列 2：
x2 = (1-r_c)*noise2 + r_c*x_common
ここで r_c は「共通成分の混合比」で、r_c が小さいほど共通成分は埋もれ、単独解析では見えにくくなります。
2.2 DMCA で評価する量
DMCAでは、スケール（窓長） を変えながらトレンド除去後の「2 本の時系列の共通の変動の大きさ」を
により評価し、そのスケーリング
特性から、共通成分の長時間相互相関特性を評価します。また、共通成分の相互相関の強さを 多重スケール相互相関係数
で評価します。
したがって、DMCA では以下の統計量に注目します。
(1) 多重スケール相互相関係数 
多重スケール相互相関係数
はスケール ごとに「相互相関がどれだけ強いか」を表す量です。
なら「そのスケールでは相互相関がない」ということですので、共通成分の存在は確認できません。
が大きく、１に近いほど「そのスケールで共通成分が存在している」という解釈になります。
(2) ゆらぎ関数のスケーリング
F1(s)：時系列 1 のゆらぎの大きさ（自己相関と関連）F2(s)：時系列２のゆらぎの大きさ（自己相関と関連）F12(s)：２本の時系列に共通のゆらぎの大きさ（相互相関と関連）
を対数スケールでプロットし、log10(s) に対する傾きを推定しています。
この傾きは「長時間自己相関と長時間相互相関」を特徴づける指標になります。
2.3 Rスクリプト
############################################################
# Detrending Moving-Average Cross-Correlation (DMCA) demo
# - Simulate two series sharing a common long-memory component
# - Compute DMA/DMCA-like fluctuation functions using SG detrending
# - Add axis labels + legends
# - In bottom-right panel: estimate slopes in log10(s) in [1,2]
# and overlay fitted line segments
############################################################
## ---- utilities ----
check_pkg <- function(pkg) {
if (!requireNamespace(pkg, quietly = TRUE)) {
stop(sprintf("Package '%s' is required. Install it via install.packages('%s').", pkg, pkg),
call. = FALSE)
}
}
# fGn-like Gaussian noise via ARFIMA(0,d,0) simulation using fracdiff
# d = H - 0.5
sim_fgn_fracdiff <- function(n, H, sd = 1, mu = 0, seed = NULL) {
if (!is.null(seed)) set.seed(seed)
stopifnot(n >= 2, is.finite(H), H > 0, H < 1, is.finite(sd), sd > 0)
check_pkg("fracdiff")
d <- H - 0.5
x <- fracdiff::fracdiff.sim(n = n, d = d)$series
x <- as.numeric(scale(x)) # standardize
x <- mu + sd * x # set mean/sd
x
}
# compute fluctuation functions for given scales s (odd integers)
# SG detrending on integrated series (y = cumsum(x))
compute_dma_dmca <- function(x1, x2, s, sg_order = 1) {
check_pkg("signal")
n <- length(x1)
stopifnot(length(x2) == n, all(s >= 3), all(s %% 2 == 1))
y1 <- cumsum(x1)
y2 <- cumsum(x2)
F1 <- numeric(length(s))
F2 <- numeric(length(s))
F12 <- numeric(length(s)) # signed cross-deviation
for (i in seq_along(s)) {
si <- s[i]
# SG smoothing of integrated series, then detrend
y1_tr <- signal::sgolayfilt(y1, p = sg_order, n = si, m = 0, ts = 1)
y2_tr <- signal::sgolayfilt(y2, p = sg_order, n = si, m = 0, ts = 1)
r1 <- y1 - y1_tr
r2 <- y2 - y2_tr
F1[i] <- sqrt(mean(r1^2))
F2[i] <- sqrt(mean(r2^2))
F12[i] <- mean(r1 * r2) # can be negative
}
rho <- F12 / (F1 * F2)
list(
s = s,
log10s = log10(s),
F1 = F1,
F2 = F2,
F12 = F12,
rho = rho
)
}
# fit slope in log-log space in a given log10(s) range and return model + segment coords
fit_slope_segment <- function(x, y, xlim = c(1, 2)) {
ok <- is.finite(x) & is.finite(y) & (x >= xlim[1]) & (x <= xlim[2])
if (sum(ok) < 2) return(NULL)
# avoid "x[ok]" inside formula; use a data.frame with clean names
dat <- data.frame(x = x[ok], y = y[ok])
fit <- lm(y ~ x, data = dat)
xx <- range(dat$x)
yy <- predict(fit, newdata = data.frame(x = xx))
list(
fit = fit,
xseg = xx,
yseg = yy,
slope = unname(coef(fit)[["x"]]),
intercept = unname(coef(fit)[["(Intercept)"]])
)
}
## ---- simulate example ----
n <- 1000
Hc <- 0.7
r_c <- 0.5
x_common <- sim_fgn_fracdiff(n, Hc, sd = 1, mu = 0, seed = 1)
H1 <- 0.85
H2 <- 0.45
x1 <- as.numeric(scale(sim_fgn_fracdiff(n, H1, seed = 13) * (1 - r_c) + r_c * x_common))
x2 <- as.numeric(scale(sim_fgn_fracdiff(n, H2, seed = 15) * (1 - r_c) + r_c * x_common))
## ---- choose scales (odd only) ----
n_s <- 30
s <- unique(round(exp(seq(log(5), log(n / 10), length.out = n_s)) / 2) * 2 + 1)
s <- s[s %% 2 == 1]
log10s <- log10(s)
## ---- compute DMA/DMCA-like measures ----
res <- compute_dma_dmca(x1, x2, s = s, sg_order = 1)
log10F1 <- log10(res$F1)
log10F2 <- log10(res$F2)
log10F12 <- 0.5 * log10(abs(res$F12)) # follows your original definition
y_min <- floor(min(log10F1, log10F2, log10F12, na.rm = TRUE) * 5) / 5
y_max <- ceiling(max(log10F1, log10F2, log10F12, na.rm = TRUE) * 5) / 5
## ---- plot ----
op <- par(no.readonly = TRUE)
on.exit(par(op), add = TRUE)
par(mfrow = c(2, 2), pty = "m", xaxs = "i", yaxs = "i", mar = c(4, 4, 2.5, 1))
ii <- seq_len(n)
# (1) x1 with common component
plot(ii, x1, type = "l", lwd = 2, col = 3,
xlab = "Time index i", ylab = "Amplitude (standardized)",
main = "Series 1 and common component",ylim=c(-4,4))
lines(ii, r_c * x_common, lwd = 2, col = rgb(1, 0.1, 0.3, alpha = 0.4))
legend("topright",
legend = c("x1", "x_common"),
col = c(3, rgb(1, 0.1, 0.3, alpha = 0.4)),
lwd = c(2, 2), bty = "n")
# (2) x2 with common component
plot(ii, x2, type = "l", lwd = 2, col = 4,
xlab = "Time index i", ylab = "Amplitude (standardized)",
main = "Series 2 and common component",ylim=c(-4,4))
lines(ii, r_c * x_common, lwd = 2, col = rgb(1, 0.1, 0.3, alpha = 0.4))
legend("topright",
legend = c("x2", "x_common"),
col = c(4, rgb(1, 0.1, 0.3, alpha = 0.4)),
lwd = c(2, 2), bty = "n")
# (3) multiscale cross-correlation coefficient rho(s)
par(xaxs = "r") # allow a little margin for readability here
plot(res$log10s, res$rho, type = "p", pch = 16, col=2,
xlab = expression(log[10](s)),
ylab = expression(rho[12](s)),
ylim = c(-1, 1),
main = "Multiscale cross-correlation")
abline(h = c(0), lty = 2, col = gray(0.6))
legend("topright",
legend = expression(rho[12](s)),
pch = 16, bty = "n")
# (4) log-log fluctuation functions with slope estimates on x in [1,2]
par(xaxs = "r")
plot(res$log10s, log10F1, type = "p", pch = 1, col = 3,
xlab = expression(log[10](s)),
ylab = expression(log[10](F)),
ylim = c(y_min, y_max),
main = expression("Fluctuation functions: " * log[10](F) * " vs " * log[10](s)))
points(res$log10s, log10F2, pch = 1, col = 4)
points(res$log10s, log10F12, pch = 16, col = 2)
# fit slopes in xlim = c(1,2) and overlay segments
xfit_rng <- c(1, 2)
fit1 <- fit_slope_segment(res$log10s, log10F1, xlim = xfit_rng)
fit2 <- fit_slope_segment(res$log10s, log10F2, xlim = xfit_rng)
fit12 <- fit_slope_segment(res$log10s, log10F12, xlim = xfit_rng)
if (!is.null(fit1)) segments(fit1$xseg[1], fit1$yseg[1], fit1$xseg[2], fit1$yseg[2], col = 3, lwd = 2)
if (!is.null(fit2)) segments(fit2$xseg[1], fit2$yseg[1], fit2$xseg[2], fit2$yseg[2], col = 4, lwd = 2)
if (!is.null(fit12)) segments(fit12$xseg[1], fit12$yseg[1], fit12$xseg[2], fit12$yseg[2], col = 2, lwd = 2)
# draw the fitting window for clarity
abline(v = xfit_rng, lty = 3, col = gray(0.6))
# legend with slope values
leg_txt <- c(
sprintf("F1 (slope=%.3f)", ifelse(is.null(fit1), NA, fit1$slope)),
sprintf("F2 (slope=%.3f)", ifelse(is.null(fit2), NA, fit2$slope)),
sprintf("F12 (slope=%.3f)", ifelse(is.null(fit12), NA, fit12$slope))
)
legend("topleft",
legend = leg_txt,
col = c(3, 4, 2),
pch = c(1, 1, 16),
lty = c(NA, NA, NA),
lwd = c(NA, NA, NA),
bty = "n")
Rスクリプトについて
ここからは、Rスクリプトのブロックごとに「何をしているか」を説明します。
(1) 依存パッケージのチェック：check_pkg()
check_pkg <- function(pkg) { if (!requireNamespace(pkg, quietly = TRUE)) { stop(sprintf("Package '%s' is required. ...", pkg), call. = FALSE) } }
この関数は、「必要なパッケージが入っていなければ、分かりやすいエラーメッセージで止める」ための安全装置です。 このスクリプトでは
fracdiff（ARFIMA のシミュレーション）signal（Savitzky–Golay フィルタ）
を使うので、ここでチェックしています。
(2) 非整数ガウスノイズのサンプル時系列の生成：sim_fgn_fracdiff()
sim_fgn_fracdiff <- function(n, H, ...) { d <- H - 0.5 x <- fracdiff::fracdiff.sim(n = n, d = d)$series x <- as.numeric(scale(x)) x <- mu + sd * x x }
ここは「長時間相関をもつノイズ」を作る部分です。
- ARFIMA(0,d,0) を用いて系列を生成
d = H - 0.5によって、Hurst 指数を直感的なパラメタとして指定
- 生成後に
scale()で標準化して平均・分散を整える
という処理をしています。
(3) DMCA の核：compute_dma_dmca()
y1 <- cumsum(x1) y2 <- cumsum(x2)
まず積分時系列（profile）を作ります。
次にスケール ごとに
y1_tr <- signal::sgolayfilt(y1, p = sg_order, n = si, ...) r1 <- y1 - y1_tr
として、Savitzky–Golay 平滑化をトレンド推定とみなし、残差 r1 を作ります（時系列 2 も同様）。
最後に
F1[i] <- sqrt(mean(r1^2)) F2[i] <- sqrt(mean(r2^2)) F12[i] <- mean(r1 * r2) rho <- F12 / (F1 * F2)
を計算します。
F1,F2：自己相関を評価するためのゆらぎ関度F12：相互相関を評価するためのゆらぎ関度rho：多重スケール相互相関
ここが本スクリプトの中心です。
注意： 論文の厳密な処理では、２本の時系列間のラグ (\kappa) の推定が必要ですが、ここでは、その部分を省略し、簡略化したDMCAのデモになっています。
(4) 傾き推定：fit_slope_segment()
ok <- (x >= 1) & (x <= 2) fit <- lm(y ~ x, data = dat)
右下パネルで
log10(s)の範囲を [1,2] に限定し- その範囲の点だけで直線回帰して傾きを推定
- 推定した直線を「線分」として重ね描き
します。
これは「スケーリング指数を目視で読む」のではなく、指定したスケール領域で定量化するためです。
(5) 可視化
左上・右上のパネルについては、時系列と共通成分
x1,x2を描き- 半透明で
r_c*x_commonを重ねています。
左下のパネルについては、 を描いています。この結果から、
- どのスケールで相互相関が強いか
- スケールに依存して相互相関が変化するか
を見ます。水平線 h=0 は「相互相関なし」の基準です。
右下のパネルは、 のスケーリング（べき則）を評価しています。ここでは
log10(F) vs log10(s) を描き、さらに [1,2] の区間で傾きを推定します。
F1,F2の傾き：それぞれの時系列の長時間自己相関特性F12の傾き：共通成分（相互相関）の長時間相互相関特性
パラメタを変えて数値実験
まずは、DMCAの挙動を理解するために、上のRスクリプト中の
r_c <- 0.5
の部分を
r_c <- 0
として、共通成分が存在しない場合の結果と比較してみてください。r_c の値を大きくすると、共通成分の寄与が大きくなりますが、時系列が短いと検出精度は安定しません。時系列の長さ指定
n <- 1000
の値を、小さくしたり、大きくしたりして、時系列の長さの影響も考察してください。
r_c <- 0 とした実行例
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
