時系列データの長期相関（long-range correlation）やフラクタル性を調べる代表的な手法として、DFA（Detrended Fluctuation Analysis） と DMA（Detrending Moving Average Analysis） が広く用いられています。一般には DFA の方が圧倒的に有名ですが、DFA にはいくつかの弱点があり、その点を補う方法として DMA を位置づけることができます。
今回は、この DMA を題材として、これまで解説してきた FIR フィルタ の知識が DMA の理解にどのように役立つのかを紹介します。
- 1. DMAの大部分はFIRフィルタ処理
- 2. 各ステップをFIRフィルタとして書いてみる
- 3. すべてを「1 本の FIR 畳み込み」にまとめる
- 4.【Ｒで体験】DMAをデジタルフィルタで実現
- 5. まとめ
1. DMAの大部分はFIRフィルタ処理
DMAでは、次のような処理をします。
(1) 分析したい時系列 からその平均を引く (必須ではありません)。ここでは、平均 0 の時系列を とします。
(2) の 累積和（積分）を計算する：
(3) 各スケール（窓長 ）で、積分時系列
に対し、フィルタ長
のSavitzky–Golayフィルタをかけて、その成分を、
をから引いた残差を求める。
(4) スケール ごとに、残差の二乗平均の平方根であるゆらぎ関数
を計算する。
(5) スケール に対する、
のべき則（スケーリング関係）、
から を推定する (両対数プロットに直線をフィットし傾きを推定)。
ここで、これまでデジタルフィルタの基礎に興味をもって学んでくれた皆さんには、上の太字の部分の処理はすべて、FIRフィルタで書けることに気づくと思います。つまり、
- 累積和（プロファイル）をとる
- Savitzky–Golay フィルター（多項式フィット）で局所トレンドを推定
- プロファイルからそのトレンドを引く
は、すべて畳み込み演算の形の線形処理であり、FIR フィルタとして表現できるのです。
したがって、
DMA の前処理全体を「1 本の畳み込み（FIR フィルタ）」としてまとめることができる
のです。
2. 各ステップをFIRフィルタとして書いてみる
ということで、実際に各ステップのインパルス応答を書き下してみましょう。
■ 累積和（プロファイル）
累積和
は、「今までの値を全部足す」という操作なので、一見すると IIR の無限和にも見えますが、DMA では 有限の窓幅の中で この累積を扱います。
たとえば、窓
の中だけで積分を考えると、時系列ベクトル
に対し、累積ベクトル
は
と書けます。ここで、 は
- 対角を含めた 下三角部分がすべて 1
- その他は 0
という行列
です。
■ Savitzky–GolayフィルタもFIRフィルタ
Savitzky–Golayフィルタは、窓内のデータ に対して
という 次多項式を最小二乗でフィットし、中心点
の値
を、「平滑化された値（局所トレンド）」として出力します。
最小二乗解は
で、中心値
は
と書けます。
ここで、
：デザイン行列（
）
：定数項を取り出すベクトル
です。
たとえば、窓長 （つまり
）、多項式次数
の場合、
の行列要素は
であり、 とすれば、
となります。これが Savitzky–Golayフィルタの「多項式フィッティング」を表すデザイン行列です。
さらに、多項式係数ベクトルを
とすると、（定数項）を取り出したいので、ベクトル
は、
となります。こうすることで、
の形で、定数項だけを抽出できます。
ここからがポイントです。窓の大きさを とし、インデックス
の順に並べたとします。このとき、窓内のデータを
Savitzky–Golayフィルタの係数ベクトルを
と書くことにします。
上の式
は、ある行ベクトルを、 にかけているだけなので、
とおけば、
と書けます。
つまり、Savitzky–Golayの出力 は、窓内のデータ
に対する重み付き和（FIRフィルタ）になっています。係数
は、上の式から
として、行列計算から一意に決まります。
■ 残差（プロファイルからトレンドを引く）も FIR フィルタ
DMA で使用するのは、平滑化された値そのものではなく、
という 残差 です。
中心付き窓 の中心を 0 としたとき、中心点の残差
をベクトルで書くと：
です。これは
と書き直すことができて、畳み込みの形
で表すことができます。ここで、は、
の 番目の成分です。
3. すべてを「1 本の FIR 畳み込み」にまとめる
ここまでをまとめると：
- 累積和：下三角 1 の行列
による線形変換
- Savitzky–Golay 多項式フィット：
による線形変換
- プロファイルからトレンドを引く：
による線形結合
これらはすべて 線形演算 なので、行列を掛け合わせれば
という形にまとめられます。
つまり、係数ベクトル を用いて、
と書くことができ、一般の時刻 についても
と書けます。これはまさに、インパルス応答 をもつ 1 本の FIR フィルタによる畳み込みとして、DMA の前処理全体が記述できることを意味しています。
次の節に、この係数 を具体的に計算するRスクリプトを掲載しておきます。
を求めることで、DMAが時系列の何を評価しているのかを理解することができます。
4.【Ｒで体験】DMAをデジタルフィルタで実現
■ FIRフィルタの係数を求める関数とDMAの計算を行う関数
まずは、以下のＲスクリプトを実行してください。
############################################################
# 1. FIR 設計関数
# Savitzky–Golay 型 DMA 前処理の「1ウィンドウ分」の FIR を求める
#
# 目的：
# - 時系列 x[n] から平均を引き、累積和（プロファイル）をとった後、
# 窓 [-M,...,M] 上の累積系列に p 次多項式をフィットし、
# 中央点の値からトレンド（多項式）を引く、という処理を
# 1 本の FIR 畳み込みとして表現する。
############################################################
sg_dma_local_fir <- function(M, p) {
#---------------------------------------------
# 引数チェック
#---------------------------------------------
if (M < 1) stop("M must be >= 1.")
if (p < 0) stop("p must be >= 0.")
# 窓長 L = 2M+1, インデックス j = -M, ..., M
L <- 2 * M + 1
j <- -M:M
#---------------------------------------------
# デザイン行列 A (L x (p+1))
# - 行 : インデックス j
# - 列 : 多項式の次数 0,...,p
# - 要素 : A_{j,m} = j^m
# これで、ベクトル a = (a_0,...,a_p)^T に対して
# q(j) = A %*% a
# が、j 上の多項式値を与える。
#---------------------------------------------
A <- outer(j, 0:p, function(jj, mm) jj^mm)
#---------------------------------------------
# 累積和の行列 B (L x L)
# - y = B x とすると、
# y[k] = sum_{r <= k} x[r]
# を実現する。
# - 行列としては「下三角が 1 の行列」。
#---------------------------------------------
B <- matrix(0, nrow = L, ncol = L)
B[lower.tri(B, diag = TRUE)] <- 1
#---------------------------------------------
# 中央点 j = 0 を取り出すためのベクトル c
# - j = -M,...,M を並べたとき、
# j = 0 は位置 M+1 に対応する。
# - y は長さ L のベクトルなので、
# c^T y = y[ j=0 ] を実現する。
#---------------------------------------------
c_vec <- rep(0, L)
c_vec[M + 1] <- 1 # j = 0 に対応
#---------------------------------------------
# 多項式係数 a から定数項 a_0 を取り出すベクトル e
# - a = (a_0, a_1, ..., a_p)^T に対して、
# e^T a = a_0 となる。
#---------------------------------------------
e_vec <- c(1, rep(0, p))
#---------------------------------------------
# 最小二乗解 a = (A^T A)^{-1} A^T y
# を使うための行列 (A^T A)^{-1} A^T を前計算する。
# - これは (p+1) x L の行列。
# - 入力 y (長さ L) に対して、
# a = ATA_inv_AT %*% y
# となる。
#---------------------------------------------
ATA_inv_AT <- solve(t(A) %*% A) %*% t(A)
#---------------------------------------------
# 中央の残差 eps(0) を x に対する線形結合で表す。
#
# y = B x
# a = (A^T A)^{-1} A^T y = ATA_inv_AT %*% y
# y[0] = c^T y
# q(0) = a_0 = e^T a = e^T (A^T A)^{-1} A^T y
#
# eps(0) = y[0] - q(0)
# = (c^T - e^T (A^T A)^{-1} A^T) y
# = (c^T - e^T ATA_inv_AT) B x
#
# 従って、w^T = (c^T - e^T ATA_inv_AT) B
#---------------------------------------------
w <- as.numeric(c_vec %*% B - e_vec %*% ATA_inv_AT %*% B)
list(
k = j, # インデックス -M,...,M
w = w # FIR係数（元の時系列 x[n+j] への重み）
)
}
############################################################
# 2. 左端用 FIR 設計関数
#
# 区間 [1, n0+M] のみが利用可能な場合
# （右側には M 点確保できるが左側が足りない端点）に対して、
# そこだけで累積和 → p 次多項式フィット → 残差を求める
# FIR を計算する。
#
# eps(n0) = c[n0] - q(n0) を
# eps(n0) = sum_k w_k x[k]
# の形に落とす。
############################################################
sg_dma_left_edge_fir <- function(M, p, n0) {
#---------------------------------------------
# 引数チェック
#---------------------------------------------
if (M < 1) stop("M must be >= 1.")
if (p < 0) stop("p must be >= 0.")
if (n0 < 1) stop("n0 must be >= 1.")
if (n0 >= M + 1) warning("この関数は n0 < M+1 の左端用として想定しています。")
# 利用する区間は [1, T]（T = n0 + M）
# 例: n0=3, M=5 のとき [1,...,8] を使って n0=3 の残差を出す。
T <- n0 + M
t_vec <- 1:T
#---------------------------------------------
# 累積和行列 S (T x T)
# c = S x とすると、c[t] = sum_{r <= t} x[r]
#---------------------------------------------
S <- matrix(0, nrow = T, ncol = T)
S[lower.tri(S, diag = TRUE)] <- 1
#---------------------------------------------
# 多項式デザイン行列 A (T x (p+1))
# A_{t,m} = t^m
# a = (a_0,...,a_p)^T に対して q(t) = A %*% a
#---------------------------------------------
A <- outer(t_vec, 0:p, function(tt, mm) tt^mm)
#---------------------------------------------
# c[n0] を取り出すベクトル e_{n0}
# e_{n0}^T c = c[n0]
#---------------------------------------------
e_n0 <- rep(0, T)
e_n0[n0] <- 1
#---------------------------------------------
# q(n0) = r0 a, ここで r0 は A の n0 行
# r0 %*% a = q(n0)
#---------------------------------------------
r0 <- A[n0, ]
#---------------------------------------------
# a = (A^T A)^{-1} A^T c = ATA_inv_AT %*% c
#---------------------------------------------
ATA_inv_AT <- solve(t(A) %*% A) %*% t(A)
#---------------------------------------------
# eps(n0) = c[n0] - q(n0)
# = e_n0^T c - r0 a
# = e_n0^T c - r0 (A^T A)^{-1} A^T c
# = (e_n0^T - r0 ATA_inv_AT) c
# = (e_n0^T - r0 ATA_inv_AT) S x
#
# よって w^T = (e_n0^T - r0 ATA_inv_AT) S
#---------------------------------------------
w <- as.numeric(e_n0 %*% S - r0 %*% ATA_inv_AT %*% S)
list(
k = t_vec, # インデックス 1,...,T
w = w # FIR係数（x[1],...,x[T] への重み）
)
}
############################################################
# 3. SG-DMA 本体
#
# sg_dma_sg_fast(x, M, p)
# - x : 時系列ベクトル
# - M : 半窓幅（整数, 窓長 s=2M+1）
# - p : 多項式次数
#
# 処理の流れ：
# 1) x の平均を引く
# 2) プロファイル (累積和) 上での SG-DMA 前処理に対応する
# FIR（sg_dma_local_fir）を計算
# 3) 中央領域 (n = M+1,...,N-M) は filter() で一括計算
# 4) N <= 1024 のときだけ、左端・右端も精密に処理
# 5) eps(n)^2 の平均値を F^2(s)、その平方根を F(s) として返す
############################################################
sg_dma_sg_fast <- function(x, M, p) {
x <- as.numeric(x)
N <- length(x)
#---------------------------------------------
# 引数チェック
#---------------------------------------------
if (N < 2) stop("x の長さが短すぎます。")
if (M < 1) stop("M must be >= 1.")
if (p < 0) stop("p must be >= 0.")
# 実務的には N >= 2*M+1 を推奨（窓をきちんと取れるように）
if (N < 2 * M + 1) {
warning("N < 2*M+1 です。中央部の有効な窓がほとんどありません。")
}
#---------------------------------------------
# 1. 平均を引く（オフセット除去）
# SG-DMA の前処理として標準的な操作。
# この後の操作は線形なので、平均を引いても引かなくても
# DC 成分には感度を持たないが、定義上は引いておく。
#---------------------------------------------
x0 <- x - mean(x, na.rm = TRUE)
#---------------------------------------------
# 2. ローカル窓 [-M, M] 用 FIR（中央領域）
# sg_dma_local_fir は 1 ウィンドウ [-M,...,M] に対する
# eps(0) の FIR 係数を返す。
# これを畳み込みとして全 n に適用する。
#---------------------------------------------
local_fir <- sg_dma_local_fir(M = M, p = p)
w_local <- local_fir$w # 長さ 2M+1
# 残差 eps(n) を格納するベクトル（端点は NA の可能性あり）
eps <- rep(NA_real_, N)
#---------------------------------------------
# 3. 中央部: n = M+1,...,N-M
# - 左右に M 点ずつデータがある領域。
# - stats::filter() を使うことで for ループなしに
# 全点を一括計算する。
# - sides=2 を指定することで、フィルタ中心をなかほどに置く。
#---------------------------------------------
if (N >= 2 * M + 1) {
eps_mid <- as.numeric(stats::filter(x0, w_local, sides = 2))
eps[(M + 1):(N - M)] <- eps_mid[(M + 1):(N - M)]
}
#---------------------------------------------
# 4. 端点処理：
# N が十分大きいときは計算コスト優先で端点無視でも影響が小さい。
# ここでは N <= 1024 の場合のみ、左端・右端を
# sg_dma_left_edge_fir を用いて「きちんと」補正する。
#---------------------------------------------
if (N <= 1024) {
# 左端で処理する点の数（M もしくは N の小さい方まで）
max_left <- min(M, N)
if (max_left >= 1) {
# 左端 n0 = 1,...,max_left 用の FIR を事前にすべて計算
left_firs <- lapply(1:max_left, function(n0) {
sg_dma_left_edge_fir(M = M, p = p, n0 = n0)
})
# 左端: n0 = 1,...,max_left に対して eps(n0) を計算
for (n0 in 1:max_left) {
firL <- left_firs[[n0]]
T_use <- min(length(firL$w), N) # 念のため N でクリップ
eps[n0] <- sum(firL$w[1:T_use] * x0[1:T_use])
}
#-----------------------------------------
# 右端: n = N-M+1,...,N
# 反転系列 x_rev に対する「左端問題」として同じ FIR を再利用。
# - x_rev[k] = x0[N - k + 1]
# - n に対応する反転位置は n0_rev = N - n + 1
#-----------------------------------------
x0_rev <- rev(x0)
max_left_rev <- max_left
for (n in max(N - M + 1, 1):N) {
n0_rev <- N - n + 1 # 反転系列での位置
if (n0_rev >= 1 && n0_rev <= max_left_rev) {
firL_rev <- left_firs[[n0_rev]] # 左端 FIR を再利用
T_use <- min(length(firL_rev$w), N)
eps_rev <- sum(firL_rev$w[1:T_use] * x0_rev[1:T_use])
eps[n] <- eps_rev
}
}
}
}
#---------------------------------------------
# 5. F^2(s), F(s) の計算
# - eps(n) のうち NA でない点について eps(n)^2 の平均をとる。
# - これが DMA の F^2(s) に対応。
# - その平方根を F(s) として返す。
#---------------------------------------------
F2 <- mean(eps^2, na.rm = TRUE)
F <- sqrt(F2)
list(
eps = eps, # 各 n における残差（端点は NA の場合あり）
F2 = F2, # DMA fluctuation (二乗平均)
F = F # その平方根
)
}
■ RでDMA
さらに、サンプル時系列を生成してDMAを実際にやってみるＲスクリプトを掲載しておきます。事前に、上のRスクリプトを一度実行するのを忘れないでください。
############################################################
# Detrending Moving Average Analysis (DMA)で長時間相関時系列を解析する例
#
# - longmemo::simFGN0() を用いて Hurst 指数 H の
# Fractional Gaussian Noise (FGN) を生成する。
# - FGN は 1/f^α 型スペクトルをもつ長期相関ノイズであり、
# DMA（Detrending Moving Average）で H を推定する典型的な教材。
#
# 【パッケージの読み込みについて】
# longmemo パッケージには “simFGN0()” が含まれており、
# 時系列長 N と Hurst 指数 H を指定するだけで
# FGN（自己相似過程の差分系列）を簡単に生成できる。
#
# 【longmemo のインストール方法】
# longmemo は CRAN から通常通りインストールできる：
# install.packages("longmemo")
#
# ここでは FGN の生成に利用するため、
# 必ず longmemo を読み込んでおく必要がある。
############################################################
require(longmemo)
############################################################
# シミュレーション設定
H <- 0.8 # 真の Hurst 指数
N <- 1000 # 時系列長
x <- simFGN0(N, H)
# Savitzky–Golay の多項式次数
# p = 0: 移動平均に近い（定数トレンドのみ）
# p = 2: 2次多項式トレンドまで除去
p <- 0
# スケール（半窓幅 M）の系列を対数スケールで生成
M <- unique(round(exp(seq(
log(ceiling(p/2) + 1), # 最小 M（p に応じて少し大きめから）
log(round(length(x) / 10)), # 最大 M（N/10 程度まで）
length.out = 20 # スケールの点数
))))
s <- 2 * M + 1 # 実際の窓長 s
# 各スケールで F(s) を計算
get_F_for_M <- function(Mi) {
res <- sg_dma_sg_fast(x, Mi, p)
res$F
}
F_vec <- vapply(M, get_F_for_M, numeric(1L))
log10.F.s <- log10(F_vec)
log10.s <- log10(s)
# log10(s) vs log10 F(s) の線形回帰
plot(log10.s, log10.F.s,
xlab = "log10(s)",
ylab = "log10 F(s)",
main = "SG-DMA fluctuation vs scale")
lfit <- lm(log10.F.s ~ log10.s)
abline(lfit, col = 2, lty = 2)
# 傾きを図中に表示
slope <- coef(lfit)[2]
# 描画領域の座標を取得して、右上に表示
usr <- par("usr") # c(xmin, xmax, ymin, ymax)
text_x <- usr[1] + 0.05 * (usr[2] - usr[1])
text_y <- usr[4] - 0.05 * (usr[4] - usr[3])
text(text_x, text_y,
labels = sprintf("slope = %.3f", slope),
adj = c(0, 1)) # 左寄せ・上揃え
# 回帰の結果（必要ならコンソールにも表示）
lfit
■ インパルス応答と周波数応答の計算・プロット例
# フィルタ長 s = 2 * M + 1
M <- 50
# フィルタ次数
p <- 2
# インパルス応答の計算
FIR.DMA <- sg_dma_local_fir(M, p)
par(mfrow = c(1, 2), cex.lab = 1.5)
## 1. インパルス応答 w_k ------------------------------------------
plot(
FIR.DMA$k, FIR.DMA$w,
col = 4, pch = 16,
xlab = "k",
ylab = expression(w[k]),
main = sprintf("FIR length = %d, order = %d", length(FIR.DMA$w), p)
)
## 2. 周波数応答 ---------------------------------------------------
k <- FIR.DMA$k
w <- FIR.DMA$w
# 周波数軸（0は禁止なので少しずらす）
Nf <- 2000
f_min <- 1e-4
f <- seq(f_min, 0.5, length.out = Nf)
# H(f) = Σ w_k exp(-i2πfk)
E <- exp(-1i * 2 * pi * outer(k, f))
H <- as.vector(t(w) %*% E)
H_mag <- Mod(H)
# ---- ここで log10 の範囲を自動決定 ----
H_min_log <- floor(log10(min(H_mag)))
H_max_log <- ceiling(log10(max(H_mag)))
major_yticks <- 10^(H_min_log:H_max_log)
minor_yticks <- outer(major_yticks, 1:9)
# プロット（軸は後で描く）
plot(
f, H_mag,
log = "xy",
type = "l",
xlab = "f (cycles/sample)",
ylab = expression("|H(f)|"),
col = 4, lwd = 3,
xaxt = "n", yaxt = "n",
main = "Log–log frequency response"
)
###########################################
# 対数軸の描画：X軸（固定でOK）
###########################################
major_xticks <- 10^(-4:0)
minor_xticks <- outer(major_xticks, 1:9)
axis(1, at = minor_xticks, labels = FALSE, tck = -0.01)
axis(
1, at = major_xticks,
labels = parse(text = paste0("10^", -4:0)),
las = 1, cex.axis = 1.2, tck = -0.02
)
###########################################
# 対数軸の描画：Y軸（自動スケール）
###########################################
# 補助目盛
axis(2, at = minor_yticks, labels = FALSE, tck = -0.01)
# 主目盛
axis(
2, at = major_yticks,
labels = parse(text = paste0("10^", H_min_log:H_max_log)),
las = 1, cex.axis = 1.2, tck = -0.02
)
5. まとめ
本稿では、DMA（Detrending Moving Average Analysis）に含まれる一連の処理が、すべて 線形演算（FIR フィルタ）として書き表せることを示しました。
DMA は、見た目こそ「積分して、多項式を当てて、引き算をする」という複雑な処理に見えますが、本質的には 1 本の畳み込みフィルタの応答を求めているだけです。この視点に立てば、DMA はデジタル信号処理の基本原理と完全に整合し、周波数特性・安定性・高速化などの議論も統一的に理解できます。
それでも現実には、多くの研究者が有名なDFA を使い続けます。歴史的な定着や「権威」の力は強く、数学的に洗練されたDMA が正当に評価されるには、もう少し時間がかかるのかもしれません。
※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。それと、下のはてなブログランキングはあまり信用できる指標ではなさそうですが (私のブログを読んでいる人は、実際とても少ないです)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。
