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

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

【長時間相関解析】DFAでホワイトノイズを分析するとどうなるか

Detrended Fluctuation Analysis (DFA)で,無相関な信号,つまり,ホワイトノイズを解析したらどうなるかのお話です.DFAに詳しい人の多くは,そんなのスケールをsとして,ゆらぎ関数 (fluctuation function)は,

 F(s) \propto s^{1/2}

だろというかもしれません.sが大きい,漸近的な振舞を考えれば,これは正しいです.でも,全スケールにわたって正しいわけでもないし,具体的な関数形は分かりません.ということで,今回は,F(s)を解析に求めると,という話です.

 先に結論をまとめておくと,解析する観測時系列\{x_i\}の標本分散を\hat{\sigma}^2とすれば,fluctuation function F(s)は,以下のようになります.

  • 1次DFAでは,

  \begin{eqnarray} F(s) &=& \sqrt{\frac{s^2-4}{15 s}} \hat{\sigma}  \end{eqnarray}

  • 2次DFAでは,

  \begin{eqnarray} F(s) &=& \sqrt{\frac{3 s^2-27}{70 s}} \hat{\sigma}  \end{eqnarray}

  • 3次DFAでは,

  \begin{eqnarray} F(s) &=& \sqrt{\frac{2 s^2- 32}{63 s}} \hat{\sigma}  \end{eqnarray}

  • 4次DFAでは,

  \begin{eqnarray} F(s) &=& \sqrt{\frac{5 s^2-125}{198 s}}\hat{\sigma}  \end{eqnarray}

  • 5次DFAでは,

  \begin{eqnarray} F(s) &=& \sqrt{\frac{3 s^2 - 108}{143 s}} \hat{\sigma}  \end{eqnarray}

 上のどの場合でも,sが大きくなると,

 F(s) \propto s^{1/2}

になることが分かります.

 どうやって導出するかは,MarcとHolgerの論文
link.springer.com
に説明してあります.この論文中の式(12)とか、式(17)のW(s)の部分にタイプミスがありますが,計算間違いやタイプミスは誰にでもあるので,読む側が気をつければいいだけです.

解析的結果の導出

 詳しく説明するのは面倒なので,ポイントだけ言います.要は,

  • 時系列の積分は線形演算 (畳み込みで表現可能)
  • 最小二乗法でえられる回帰曲線の計算も線形演算 (畳み込みで表現可能)
  • 定義通りF^2(s)を計算すれば,自己共分散 (自己相関)との関係がでてくる
  • 自己共分散 (自己相関)を,単位インパルス関数 (ラグ0以外は全部ゼロ)として計算

ということです.つまり,面倒くさがらずに,真面目に計算するだけです.真面目に計算といっても,解いただけでは記憶に残らないので,そこにある構造や解釈のイメージをつかむように意識してください.

 n次DFAの場合について,計算の流れだけ具体的に書いておきます.

 まずは,以下の行列B(n)逆行列B^{-1}(n)を計算します.これは,最小二乗法フィットの畳み込み係数を計算しています.

 \begin{eqnarray} B(n) &=& \left(
\begin{array}{cccc}
 \sum _{j=1}^s 1 & \sum _{j=1}^s j & \cdots & \sum _{j=1}^s j^n \\
 \sum _{j=1}^s j & \sum _{j=1}^s j^2 & \cdots & \sum _{j=1}^s j^{n+1} \\
\vdots & \vdots & \ddots & \vdots  \\
 \sum _{j=1}^s j^{n} & \sum _{j=1}^s j^{n+1} & \cdots & \sum _{j=1}^s j^{2 n} \\
\end{array}
\right) \end{eqnarray}

 逆行列B^{-1}(n)の計算を手ですると大変ですが,数式処理ソフトのMathematicaとか,Maximaとかを使えば,屁くらいの手間です.以下では,B^{-1}(n)の要素を\{b^{-1}_{ij}\}とします.

 ここで,以下の関数を定義します.sは部分区間の長さ (スケール)です.

  \begin{eqnarray} P(t, k) &=& \sum_{i=1}^{n+1} t^{i-1} \sum_{j=1}^{n+1} b^{-1}_{ij}  \sum_{r=k}^{s} r^{j-1} \\
&=& \left( \begin{array}{cccc}
 t^0 & t^1 & \cdots & t^{n} \\
\end{array}
\right)  \left(
\begin{array}{cccc}
b^{-1}_{11} & b^{-1}_{12} & \cdots & b^{-1}_{1n} \\
b^{-1}_{21} & b^{-1}_{22} & \cdots & b^{-1}_{2n} \\
\vdots & \vdots & \ddots & \vdots  \\
b^{-1}_{n1} & b^{-1}_{n2} & \cdots & b^{-1}_{nn} \\
\end{array}
\right) \left( \begin{array}{c}
 \sum_{r=k}^{s} r^{0} \\
 \sum_{r=k}^{s} r^{1} \\
\vdots \\
 \sum_{r=k}^{s} r^{n} \\
\end{array}
\right) \end{eqnarray}

 ということで,標本分散\hat{\sigma}^2のホワイトノイズのサンプル時系列のfluctuation functionは,

 \begin{eqnarray} F^2 (s) = \frac{1}{s} \sum_{p=1}^{s} \left\{ \hat{\sigma}^2 \sum_{k=p}^{s} \left(1 - 2 P(k,p) \right) + \sum_{k=1}^{s} P(k,p)^2 \right\} \end{eqnarray}

です.こんな知識,日本で知りたい人は5人もいないと思います.1人でも役に立てばうれしいので書いておきました.

Fluctuation functionのグラフ

 何でF(s)の全体像を考えるのかといえば,ときどき,勘違いで短時間のスケーリング指数や相関特性を議論しちゃう人がいるからです.

 上のF(s)のグラフは下の図のようになります.

ホワイトノイズをDFAで解析したときのゆらぎ関数 (fluctuation function)

 上の図では,sが小さいスケールで直線的になってません.この非直線的な構造は,時系列の性質ではなく,DFAの方法の都合で生じているのです.ですので,短いスケールの傾きから,相関の性質は議論できません.