【時系列モデリング】長時間相関過程をMA過程・AR過程として考えてみる

　このブログでは過去にも説明しましたが、ここで改めて、「長時間相関過程」（long-range dependent process、あるいは、long-range correlated process）を、時系列モデリングで古典的に用いられてきた移動平均過程（moving average process; MA process）と、自己回帰過程（autoregressive process; AR process）を使って考えてみます。

１．MA過程とAR過程の関係

　時系列の特性を数式に表すときの基本モデルとして、MA（移動平均）過程とAR（自己回帰）過程があります。

MA（Moving Average）過程では、現在の値 $\displaystyle{ x _ t }$ は、時々刻々生じるランダムなゆらぎ（ノイズ） $\displaystyle{ \varepsilon _ t }$ の加重平均で表されます。つまり、

$\displaystyle{ x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \cdots + \theta_q \varepsilon_{t-q}. }$

ここで $\displaystyle{ \varepsilon _ t }$ は平均0、分散一定の白色雑音（white noise）です。白色雑音とは、時間的な相関をまったく持たない、純粋なランダムなゆらぎのことです。上の式の $\displaystyle{ q }$ （係数の長さ）は、次数と呼ばれ、 $\displaystyle{ q }$ 次のMA過程は、MA(q)過程と表現します。

AR（AutoRegressive）過程では、現在の値 $\displaystyle{ x _ t }$ は、過去の各時点の観測値 $\displaystyle{ x _ {t-1}, x _ {t-2}, \cdots }$ の線形結合とノイズの和で表されます。つまり、

$\displaystyle{ x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_p x_{t-p} + \varepsilon_t. }$

これは、過去の値と線形の関係で決定論的に依存する部分 $\displaystyle{ \phi _ 1 x _ {t-1} + \cdots + \phi _ p x _ {t-p} }$ と、予測不可能なランダムな要因で駆動される部分 $\displaystyle{ \varepsilon _ t }$ から構成されています。上の式の $\displaystyle{ p }$ （係数の長さ）は、次数と呼ばれ、 $\displaystyle{ p }$ 次のAR過程は、AR(p)過程と表現します。

　式を見れば、MA過程は「過去のノイズの影響」を記述するモデルであり、AR過程は「過去のデータ自体の影響」を記述するモデルのように感じられます。両者は一見すると異なる種類の過程のように見えますが、実際には数学的に密接な関係があり、無限級数を用いることで MA過程をAR表現に変換したり、逆に AR過程をMA表現に変換することも可能です。

　例として、AR(1)過程をMA過程に書き換えてみます。AR(1)過程は、

$\displaystyle{ x_t = \phi x_{t-1} + \varepsilon_t }$

と表され、 $\varepsilon_t$ は平均0、分散 $\displaystyle{ \sigma^ 2 }$ の白色雑音、 $\phi$ は自己回帰係数です。安定な過程である（発散しない）ためには $\displaystyle{ |\phi| \lt 1 }$ である必要があります。この式を後退演算子 $B$ を用いて書くと

$\displaystyle{ (1 - \phi B)x_t = \varepsilon_t }$

となります。

　 $\displaystyle{ |\phi| \lt 1 }$ の条件の下で、

$\displaystyle{ (1 - \phi B)^{-1} = 1 + \phi B + \phi^2 B^2 + \phi^3 B^3 + \cdots }$

が成り立つため（等比級数の和の公式を思い出してください）、上の式の両辺に左から $\displaystyle{ (1 - \phi B)^ {-1} }$ をかけると

$\displaystyle{ x_t = (1 - \phi B)^{-1}\varepsilon_t = (1 + \phi B + \phi^2 B^2 + \cdots)\varepsilon_t }$

となります。これを後退演算子 $B$ を使わずに表現すれば

$\displaystyle{ x_t = \varepsilon_t + \phi \varepsilon_{t-1} + \phi^2 \varepsilon_{t-2} + \phi^3 \varepsilon_{t-3} + \cdots }$

という無限次の移動平均過程（MA(∞)過程）が得られます。したがって、AR(1)過程のMA表現における係数は

$\displaystyle{ a_k = \phi^k }$

であり、これらは指数関数的（exponential）に減衰します。つまり、古いノイズの影響は急速に小さくなり、AR(1)過程は短時間相関しか持ちません。指数的減衰とべき的減衰の違いについては、過去の記事（下記のリンク）も参考にしてください。

長期記憶過程は「ちりも積もれば無限大」ってこと - ケィオスの時系列解析メモランダム

　一般に、AR過程は、無限次のMA過程に書き直せます。差分方程式を使った時系列モデリングの世界では、MA過程として表現できるものが線形過程ですので、AR過程も、今回扱うARFIMA過程も、線形過程です。「長時間相関は非線形」などという記述をたまに見ますが、そうとは言い切れないので注意してください。

２．MA過程とAR過程のパワースペクトルの計算

　MA過程とAR過程は、線形な確率過程ですので、その特性は、パワースペクトル（power spectrum）により完全に特徴付けられます。パワースペクトルは、「各周波数成分がどれくらいの強さ（パワー）を持っているか」を示す関数です。確率過程と、周期成分の関係がイメージできない人もいるかもしれませんが、その場合は過去の記事（下記のリンク）を参考にしてください。

《過去記事まとめ》フーリエ変換とパワースペクトル推定：パワースペクトルは時系列の「DNA」 - ケィオスの時系列解析メモランダム

　ここでは、MA過程とAR過程のパワースペクトルが簡単に（知識がないと難しく感じるかも）、あるいは、公式的に計算できることを説明します。

　まず、MA(1)過程を例にとると、

$\displaystyle{ x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} }$

のフーリエ変換を考えると、そのスペクトル密度関数（パワースペクトル）は

$\displaystyle{ S(\omega) = \sigma^2 |1 + \theta_1 e^{-i\omega}|^2 }$

となります。ここで、 $\displaystyle{ \sigma^ 2 }$ はノイズ成分 $\displaystyle{ \varepsilon _ t }$ の分散です。

　一般に、MA(q)過程のスペクトルは

$\displaystyle{ S(\omega) = \sigma^2 \left| 1 + \theta_1 e^{-i\omega} + \cdots + \theta_q e^{-iq\omega} \right|^2 }$

と表されます。つまり、MA過程の係数 $\displaystyle{ \left\{ \theta _ 1, \theta _ 2, \cdots, \theta _ q \right\} }$ とノイズ成分 $\displaystyle{ \varepsilon _ t }$ の分散 $\displaystyle{ \sigma^ 2 }$ が決まれば、パワースペクトルが計算できます。時系列からこれらのパラメタを推定することが、時系列モデリングを使った時系列解析の目的です。

　一方、AR(p) 過程のスペクトルは、分母に多項式が現れます：

$\displaystyle{ S(\omega) = \frac{\sigma^2}{\left|1 - \phi_1 e^{-i\omega} - \cdots - \phi_p e^{-ip\omega}\right|^2}. }$

　以上からわかるように、MAモデルは分子側に多項式があり、ARモデルは分母側に多項式があります。両者には、分子と分母の違いがあるので、何か全然違うもののような印象を与えます。

後退演算子Bを用いた表現

　ここでは、パワースペクトルを公式的に計算しやすくするために、後退演算子（backward shift operator）を導入します。後退演算子 $\displaystyle{ B }$ は、時系列を1ステップ前にずらす作用を持ち、次のように定義されます：

$\displaystyle{ B x_t = x_{t-1}. }$

この演算子を用いると、MA(q)過程とAR(p)過程はそれぞれ次のように書けます。

MA(q)過程（Moving Average process）

$\displaystyle{ x_t = \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \cdots + \theta_q \varepsilon_{t-q}. }$

は、後退演算子を用いて、

$\displaystyle{ x_t = (1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q)\varepsilon_t }$

と表現することができます。さらに、

$\displaystyle{ \Theta(B) = 1 + \theta_1 B + \cdots + \theta_q B^q }$

のように置けば、MA(q)過程は、

$\displaystyle{ x_t = \Theta(B)\varepsilon_t, }$

と表すことができます。

AR(p)過程（AutoRegressive process）

$\displaystyle{ x_t = \phi_1 x_{t-1} + \phi_2 x_{t-2} + \cdots + \phi_p x_{t-p} + \varepsilon_t. }$

は後退演算子を用いて、

$\displaystyle{ (1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p)x_t = \varepsilon_t, }$

と表現することができます。さらに、

$\displaystyle{ \Phi(B) = 1 - \phi_1 B - \cdots - \phi_p B^p }$

のように置けば、AR(p)過程は、

$\displaystyle{ \Phi(B)x_t = \varepsilon_t, }$

と表すことができます。

　このように、MA過程とAR過程はともに「多項式演算子」を使って統一的に表現できます。MA過程ではノイズ側に多項式演算子 $\displaystyle{ \Theta(B) }$ が作用し、AR過程ではデータ側に多項式演算子 $\displaystyle{ \Phi(B) }$ が作用しています。

　この記法を使うと、パワースペクトルも非常に簡潔に、公式的に書けます。つまり、

$\displaystyle{ S(\omega) = \sigma^2 \frac{|\Theta(e^{-i\omega})|^2}{|\Phi(e^{-i\omega})|^2}. }$

となります。

　この式は、ARMA過程（Autoregressive Moving Average process）と呼ばれる過程の一般形です。MAモデルでは $\displaystyle{ \Phi(B)=1 }$ 、ARモデルでは $\displaystyle{ \Theta(B)=1 }$ となり、どちらもこの式に含まれることがわかります。

３．長時間相関を示す ARFIMA(0,d,0) 過程のパワースペクトルの漸近形

　ARFIMA（Autoregressive Fractionally Integrated Moving Average）過程は、ARやMAを一般化したモデルで、「分数階差分」という考えを取り入れています。ARFIMA(0,d,0) 過程は次のように定義されます：

$\displaystyle{ (1 - B)^d x_t = \varepsilon_t, }$

ここで、 $\displaystyle{ B }$ は先ほど導入した後退演算子であり、 $\displaystyle{ d }$ は「分数階の差分」を表すパラメータで、0 < d < 0.5 の範囲では定常で、長時間相関（long-range dependence）を示します。

　このとき、上で紹介した公式を使えば、パワースペクトルは、

$\displaystyle{ S(\omega) = \frac{\sigma^2}{|1 - e^{-i\omega}|^{2d}}. }$

となります。

0に近い小さな周波数 $\displaystyle{ \omega \to 0 }$ のとき、

$\displaystyle{ |1 - e^{-i\omega}| \approx |\omega| }$

なので、

$\displaystyle{ S(\omega) \propto \omega^{-2d}. }$

となることがわかります。

　つまり、この過程は、低周波領域で $\displaystyle{ 1/f^ {-2d} }$ 型のパワースペクトルを示す長時間相関過程になっています。

４．ARFIMA(0, d, 0) 過程の MA表現と AR表現

　ARFIMA(0,d,0) 過程は、MA過程、あるいは、AR過程としての表現することができます。これらの２つの表現のどちらかを使えば、ARFIMA(0,d,0) 過程のサンプル時系列を生成することができます。

MA表現（移動平均形）　ARFIMA(0, d, 0)過程は、次のような無限次の移動平均過程で表すことができます。

$\displaystyle{ x_t = \sum_{k=0}^{\infty} a_k \varepsilon_{t-k}, \quad a_k = \frac{\Gamma(k + d)}{\Gamma(k + 1)\Gamma(d)}. }$

ここで、 $\displaystyle{ \Gamma(\cdot) }$ はガンマ関数（階乗を一般化した関数）です

　係数 $\displaystyle{a _ k}$ の大きさは $\displaystyle{ k }$ が大きくなるにつれて次のように漸近的に減少します。

$\displaystyle{ a_k \sim \frac{1}{\Gamma(d)}\,k^{\,d-1}, \quad (k \to \infty). }$

AR表現（自己回帰形）　ARFIMA(0, d, 0)過程は、次のような無限次の自己回帰過程としても表すことができます。

$\displaystyle{ x_t = \sum_{k=1}^{\infty} \phi_k\, x_{t-k} + \varepsilon_t, \quad \phi_k = (-1)^{k+1}\frac{\Gamma(k - d)}{\Gamma(k + 1)\,\Gamma(-d)}. }$

　係数 $\displaystyle{ \phi _ k }$ は $\displaystyle{ k }$ が大きいとき、次のように漸近近似できます。

$\displaystyle{ \phi_k \sim \frac{(-1)^{k+1}}{\Gamma(-d)}\,k^{-(d+1)},\qquad (k\to\infty). }$

５．まとめ

　今回の記事では、時系列モデルの基本であるMA（移動平均）過程とAR（自己回帰）過程を出発点として、それらを統一的に理解するための数式表現やパワースペクトルの関係を確認しました。両者は一見まったく異なるモデルのように見えますが、後退演算子（backward shift operator; B）を導入すると、どちらも多項式演算子で表される同一の枠組みに含まれることがわかります。この演算子表現により、ARとMAを組み合わせたARMA過程や、その一般化であるARFIMA過程を自然に導入できます。

　特に、ARFIMA(0, d, 0)過程は、分数階差分 $\displaystyle{ (1 - B)^ d }$ を導入することで、過去の影響が指数関数的ではなく、べき乗的にゆっくりと減衰することを表現できるモデルです。MA表現の係数 $\displaystyle{a _ k}$ は漸近的に $\displaystyle{a _ k \sim k ^ {d-1}}$ に従い、過去の影響の蓄積が無限になる効果を持ちます。この「無限に蓄積する記憶」が、長時間相関（long-range dependence; LRD）の本質です。

　パワースペクトルの観点から見ると、ARFIMA(0, d, 0) 過程は低周波数領域で $\displaystyle{S(\omega) \propto \omega^{-2d}}$ という 1/f 型スペクトルを示し、時間スケールを超えてゆらぎが自己相似的に現れることを意味します。これは、自然界や生体信号、経済データなど、多くの現象に見られる「1/fゆらぎ（1/f fluctuation）」と深く関係しています。

(おまけ) Rで係数をプロット

　長時間相関過程では、過去の値がどのように現在に影響を与えるかを見るために、ARFIMA(0,d,0)のMA表現とAR表現の係数をRでプロットしてみます。

############################################################
## ARFIMA(0,d,0) の MA係数 a_k と AR係数 φ_k をプロット
## それぞれに漸近形を点線で重ねる
############################################################
## ---- パラメータ設定 ----
d <- 0.30          # 0 < d < 0.5（定常・長時間相関）
K <- 1000          # 最大ラグ（大きめでも安定）

## ---- 係数の計算 ----
## a_k = Gamma(k+d) / [Gamma(k+1) * Gamma(d)]
ak_lgamma <- function(d, K) {
  k <- 0:K
  a <- numeric(K + 1)
  a[1] <- 1.0
  if (K >= 1) {
    kk <- k[-1]
    a[-1] <- exp(lgamma(kk + d) - lgamma(kk + 1) - lgamma(d))
  }
  a
}

## phi_k = (-1)^{k+1} * Gamma(k-d) / [Gamma(k+1) * Gamma(-d)]
phi_lgamma <- function(d, K) {
  k <- 1:K
  const <- -lgamma(-d)
  (-1)^(k + 1) * exp(lgamma(k - d) - lgamma(k + 1) + const)
}

a   <- ak_lgamma(d, K)
phi <- phi_lgamma(d, K)

## ---- 漸近形（理論定数を用いる）----
k_pos     <- 1:K
a_asym    <- (1 / gamma(d))   * (k_pos)^(d - 1)       # a_k ~ k^{d-1}/Gamma(d)
phi_asymA <- (1 / abs(gamma(-d))) * (k_pos)^(-(d + 1)) # |phi_k| ~ k^{-(d+1)}/|Gamma(-d)|

## ---- 可視化 ----
op <- par(no.readonly = TRUE)
par(mfrow = c(1, 2), mar = c(5, 6.0, 3.0, 1))

## 左パネル：a_k の log-x プロット
plot(k_pos, a[k_pos + 1], pch = 16, log = "x", col = 2,
     xlab = expression(k),
     ylab = expression(a[k] == frac(Gamma(k + d), Gamma(k + 1) * Gamma(d))),
     main = bquote(MA~coefficients~~a[k]~~"(ARFIMA(0," * .(d) * ",0))"))
lines(k_pos, a_asym, lty = 2, lwd = 2, col = 4)  # 漸近形（破線）
legend("topright",
       legend = c(expression(a[k]),
                  expression(a[k] %~~% k^{d-1} / Gamma(d))),
       pch = c(16, NA), lty = c(NA, 2), lwd = c(1, 2),
       col = c(2, 4), bty = "n", cex = c(1, 0.9))

## 右パネル：|phi_k| の log-x プロット（符号を点種で表現）
y_abs   <- abs(phi)
pch_sig <- ifelse(phi > 0, 16, 1)
plot(k_pos, y_abs, log = "x", type = "p", pch = pch_sig, col = 2,
     xlab = expression(k),
     ylab = expression("|" * phi[k] * "|" == frac(Gamma(k - d),
                       Gamma(-d) * Gamma(k + 1)) ~ "(abs)"),
     main = bquote(AR~coefficients~~phi[k]~~"(sign-coded)"))
lines(k_pos, phi_asymA, lty = 2, lwd = 2, col = 4)  # |漸近形|（破線）
legend("topright",
       legend = c(expression("|" * phi[k] * "| (positive)"),
                  expression("|" * phi[k] * "| (negative)"),
                  expression("|" * phi[k] * "|" %~~% k^{-(d+1)} / "|" * Gamma(-d) * "|")),
       pch = c(16, 1, NA), lty = c(NA, NA, 2),
       col = c(2, 2, 4), bty = "n", cex = 1, lwd = c(1, 1, 2))

par(op)