統計数理（2023）
第 71 巻 第 1 号 65–80©2023 統計数理研究所

特集「諸科学における統計数理モデリングの拡がり」［総合報告］

トレンド除去演算に基づく
フラクタル時系列解析の数理——非定常時系列にみられる 1/f β 型ゆらぎ——清野 健†（受付 2022 年 6 月 30 日；改訂 10 月 16 日；採択 11 月 7 日）

要旨

非定常なトレンド成分を含む時系列にみられる長時間相関，1/f β 型ゆらぎ，あるいは，自己アフィンフラクタル性などの特性を評価するために，detrended fluctuation analysis（DFA）と呼ばれる解析法が，物理学や生体信号解析の分野で頻繁に用いられるようになっている．DFA では，非整数ブラウン運動のような拡散的時系列の統計的自己アフィン性（Hurst 指数）を定量化する手続きの中に，時系列のトレンド除去演算が含まれている．そのようなトレンド除去演算の導入には，非定常なトレンド成分に起因する誤った Hurst 指数の推定を避けたり，評価可能なスケーリング指数の範囲を拡張できたりする利点がある．最近では，R，Python のパッケージとして DFA が実装され，誰もが DFA を簡単に利用できるようになったため，パワースペク
トルの推定に代わり DFA が用いられることが多くなっている．本稿では，DFA とその派生版（DMA）を中心として，これらの解析法の数理的基である detrending moving average algorithm
礎を解説する．
キーワード：長時間相関，フラクタル，長期記憶，1/f ゆらぎ，時系列解析．
1.

はじめに

心拍変動，呼吸変動，脳波など，多くの生体信号時系列は，周波数を f として 1/f β 型パワー
．β = 1 に近いもの，ある
スペクトルを示す（Stanley et al., 1994; Ward and Greenwood, 2010）
いは，0 < β < 2 の範囲のものは，1/f ゆらぎと呼ばれることがある．ここでは，β の範囲を制限せず，1/f β 型パワースペクトルを示す時系列を，1/f β 型ゆらぎと呼ぶ．1/f β 型ゆらぎは，非整数ガウスノイズに見られる長時間相関，あるいは，非整数ブラウン運動に見られる自己ア
フィンフラクタル性と関連している（Mandelbrot and Van Ness, 1968）．弱定常性を仮定できる時系列の場合は，パワースペクトル密度の推定に基づき，長時間相関を評価できる．しかし，生体信号時系列をはじめとして，自然界に見られる時系列のほとんどは非定常であるため，
1/f β 型ゆらぎのスケーリング指数 β を推定するには，時系列に含まれる非定常トレンドの影
響を除く必要がある．そのような非定常トレンドに対応するために，最近の 20 年で detrended（DFA）と呼ばれる解析法が，物理学や生体信号解析の分野で頻繁に用いら
fluctuation analysis
れるようになっている（Peng et al., 1995）．DFA では，非整数ブラウン運動のような拡散的時†

大阪大学大学院 基礎工学研究科：〒 560–8531 大阪府豊中市待兼山町 1 番 3 号

統計数理

66

第 71 巻 第 1 号 2023

系列の Hurst 指数 H を推定する手続きの中に，時系列のトレンド除去演算が含まれている．そのため，トレンド成分に起因する誤った Hurst 指数の推定を避けることができる利点がある．
最近では，R，Python のパッケージとして，DFA を簡単に利用できるため，古典的なパワースペクトルの推定に代わり DFA が用いられることが多くなっている．本稿では，まず，1/f β
型ゆらぎの特性について紹介し，その後，DFA とその派生版である detrending moving average（DMA）（Carbone and Kiyono, 2016）の数理的基礎について解説する．これまで，こalgorithm
れらの解析法については，主に数値実験に基づきその特性が議論されてきた（Hu et al., 2001;．それに対し，本稿では解析的にこれらの解析法の特性を議Chen et al., 2002; Ma et al., 2010）
論する．
2.

1/f β 型ゆらぎの特性

本稿において注目する時系列の特性は，長期記憶，長時間相関，フラクタル性などと呼ばれるものである（図 1）．長期記憶，長時間相関，フラクタル性は，同じ意味のように使われることがあるが，異なる特徴を表す場合があることに注意してほしい．
長期記憶は，弱定常性の仮定の下で，離散時系列の自己共分散関数 C[k] が，
(2.1)

∞


|C(k)| = ∞k=0を満たし，C(k) ∼ k−γ (0 < γ < 1) の形をとる（Beran, 2017; Li, 2010）．ここで，∼ は比例を表す．式（2.1）の性質により，自己共分散関数は特徴的時間スケールを持たないため，その特性はスケーリング指数 γ で特徴付けられる．

図 1．ガウス過程時系列とパワースペクトル密度 S(f ) の関係．サンプル時系列の右側にそのパワースペクトルの推定結果（100 例の平均）を描いた．左側が非整数ガウスノイズの例
であり，右側が非整数ブラウン運動の例．

トレンド除去演算に基づくフラクタル時系列解析の数理

67

長時間相関の概念は，長期記憶を含むが，それ以外に長時間反相関と呼ばれる性質も含んで．長時間反相関の場合，自己共分散関数は負の値をとり，自いる（Rangarajan and Ding, 2000）
己共分散関数の絶対値はべき的減衰（1 < γ < 2）を示すが，式（2.1）を満たすほどゆっくりではない．長時間相関を示す例として，非整数ガウスノイズがある．非整数ガウスノイズは，その自己共分散関数が，0 < H < 1 をパラメタとして，
1(2.2)
C(k) = (|k + 1|2H − 2 |k|2H + |k − 1|2H )2で定義されるガウス過程である（Mandelbrot and Van Ness, 1968）．|k| が大きくなると，漸近
的に，
(2.3)

|C(k)| ∼ |k|2(H−1)となる．0.5 < H < 1 の場合が長時間相関（=長期記憶過程），0 < H < 0.5 の場合が長時間反相関の例となる．また，非整数ガウスノイズは，Hurst 指数 H の非整数ブラウン運動の差分過程になっている．
また，時系列のフラクタル性は，非整数ブラウン運動のような統計的な自己アフィン性をもつ拡散的軌道を意味し，ハースト指数 H とフラクタル次元 D の間に D = 2 − H の関係が成り立つ（Moreira et al., 1994）．非整数ブラウン運動については，その増分が長時間相関を持つ過程になっている．長期記憶，長時間相関を持つ時系列自体にフラクタル性はないが，それらの時系列の積分
（累積）
は．フラクタル性を持つようになる．
以上の点を踏まえると，長期記憶，長時間相関，フラクタル性の名称は，統一的な広い枠組みを与えるには無理がある．そこで，ここでは，それらの特性を含む枠組の呼称として，1/f β型ゆらぎを使うことにする．ここでは，1/f β 型ゆらぎを，時系列に含まれる確率的変動について，有限長時系列のパワースペクトル S(f ) が，S(f ) ∼ f −β を示すものとして定義する．ここで，確率的変動とは，その変動成分の n 階差分が弱定常過程になる実数値 n が存在することを意味する．n 階差分の分散が有限であれば，確率的変動が従う確率分布についての制限はな
い．したがって，本稿の議論は，有限な分散をもつ非ガウス過程，および，n 階差分が弱定常となる非ガウス過程に対して適用可能である．
弱定常過程の場合，長期記憶は 0 < β < 1 の領域に対応し，γ + β = 1 が成り立つ．また，長時間相関は −1 < β < 0 または 0 < β < 1 の領域に対応し，0 < β < 1 の場合を長時間正相関，−1 < β < 0 の場合を長時間反相関と呼ぶ．β = 0 は，無相関な白色ノイズに対応する．自然界には β = 1 に近いゆらぎがしばしば観測され，そのようなゆらぎは 1/f ゆらぎと呼ばれる．
1/f ゆらぎの起源については，不明なことが多い．フラクタル時系列は，1 < β < 3 の領域に対応し，ハースト指数 H の非整数ブラウン運動と (β − 1)/2 = H の関係が成り立つ．この場合，時系列は非定常な拡散過程となる．つまり，β ≥ 1 の場合，数学的には時系列の分散が発散す
るため，そのような過程のパワースペクトル密度が存在するとはいえない．そのような場合でも，Wigner-Ville スペクトルを代用して β を定義することができる（Flandrin, 1989）．非定常過程を含んで，スケーリング指数 β は多くの場合，意味のあるパラメタである．
β はパワースペクトル解析におけるスケーリング指数として推定されるものである．実用的なパワースペクトルの推定では，観測時系列を複数個の部分区間に分割してパワースペクトルの平均を求めたり，窓関数を適用したり，パワースペクトルの平滑化を行ったりと，様々なテク．
ニックが必要となり，このことが非専門家に扱いづらい印象を与えるようである（日野, 1996）
そのような煩わしさを軽減しながら，非定常時系列の解析にも適用可能な方法として登場した（DMA）
解析法が，detrended fluctuation analysis（DFA）や detrending moving average algorithm
である（Shao et al., 2012）．次節ではまず，これらの解析法の基礎となるランダムウォーク解

統計数理

68

第 71 巻 第 1 号 2023

析について紹介する（Peng et al., 1992）．
3.

ランダムウォーク解析

ここでは，γ と β に加えて，ランダムウォーク解析に基づく新たなスケーリング指数 α を導入する．ランダムウォーク解析では，まず，弱定常過程の観測時系列 {x[i]} について，事前に標本平均を 0 にして積分時系列を求める：
y[i] =

(3.1)

i


x[j].

j=1

このように観測時系列をランダムウォーク軌道に変換することで，非整数ブラウン運動のような拡散的な時系列を生成する．
ランダムウォーク解析では，積分時系列の増分(3.2)

Δs y[i] = y[i + s] − y[i] =s

x[i + j]

j=1の 2 乗平均の期待値（ゆらぎ関数:fluctuation function）(3.3)

F 2 (s) = Δs y 2 [i]を考える（図 2）．ここで， ·  は期待値を表す（本稿では，確率変数も小文字で表す）．1/f β 型
ゆらぎでは，
F 2 (s) ∼ s2α(3.4)となり，スケーリング指数 α と β には
(3.5)

α=

β+1
2

(0 < α < 1)

図 2．非整数ガウスノイズのランダムウォーク解析．（a）H = 0.8 のサンプル時系列．（b（）a）の積分時系列．
（c）ランダムウォーク解析の結果（100 例の平均）．両対数プロットの傾（100 例のきからスケーリング指数 α が推定される．（d）パワースペクトルの推定結果
平均）．両対数プロットの傾きからスケーリング指数 β が推定される．β = 2α − 1 が成り立つ．

トレンド除去演算に基づくフラクタル時系列解析の数理

69の関係が成り立つ．以下では，このようなスケーリング指数の関係式の導出について説明する．
弱定常過程のランダムウォーク解析において，F 2 (s) と C(k) の関係は
F 2 (s) =

(3.6)

s−1


(s − |k|)C(k)k=−(s−1)で与えられる．したがって，C(k) ∼ |k|−γ を仮定することで，s  1 のとき，
F 2 (s) ∼ s2−γ

(3.7)

となる．よって，α = 1 − γ/2 が成り立つ（1/2 < α < 1）．
2
また，F (s) と S(f ) の関係は，
 1/2
sin2 (πf s)(3.8)
df
F 2 (s) =S(f )sin2 (πf )−1/2となる．式（3.6）にある三角窓関数（テント関数，三角形関数）
は，2 つの矩形窓関数の畳み込みには，矩形窓関数のフーリエ変換である sinc 関数の 2 乗として表すことができるので，式（3.8）
が入る．
α = (β + 1)/2 の関係式は，式
（3.8）
において，
⎧
1
⎪
⎨ s2
for 0 < |f | <sin2 (πf s)
πs
≈
(3.9)
1
1
⎪
sin2 (πf )⎩ (πf )−2
for
< |f | ≤πs21と近似することで，導くことができる．式（3.9）は，0 < |f | < πsの領域において，f のマク
1< |f | ≤ 12 の領域では，s  1 のとき，分子ローリン展開の 1 次の項までを残した．また， πsが分母に比べて非常に早く振動するので，分子にある振動の上側のエンベロープのみをとり（1と近似），展開した．
（3.8）は，
時系列のパワースペクトルが，S(f ) ∼ f −β に従うとき，式 1/(πs) 1/2(3.10)
F 2 (s) ∼f −β s2 df +f −β (πf )−2 df
01/(πs)と，近似的に評価できる．ここで，−1 < β < 1 のとき，
 1/(πs)f −β s2 df ∼ sβ+1
 1/2

0

1/(πs)22α

f −β (πf )−2 df ∼ sβ+1であるので，F (s) ∼ s と比較して，α = (β + 1)/2 となる．また，β < −1 では，|s|  1 のとき，F 2 (s) が一定値に漸近するため，α = 0 となる．β > 1 では，式（3.10）の第一項の積分において，f = 0 に近い領域の寄与が支配的になるため，F 2 (s) ∼ s2 ，すなわち，α = 1 となる．
したがって，ランダムウォーク解析では，0 < α < 1 (−1 < β < 1) の領域のみ，β との関係式のように，F 2 (s) と S(f ) の関係式を導くことで，スケーリング指数の検が成り立つ．式（3.8）
出範囲など，解析手法の基本特性を明らかにすることができる（Kiyono, 2017）．
（FA）
これまで，ランダムウォーク解析は，乱流統計の構造関数解析や，fluctuation analysis
として使われてきた（Peng et al., 1992; Shao et al., 2012）．ここでは，ランダムウォーク解析の明確な欠点を 2 つ指摘する．一つは，評価可能なスケーリング指数 α の領域がスペクトル解

70

統計数理

第 71 巻 第 1 号 2023

析に比べて狭いことである．つまり，ランダムウォーク解析では，−1 < β < 1 に対応する領域
しか分析できない．もう一つは，トレンド成分を含む非定常時系列では α の推定精度が低下することである．これらの欠点を解消する方法が，トレンド除去演算の導入である．次節では，トレンド除去演算を含む DFA および DMA を導入する．
4.

トレンド除去演算を導入したランダムウォーク解析

ここではまず，時系列解析における DFA と DMA の手順を説明する（図 3）．その後，DFAと DMA の特性をまとめる（Kiyono, 2015; Kiyono and Tsujimoto, 2016b; Carbone and Kiyono,
2016; Höll et al., 2019）．DFA も，DMA も，ランダムウォーク解析にトレンド除去演算が含まれる点のみ異なる．
4.1 DFA と DMA の手順
DFA の手順は以下である（Peng et al., 1995）．

（1）平均 0 の時系列 {x[i]} の積分時系列 y[i] = ij=1 x[j] を計算する．
{y[i]} を長さ（スケール）s の部分区間に分割し，各部分区間の変動に対して m 次多項式の（2）．以下では，l 番目の部分区間の回帰曲線を
回帰曲線を最小 2 乗法であてはめる（図 3 左）p(l) [i] とする．スケール s は，m + 1 より大きく，N/10 程度より小さい範囲に設定する．
y[i] と p(l) [i] の差の 2 乗平均の平方根（fluctuation function: ゆらぎ関数）（3）

図 3．DFA（左列）と DMA（右列）における積分時系列に対するトレンド除去演算（s = 101とした）．（a）DFA では，部分区間に回帰直線（あるいは曲線）をフィットする（破線）．
（破線は 0（b）DMA では，m 次 Savitzky-Golay フィルタを用いてトレンドを推定する
次の例）．（c,d）トレンド除去後の積分時系列．（e）DFA では非線形フィルタを用いるため，ガウス過程時系列を分析した場合も，トレンドまわりの分布は非ガウスになる（上に凸な破線がガウス分布）．（f）DMA では線形フィルタを用いるため，ガウス過程の分
布特性は維持される（上に凸な破線がガウス分布）．

トレンド除去演算に基づくフラクタル時系列解析の数理

⎡
(4.1)

1
F (s) = ⎣N/s⎧

N/s ⎨


l=1

1
⎩s

ls


(l)s

(y[k] − pk=(l−1)s+1

71

⎫⎤1/2
⎬
[k]) ⎦ ,⎭2を計算する．ここで， · は，床関数である．
s の値を変化させ，F (s) を計算する．
（4）
log10 s に対して log10 F (s) をプロットし，直線（スケーリング）領域の傾きをスケーリン（5）
グ指数 α として推定する．
m 次多項式をトレンド除去に用いる方法は，m 次 DFA と呼ばれる．DFA では，検出可能なスケーリング指数 α の範囲が，0 < α < m + 1 となる．さらに，解析する時系列 {x[i]} に含まれる m − 1 次の多項式トレンドを完全に除くことができる（Hu et al., 2001）．
DMA では，DFA の区分的曲線フィットの代わりに，Savitzky-Golay フィルタを用いる（Tsujimoto et al., 2016）．Savitzky-Golay フィルタは，移動平均の定数フィットを多項式に一般化したフィルタである．m 次 Savitzky-Golay フィルタでは，スケール （s 奇数）の部分区間にm 次多項式を当てはめ，その中央の点をフィルタリングされた値として採用する．部分区間
は，一点ずつスライドさせる．なお，m は偶数とする（m 次と m + 1 次のフィルタリング結果（central DMA）では，移動平均フィルタが用いられは一致するため）．初期に提案された DMA）3）のており，これは 0 次の Savitzky-Golay フィルタと一致する．DMA の手順は，DFA の（2（部分のみ異なるのでその箇所のみ説明を与える．
{y[i]} に対し，スケール s，次数 m の Savitzky-Golay フィルタを適用する（図 3 右）．その（2）
結果を ps [i] とする．スケール s は，m + 3 より大きく，N/10 程度より小さい範囲の奇数に設定する．
y[i] と ps [i] の差の 2 乗平均の平方根
（3）
⎤1/2
⎡
N −(s−1)/2

1
(4.2)
F (s) = ⎣(y[k] − ps [k])2 ⎦ ,N −s+1k=(s+1)/2を計算する．
DMA の手順は，DFA と対比することで，説明としてはわかりやすいが，実用上は，DMA の高速アルゴリズムを使用することを薦める（Tsujimoto et al., 2016）．上記の手順では，計算時間が時系列の長さ N 2 に比例するが，高速アルゴリズムでは N に比例した計算速度が実現されている．
m 次 Savitzky-Golay フィルタをトレンド除去に用いる方法は，m 次 DMA と呼ばれる．DMA
では，検出可能なスケーリング指数 α の範囲が，0 < α < m + 2 となる．さらに，解析する時
系列 {x[i]} に含まれる m 次の多項式トレンドを完全に除くことができる．これらの特性は解析的に示されている（Carbone and Kiyono, 2016）．
DFA も，DMA も，ランダムウォーク解析と類似の方法のように見える．DFA および DMA
では，トレンド除去の方法がスケールに依存しているので，ランダムウォーク解析と同じス
ケーリング指数 α が推定できると予想できる．次節では，この点の明確にするために，DFA とDMA の数理的基礎をまとめる．
4.2 自己共分散関数 C(k) およびパワースペクトル S(f ) との関係
弱定常の仮定の下で，DFA と DMA の F 2 (s) を，自己共分散関数 C(k) およびパワースペク
トル密度 S(f ) を用いて以下の形で表すことができる（Kiyono, 2017）．

72

統計数理

第 71 巻 第 1 号 2023

(4.3)

F 2 (s) =sC(k) L(k, s)

k=−s

F 2 (s) =

(4.4)

 1/2
−1/2
s


|Gs (f )|2 =(4.5)

S(f ) |Gs (f )|2 df

L(k, s) cos(2πf k)k=−s以下では，これらの式に含まれる L(k, s) と |Gs (f )|2 の求め方について説明する．
ま ず は ，式 表 現 が 簡 単 な DMA を 考 え る ．ラ ン ダ ム ウ ォ ー ク 解 析 に お け る Δs y[i] =s
j=1 x[i + j] を，重み {ws (j)} を導入して，
Δs y[i] =

(4.6)

s


ws (j) x[i + j]j=1と一般化し，
F 2 (s) = Δs y 2 [i]

(4.7)

とする．DMA では，積分演算も，Savitzky-Golay フィルタを用いたトレンド除去演算も，畳の {ws (j)} を求めることができる．具体的には，DMA のみ込みとしてできるので，式（4.6）{ws (j)} は，
⎧

 
m
s

⎪
⎨ Θ s+1 −k −b1,j+1(r − s − 1)j for 1 ≤ k ≤ s
2ws (k) =
(4.8)
j=0
r=k
⎪
⎩
0
for otherwiseとなる（Kiyono and Tsujimoto, 2016b）．ここで，Θ(x) はユニットステップ関数
1(x ≥ 0)(4.9)
Θ(x) =
0(x < 0)であり，bi,j (s) は，行列(4.10)⎡

1⎢⎢0(s−1)/2 ⎢⎢ j2
Bm (s) =⎢⎢ .
j=−(s−1)/2 ⎢
.
⎣ .
jm

0
j2
0
..
.
0

j2
0
j4
..
.

···
···
···
..
.

j m+2

···

jm
0

⎤

⎥
⎥
⎥
m+2 ⎥
j
⎥
⎥
..
⎥
.
⎦
j 2m

の逆行列の (i, j) 成分である．図 4 は，0 次，2 次，4 次 DMA の {ws (k)} を描いたものである．
式（4.6）の表現を用いると，
s−|k|(4.11)L(k, s) =

ws (j) ws (j + |k|)j=1

(4.12)Gs (f ) =sws (k)e−i2πf kk=1となる．具体例として，0 次 DMA の場合は，

トレンド除去演算に基づくフラクタル時系列解析の数理

73

図 4．DMA の重み係数 {ws (k)}．m は，トレンド除去に用いる Savitzky-Golay フィルタの次数．

(4.13)⎧
2k3 + 6sk2 − 2(3s2 + 1)k + s3 − s⎪⎨
12s2L(k, s) =⎪⎩ (k − s)(k − s + 1)(k − s − 1)
6s2

(4.14)|Gs (f )| =1
2 sin(πf )

s−1
2
s−1
for |k| >2
for |k| ≤sin(πsf )1−
s sin(πf )となる．これらの関数の漸近形を考えることで，スケーリング指数 α，β ，γ の間の関係式を解析的に導くことができる．加えて，α の推定上限を |Gs (f )|2 から求めることができる（Kiyono．
and Tsujimoto, 2016b）
DFA の場合，トレンド除去演算を行う長さ s の部分区間全体の期待値を考える必要があるため，重み {ws (i, k)} を導入して(4.15)

Δs y[i] =sws (i, k) x[i + k]k=1を定義すれば，(4.16)

2

F (s) =1Δs y 2 [i + j]
s j=1
s



となる．DMA の場合，式（4.7）のように，F 2 (s) の期待値は各点 i で評価されている．それに（4.16）のように，長さ s の部分区間 [i + 1, i + s] の期待値が評価されて対し．DFA の場合，式いる．DFA では，この部分区間内で，各点が F 2 (s) に与える寄与が非一様になっているため，部分区間全体で平均値をとる必要がある（Kiyono and Tsujimoto, 2016a）．また，DMA では，
Δs y[i] が，式（4.6）のように，畳み込みで表現されているため，ガウス過程時系列の Δs y[i] は，
ガウス分布に従う．それに対し，DFA では，部分区間内で，Δs y[i] の重みが非一様であるため，ガウス過程時系列の Δs y[i] は，非ガウス分布になる（図 3）．DMA と比べた場合に，これらの特性は DFA の推定量としての欠点である．
DFA の場合，
s−|k|(4.17)L(k, s) =1  
ws (i, j) ws (i, j + |k|)s j=1 i=1sとなる．|Gs (f )|2 は，式（4.5）
を用いて計算できる．しかし，DFA では，式（4.11）および式（4.12）
を，使用することができず，計算の手間は DMA とくらべて非常に大きい．

統計数理

74

第 71 巻 第 1 号 2023

4.3 白色ノイズの解析結果
前節で，定義した L(k, s) を用いると，白色ノイズを DMA あるいは DFA で分析した結果を解析的に求めることができる．白色ノイズの自己共分散関数 C(k) は，時系列の分散を σ 2 として，
σ2for k = 0
(4.18)
C(k) =
0for k = 0である．式（4.3）
より，(4.19)

F 2 (s) =sC(k) L(k, s) = σ 2 L(0, s)k=−sとなる．ランダムウォーク解析（FA），1 次 DFA（DFA1），0 次 DMA（DMA0）では，それぞれ，
以下のようになる．
(4.20)2
FFA(s) = σ 2 s(4.21)

2(s) = σ 2
FDMA0

s2 − 112s(4.22)

2(s) = σ 2FDFA1s2 − 415s

ランダムウォーク解析では，s の全領域において，FFA (s) が s1/2 に比例するが，DMA および
DFA では，短いスケールでずれが生じることが分かる（図 5 参照）．このような解析的結果は，実際の時系列解析の結果を解釈する際に役に立つ．
4.4 スケーリング指数 α と β の関係（4.5）
を用いることで，スケーリング指数 α，β ，γ の間の関係を解析的に議式（4.3）および式論することができる．例えば m 次 DMA の場合，図 6 に示したように，|Gs (f )|2 の低周波数側は s2m+4 f 2m+2 に比例し，高周波数側は f −2 に比例する．0 次 DMA の場合，
⎧ √2π(s2 − 1)
6
⎪
⎪
⎨
f2
for 0 < f < fc =2
12πs
(4.23)|Gs (f )| =
⎪
1
1
⎪
⎩
for fc < f <4π 2 f 22と近似することができる．したがって，時系列が S(f ) ∼ f −β に従うとき，

図 5．DMA および DFA を用いて，白色ノイズを解析した場合の F 2 (s) の解析的結果．破線は傾き 1 となる理想的な直線．

トレンド除去演算に基づくフラクタル時系列解析の数理75

図 6．DMA の |Gs (f )|2 ．
（a）0 次 DMA．
（b）2 次 DMA．
（c）4 次 DMA．

F 2 (s) =∼となる．fc =

√

 1/2
−1/2

 fc
0

S(f ) |Gs (f )|2 df

f −β



π(s2 − 1)12

2

f 2 df +

 1/2
fc

f −β

1
df
4π 2 f 2

6/(πs) であるので，s  1，−1 < β < 3 のとき，各積分は，それぞれ，
2 √6/(πs)
sβ+1
π(s2 − 1)f −β
f 2 df ∼
12
3−β
0
 1/2
1
sβ+1
f −β
df ∼
√
2
2
4π f
β+1
6/(πs)となる．したがって，−1 < β < 3 のとき，
F 2 (s) ∼ sβ+1となる．F 2 (s) ∼ s2α と比較して，(4.24)α=β+12が成り立つ．
また，β > 3 の場合，上の第 1 項の積分は発散する．有限長時系列を考える場合，積分において第 1 項の f = 0 に近い領域の面積が支配的になり，s に依存するカットオフ周波数 fc の影，
響が無視できるため（非積分関数にある s 以外に s に依存する部分がなくなるため）(4.25)

F 2 (s) ∼ s4となる ．したがって．β > 3 の領域では，α = 2 となる．この値が，0 次 DMA における，α の検出上限である．一般に m 次 DMA では，検出上限は m + 2 になる．ランダムウォーク解析の検出上限が 1 であったことから，トレンド除去演算の導入により，検出上限が拡張できることがわかる．詳細については，論文（Kiyono, 2015; Kiyono and Tsujimoto, 2016b）を参照していただきたい．
4.5 スケーリング解析の例
数値的に生成した時系列に対して，DFA および DMA を適用した例を紹介する．解析するサ（a（）e）のように，パワースペクトル密度が S(f ) ∼ f −0.5（図 7（b））
ンプル時系列 {x[i]} は，図 7−1.0に比例し，正規分布に従う過程（ガウス過程），および，パワースペクトル密度が S(f ) ∼ f
（図 7
（f）
）
に比例し，対数正規分布に従う過程（非対称な非ガウス過程）の 2 例である．これらの

76

統計数理

第 71 巻 第 1 号 2023

図 7．1/f β 型ゆらぎのサンプル時系列（（a）（e））とそのパワースペクトル密度
（（b）
（f））．
（上
段）パワースペクトル密度が S(f ) ∼ f −0.5 に比例し，正規分布に従う過程．（下段）パワースペクトル密度が S(f ) ∼ f −1.0 に比例し，対数正規分布に従う過程. 右側は，破（（c）（g））とそのパワースペクト線で示した 2 次関数トレンドを加えたサンプル時系列
（h））．
ル密度
（
（d）

時系列は，まず，独立同分布過程のサンプル時系列を生成し，そのパワースペクトル密度をS(f ) ∼ f −β 型に変換した後に，逆フーリエ変換することで生成した．
（c（）g）のように，これらの時系列に，2 次関数で記述されるトレンド成分（図中
ここでは，図 7（d）（h）のように，の破線）を加えたものを解析した．2 次関数トレンドを含む時系列では，図 7
トレンド成分に起因する S(f ) ∼ f −2.0 型のパワースペクトル密度が低周波数側で観測される．
そのため，トレンド除去などの前処理なしのサンプル時系列のフーリエ変換を用いたパワース
ペクトル推定では，確率的変動成分のスケーリング指数を正しく推定することはできない．
2 次関数トレンドを含むサンプル時系列（N = 104 ）を DFA および DMA を用いて分析した結果が図 8 である．図 8 は，log10 s に対する log10 F (s) のプロットである．両対数プロットの傾
きから，スケーリング指数 α を推定できる．パワースペクトル密度 S(f ) ∼ f −β のスケーリン
グ指数 β は，β = 2α − 1 の関係を使って推定できる．この関係式の導出において，分布のガウ（e）のような，非対称な非ガウス過程であってもス
ス性（正規性）
は仮定していないので，図 7
ケーリング指数を推定できる．
（a）
（c）
（e）
（g）のようにス
トレンド除去の次数が低い 1 次 DFA および 0 次 DMA では，図 8（破線で示した直線からのずれ）が見られた．
ケール s が大きい領域で，2 次関数トレンドの影響
（b）
（d）
（f）
（h）の
しかし，トレンド除去の次数を上げた，2 次 DFA および 2 次 DMA では，図 8
ように，トレンド成分の影響を完全に除去したスケーリング指数の推定ができている．DFA とDMA では，トレンド除去の次数を変化させて，log10 s に対する log10 F (s) のプロットを描く
ことで，時系列に含まれるトレンド成分の影響を評価することができる．
5.

長時間相互相関解析への拡張

前 節 ま で は 単 変 量 時 系 列 {x[i]} の 1/f β 型 ゆ ら ぎ を 考 え た ．こ こ で は ，2 変 量 時 系 列

トレンド除去演算に基づくフラクタル時系列解析の数理77

図 8．図 7（c）（g）に示した過程のサンプル時系列（N = 104 ）を DFA および DMA で解析（d）2 次関数トレンドを含む f −0.5 型ゆらぎ
（図 7（c））の解析結果．
した結果．（a）–（h）2 次関数トレンドを含む f −1.0 型ゆらぎ
（図 7（g））の解析結果．（上段）1 次（e）–
DFA（DFA1）と 2 次 DFA（DFA2）の結果．（下段）0 次 DMA（DMA0）と 2 次 DMA（DMA2）の結果．両対数プロットの直線領域の傾きがスケーリング指数 α の推定値．
スケーリング指数 β は，β = 2α − 1 の関係を使って推定可能．スケーリング指数の理論的値は，左側 2 列では，α = 0.75，β = 0.5 であり，右側 2 列では，α = 1，β = 1．

{(x(1) [i], x(2) [i])} のクロススペクトルに見られる 1/f β 型特性を分析する方法を紹介する（Nakata et al., 2019）．この現象は長時間相互相関として知られるものである（Podobnik andStanley, 2008）．前提として，{x(1) [i]} と {x(2) [i]} は，それぞれが，1/f β 型ゆらぎを示す．長
時間相互相関ではさらに，{x(1) [i]} と {x(2) [i]} のクロススペクトル |S12 (f )| が，|f |  1/2 の領
域で，
1
(5.1)
|S12 (f )| ∼|f |βとなる．{x(1) [i]} と {x(2) [i]} が互いに無相関であれば，式（5.1）の関係は見られない．しかし，両者が同じ成分に駆動されている場合は，式（5.1）の関係が見られる．
式（5.1）の特性を評価する方法として，ここでは，detrending moving-average cross-correlation（DMCA）
を紹介する（Nakata et al., 2019）．出発点は，式（4.7）の代わりに，
analysis2F12(s) = Δs y (1) [i]Δs y (2) [i](5.2)

を考えることである．ここで，(5.3)

Δs y () [i] =sws (j) x() [i + j],( = 1, 2)

j=1である．{x(1) [i]} と {x(2) [i]} の相互共分散関数 C12 (k)，クロスパワースペクトルを S12 (f ) とす
れば，
(5.4)

2
F12
(s) =

s


C12 (k) L(k, s)

k=−s

(5.5)

2
F12
(s) =

 1/2

|S12 (f )| |Gs (f )|2 df

−1/2

統計数理

78

第 71 巻 第 1 号 2023が成り立つ．ここで，L(k, s)，Gs (f ) は，単変量の場合と同じく，それぞれ，式（4.11），（4.12）
で与えられる．したがって，前節と同様に，|S12 (f )| ∼ |f |−β が成り立つ場合は，F12 (s) ∼ sαとなり，α = (β + 1)/2 が成り立つ．
一般には，2 変量時系列の関係については，時間差（ラグ）κ が存在する場合があるので，式（5.2）の代わりに，
2
F12(s, κ) = Δs y (1) [i]Δs y (2) [i + κ]

(5.6)とする必要がある．この場合，式
（5.4）
は，
2
F12
(s, κ) =

(5.7)

s


C12 (k + κ) L(k, s)k=−sとなる．しかし，式（5.5）
は同じままである．基本的な形は，前節で議論した単変量時系列と同じであるので，前節の結果は，2 変量時系列の場合も成り立つ．実際の時系列解析における指を参照していただきたい．
針は論文（Nakata et al., 2019）
6.

おわりに

本稿では，DFA や DMA などのトレンド除去演算を含むスケーリング解析法の特性を，解析的に考察した．スケーリング解析手法の中では，DFA が最も知名度が高く，論文からの引用
は，現時点で 3000 を超えている．それに比べれば，DMA は少数派である．しかし，数理的な
および
（4.7）
の枠組みが素直であるので，DFA と比べて，DMA の方が記述としては，式（4.6）
数理的にきれいな印象がある．加えて，DFA では，区分的曲線フィットによるトレンド除去
が非線形フィルタであるため，元の時系列の分布構造を歪めたり，各点の重みが不均一な統計．これらの欠点
量になっていたりという，明らかな欠点がある（Kiyono and Tsujimoto, 2016a）
は，DMA にはない．とはいえ，DFA でも，DMA でも，パワースペクトル解析でも，正しく使えば，同様の結果が得られる．時系列解析を行う際には，各解析法の特性を理解し，正しく使い，正しく解釈することが重要である．

参

考

文

献

Beran, J. (2017). Statistics for Long-memory Processes, Chapman & Hall/CRC Press, New York.
Carbone, A. and Kiyono, K. (2016). Detrending moving average algorithm: Frequency response andscaling performances, Physical Review E, 93(6), 063309.
Chen, Z., Ivanov, P. C., Hu, K. and Stanley, H. E. (2002). Eﬀect of nonstationarities on detrended
ﬂuctuation analysis, Physical Review E, 65(4), 041107.
Flandrin, P. (1989). On the spectrum of fractional Brownian motions, IEEE Transactions on Information Theory, 35(1), 197–199.
日野幹雄 (1996).『スペクトル解析』, 朝倉書店, 東京.
Höll, M., Kiyono, K. and Kantz, H. (2019). Theoretical foundation of detrending methods for ﬂuctuation analysis such as detrended ﬂuctuation analysis and detrending moving average, PhysicalReview E, 99(3), 033305.
Hu, K., Ivanov, P. C., Chen, Z., Carpena, P. and Stanley, H. E. (2001). Eﬀect of trends on detrended
ﬂuctuation analysis, Physical Review E, 64(1), 011114.
Kiyono, K. (2015). Establishing a direct connection between detrended ﬂuctuation analysis and Fourier
analysis, Physical Review E, 92(4), 042925.
Kiyono, K. (2017). Theory and applications of detrending-operation-based fractal-scaling analysis, 2017
International Conference on Noise and Fluctuations (ICNF ), 1–4.

トレンド除去演算に基づくフラクタル時系列解析の数理

79

Kiyono, K. and Tsujimoto, Y. (2016a). Nonlinear ﬁltering properties of detrended ﬂuctuation analysis,
Physica A, 462, 807–815.
Kiyono, K. and Tsujimoto, Y. (2016b). Time and frequency domain characteristics of detrendingoperation-based scaling analysis: Exact DFA and DMA frequency responses, Physical Review E,
94(1), 012111.
Li, M. (2010). Fractal time series — A tutorial review, Mathematical Problems in Engineering, 2010,
157264.
Ma, Q. D., Bartsch, R. P., Bernaola-Galván, P., Yoneyama, M. and Ivanov, P. C. (2010). Eﬀect of
extreme data loss on long-range correlated and anticorrelated signals quantiﬁed by detrended
ﬂuctuation analysis, Physical Review E, 81(3), 031101.
Mandelbrot, B. B. and Van Ness, J. W. (1968). Fractional Brownian motions, fractional noises andapplications, SIAM Review, 10(4), 422–437.
Moreira, J., da Silva, J. K. L. and Kamphorst, S. O. (1994). On the fractal dimension of self-aﬃne
proﬁles, Journal of Physics A: Mathematical and General, 27(24), 8079–8089.
Nakata, A., Kaneko, M., Shigematsu, T., Nakae, S., Evans, N., Taki, C., Kimura, T. and Kiyono, K.
(2019). Generalized theory for detrending moving-average cross-correlation analysis: A practicalguide, Chaos, Solitons & Fractals: X, 3, 100022.
Peng, C.-K., Buldyrev, S. V., Goldberger, A. L., Havlin, S., Sciortino, F., Simons, M. and Stanley,
H. E. (1992). Long-range correlations in nucleotide sequences, Nature, 356(6365), 168–170.
Peng, C.-K., Havlin, S., Stanley, H. E. and Goldberger, A. L. (1995). Quantiﬁcation of scaling exponents
and crossover phenomena in nonstationary heartbeat time series, Chaos: An Interdisciplinary
Journal of Nonlinear Science, 5(1), 82–87.
Podobnik, B. and Stanley, H. E. (2008). Detrended cross-correlation analysis: A new method for analyzing two nonstationary time series, Physical Review Letters, 100(8), 084102.
Rangarajan, G. and Ding, M. (2000). Integrated approach to the assessment of long range correlationin time series data, Physical Review E, 61(5), 4991.
Shao, Y.-H., Gu, G.-F., Jiang, Z.-Q., Zhou, W.-X. and Sornette, D. (2012). Comparing the performanceof FA, DFA and DMA using diﬀerent synthetic long-range correlated time series, Scientific Reports, 2, 835.
Stanley, H., Buldyrev, S., Goldberger, A., Goldberger, Z., Havlin, S., Mantegna, R. N., Ossadnik, S.,
Peng, C.-K. and Simons, M. (1994). Statistical mechanics in biology: How ubiquitous are longrange correlations?, Physica A: Statistical Mechanics and Its Applications, 205(1-3), 214–253.
Tsujimoto, Y., Miki, Y., Shimatani, S. and Kiyono, K. (2016). Fast algorithm for scaling analysis withhigher-order detrending moving average method, Physical Review E, 93(5), 053304.
Ward, L. M. and Greenwood, P. E. (2010). The mathematical genesis of the phenomenon called “1/fnoise”, Banﬀ International Research Station, https://www.birs.ca/workshops/2010/10frg132/
report10frg132.pdf.

80Proceedings of the Institute of Statistical Mathematics Vol. 71, No. 1, 65–80 (2023)

Mathematical Foundation of Detrending-operation-based
Fractal Scaling Analysis
Ken Kiyono
Graduate School of Engineering Science, Osaka University

To characterize long-range correlation, 1/f β ﬂuctuation, and self-aﬃne fractal properties embedded in nonstationary time series, detrended ﬂuctuation analysis (DFA) has
become widely used in the ﬁelds of physics and biomedical time series analysis. In DFA,a detrending operation with piecewise regression is included in the scaling analysis procedure. Such detrending operations have the advantages of avoiding a false estimate of thescaling exponent induced by nonstationary trend components and expanding the rangeof the detectable scaling exponent. DFA has recently been implemented as a package
in R and Python and is often used instead of conventional power spectral analysis. In
this paper, we provide the mathematical basis for the detrending-operation-based scaling
analysis methods, focusing on DFA and its variant, the detrending moving average (DMA)
algorithm.

Key words: Long-range correlation, fractal, long memory, 1/f ﬂuctuation, time series analysis.