はじめに

時系列データにおいて長期記憶性という特徴は見られない場合も結構あります。ですが、見られる場合それに対して適切に対処する必要が出てきます。また、存在しているが気づかずに時系列モデリングをしてしまうと、誤ったモデルの構築につながってしまう危険性をはらんでいます。ということで今回の記事では、時系列データに長期記憶性が存在するのかを確かめる方法とその理論に関して記事にしていきたいと思います。
なんの前置きもなく長期記憶性という言葉を使ってしまいましたが、長期記憶生徒は基本的には自己相関が長期化し、遠い過去の値が現在の値または未来の値へ影響を及ぼしている状況のことを指します。では本題の内容に入っていきましょう。

長期記憶性の確認方法

主に5つ程あります。

R/S分析（Rescaled Range Analysis）

まず、概念的な説明をしましょう。時系列データが長期記憶性を持つということは、過去の変動が将来の変動に長期にわたって影響を与え続けることを意味します。R/S分析では、データの累積偏差の範囲（Range）を標準偏差（Standard deviation）で割って基準化することで、この長期的な依存関係を検出します。時間スケールを変えながらこの値を計算し、その関係性を分析することで、時系列の長期記憶性を特徴づけるHurst指数を推定することができます。

では、具体的な数式を用いて説明していきましょう。時系列データ $X_t$ $(t = 1, 2, ..., N)$ が与えられたとき、まず期間 $\tau$ における平均値を計算します

$\bar{X}(\tau) = \frac{1}{\tau}\sum_{t=1}^{\tau} X_t$

次に、累積偏差系列 $Y(t,\tau)$ を計算します。これは、時系列の値と平均値との差の累積を表します

$Y(t,\tau) = \sum_{i=1}^{t} (X_i - \bar{X}(\tau)) \quad (t = 1, 2, ..., \tau)$

この累積偏差系列から、範囲 R(τ) を計算します。これは累積偏差の最大値と最小値の差として定義されます

$R(\tau) = \max_{1 \leq t \leq \tau} Y(t,\tau) - \min_{1 \leq t \leq \tau} Y(t,\tau)$

また、期間 τ における標準偏差 S(τ) を計算します

$S(\tau) = \sqrt{\frac{1}{\tau}\sum_{t=1}^{\tau} (X_t - \bar{X}(\tau))^2}$

これらを用いて、R/S統計量が以下のように定義されます

$\frac{R(\tau)}{S(\tau)} \propto \tau^H$

ここで H は Hurst指数と呼ばれ、時系列の長期記憶性を特徴づける指標となります。実際の分析では、τ を様々な値に変化させて R/S統計量を計算し、両対数グラフにプロットします。このとき、プロットが直線的な振る舞いを示せば、その傾きから Hurst指数を推定することができます。
Hurst指数は 0 から 1 の値をとり、0.5 が無相関なランダムウォークに対応します。H > 0.5 の場合は正の長期記憶性（持続性）を、H < 0.5 の場合は負の長期記憶性（反持続性）を示します。

Detrended Fluctuation Analysis (DFA)

Detrended Fluctuation Analysis (DFA) は、時系列データから長期相関を検出するための手法で、特にトレンドの影響を受けにくい特徴を持っています。この手法は1990年代にPengらによって開発され、当初はDNAの配列における長距離相関の研究に用いられました。
概念的には、DFAは時系列データを様々な時間スケールで区切り、各区間でトレンドを除去した後の揺らぎ（fluctuation）を分析します。従来のR/S分析と比較して、局所的なトレンドの影響を効果的に除去できる点が大きな特徴です。これにより、非定常性の強いデータに対しても安定した解析が可能となります。

数式を用いて具体的な手順を説明していきましょう。まず、時系列データ $x(i)$ $(i = 1, 2, ..., N)$ が与えられたとき、累積和系列 $y(k)$ を計算します

$y(k) = \sum_{i=1}^k [x(i) - \bar{x}$ ]

ここで $\bar{x}$ は時系列の平均値です。次に、この累積和系列を長さ n の非重複の区間に分割します。各区間において、最小二乗法によって局所トレンド $y_n(k)$ を求めます。通常は線形トレンドを仮定しますが、より高次の多項式を用いることも可能です。

各区間での局所トレンドを除去した後の揺らぎ（分散）を計算します

$F^2(n) = \frac{1}{N} \sum_{k=1}^N [y(k) - y_n(k)$ ^2]]

この計算を様々な区間長 n に対して行います。長期相関が存在する場合、揺らぎの大きさ F(n) は区間長 n とべき乗則の関係を示します

$F(n) \propto n^\alpha$

ここで α はスケーリング指数と呼ばれ、時系列の長期相関を特徴づける指標となります。両対数グラフ上で F(n) を n に対してプロットすると、その傾きから α を求めることができます。　　
スケーリング指数 α は以下のような意味を持ちます
α = 0.5 の場合は、無相関な時系列（白色ノイズ）を示します。
0.5 < α < 1 の場合は、正の長期相関の存在を示します。
0 < α < 0.5 の場合は、負の長期相関（反持続性）を示します。
α = 1 の場合は、1/fノイズに対応します。
α > 1 の場合は、非定常な長期相関の存在を示します。
　　
DFAの特筆すべき利点として、以下の点が挙げられます。局所的なトレンドを除去することで、データの非定常性に対して頑健な解析が可能です。また、異なる時間スケールでの相関特性を区別することができ、クロスオーバー現象（異なるスケール領域で異なるスケーリング指数を持つ現象）の検出も可能です。

自己相関関数（ACF）の分析

自己相関関数（ACF: Autocorrelation Function）の分析は、時系列データにおける時間的な依存関係を定量化する最も基本的かつ重要な手法です。（というよりこれが一番知られている手法ですかね）この手法は、データの現在の値と過去の値との間の線形的な関係性を測定することで、時系列の記憶性や周期性を評価します。（コレログラムとかよく見たことありますよね）概念的には、自己相関は時系列データを時間軸方向にシフトさせたときの、元のデータとの類似度を表します。強い正の自己相関は、高い値の後には高い値が、低い値の後には低い値が続きやすい傾向を示します。一方、負の自己相関は、値が交互に上下する傾向を示します。長期記憶性を持つ時系列では、自己相関が時間とともにゆっくりと減衰していく特徴が見られます。

数式を用いて具体的に説明していきましょう。時系列データ $x_t$ $(t = 1, 2, ..., N)$ に対して、ラグ k の自己相関係数 $\rho(k)$ は以下のように定義されます

$\rho(k) = \frac{\mathbb{E}[(x_t - \mu)(x_{t+k} - \mu)]}{\sigma^2}$

ここで、 $\mu$ は時系列の平均値、 $\sigma^2$ は分散、 $\mathbb{E}$ は期待値演算子です。実際のデータ解析では、以下の標本自己相関係数を用います

$r(k) = \frac{\sum_{t=1}^{N-k}(x_t - \bar{x})(x_{t+k} - \bar{x})}{\sum_{t=1}^{N}(x_t - \bar{x})^2}$

ここで、 $\bar{x}$ は標本平均です。長期記憶性を持つ時系列では、自己相関関数が以下のようなべき乗則に従って減衰することが知られています

$\rho(k) \propto k^{2H-2}$

ここで H は Hurst指数です。このべき乗則的な減衰は、長期記憶性の重要な特徴の一つとなっています。一方、短期記憶性を持つ時系列（例えば、AR過程）では、自己相関関数は指数関数的に減衰します

$\rho(k) \propto e^{-\lambda k}$

ここで λ は正の定数です。また、自己相関関数とスペクトル密度関数 S(f) の間には、以下のようなフーリエ変換の関係が成り立ちます

$S(f) = \sum_{k=-\infty}^{\infty} \rho(k)e^{-2\pi i f k}$

この関係は、時系列の周波数領域での特性を理解する上で重要です。長期記憶性を持つ時系列では、低周波数領域でスペクトル密度がべき乗則に従うことが知られています

$S(f) \propto f^{1-2H}$

結局のところ、自己相関関数（ACF）の分析では、ラグを大きくしていったときの自己相関係数の減衰パターンを観察することで長期記憶性を評価できます。通常の短期記憶過程では自己相関が指数関数的に急速に減衰するのに対し、長期記憶性が存在する場合は自己相関がべき乗則に従いゆっくりと減衰します。具体的には、両対数グラフ上で自己相関関数をプロットしたとき、直線的な減衰が観察され、その減衰が ρ(k) ∝ k^2H-2 の形で表される場合、長期記憶性の存在を示唆します。ここでHは0.5より大きい場合に正の長期記憶性を、0.5より小さい場合に負の長期記憶性を示します。

スペクトル密度分析

スペクトル密度分析は、時系列データを周波数領域で解析する手法です。この手法は、時系列データに含まれる周期的な成分や長期記憶性を周波数の観点から評価することを可能にします。フーリエ解析の原理に基づいており、時間領域での自己相関構造を周波数領域で表現することで、データの特性をより明確に把握することができます。
概念的には、スペクトル密度は時系列データをさまざまな周波数の正弦波の重ね合わせとして表現したときの、各周波数成分の強さを表します。長期記憶性を持つ時系列では、低周波数領域でスペクトル密度が特徴的な振る舞いを示し、べき乗則に従うことが知られています。（先ほど出てきたやつ）これは、長期的な変動が支配的であることを示唆しています。

数式を用いて具体的に説明していきましょう。定常時系列 $x_t$ に対して、スペクトル密度関数 $S(f)$ は自己共分散関数 $\gamma(k)$ のフーリエ変換として定義されます

$S(f) = \sum_{k=-\infty}^{\infty} \gamma(k)e^{-2\pi i f k}$

ここで、自己共分散関数は以下のように定義されます

$\gamma(k) = \mathbb{E}[(x_t - \mu)(x_{t+k} - \mu)]$

実際のデータ解析では、有限長のデータに対して周期グラム（periodogram）を用いてスペクトル密度を推定します。周期グラムは以下のように定義されます

$I(f) = \frac{1}{N}\left|\sum_{t=1}^N x_t e^{-2\pi i f t}\right|^2$

長期記憶性を持つ時系列の場合、低周波数領域でスペクトル密度は以下のようなべき乗則に従います（再掲）

$S(f) \propto f^{1-2H}$

Hurst指数 H
スペクトル密度の推定において、重要な考慮事項としてスペクトル・ウィンドウの選択があります。一般的なスペクトル推定量は以下の形で与えられます

$\hat{S}(f) = \sum_{k=-(N-1)}^{N-1} w(k)\hat{\gamma}(k)e^{-2\pi i f k}$

ここで $w(k)$ はスペクトル・ウィンドウ（例えば、Bartlett、Parzen、Tukey-Hannningなど）であり、推定値の統計的な安定性を向上させる役割を果たします。

また、スペクトル密度と自己相関関数の間には、Wiener-Khinchinの定理として知られる重要な関係があります

$\gamma(k) = \int_{-1/2}^{1/2} S(f)e^{2\pi i f k}df$

この関係は、時間領域での相関構造と周波数領域でのスペクトル特性が等価な情報を含んでいることを示しています。
　　
実際の応用面では、スペクトル密度分析は以下のような特徴を持ちます。まず、複数の周期性が混在するデータにおいて、各周期成分を分離して観察することが可能です。
また、ノイズの特性を評価する上でも重要な役割を果たします。例えば、白色ノイズは全周波数帯域で一定のスペクトル密度を示すのに対し、1/fノイズは周波数に反比例するスペクトル密度を示します。
　　
結局のところ、長期記憶性を確認するには、低周波数領域でのスペクトル密度の振る舞いを観察することで長期記憶性を評価できます。両対数グラフ上でスペクトル密度をプロットしたとき、低周波数域で直線的な減衰が見られ、その傾きが -β（β は0から1の間）となる場合、長期記憶性の存在を示唆します。具体的には、この傾き β とHurst指数 H の間に β = 2H - 1 という関係があり、β が大きいほど長期記憶性が強いことを意味します。なお、白色ノイズの場合はスペクトル密度が周波数によらず一定となるため、この傾向との比較も有用です。

ARFIMA（自己回帰分数和分移動平均）モデルの適用

ARFIMA（Autoregressive Fractionally Integrated Moving Average）モデルは、長期記憶性を持つ時系列データを表現するための統計モデルです。このモデルは、従来のARIMAモデルを拡張し、差分パラメータを整数値から実数値に一般化することで、より柔軟な長期依存性の表現を可能にしています。
概念的には、ARFIMAモデルは時系列データの自己回帰（AR）成分、移動平均（MA）成分、そして分数階差分（Fractional Differencing）成分を組み合わせたものです。分数階差分を導入することで、ARIMAモデルでは表現できない長期記憶性と短期依存性を同時にモデル化することが可能となります。
数式を用いて具体的に説明していきましょう。ARFIMA(p,d,q)モデルは以下のように定義されます

$\phi(B)(1-B)^d x_t = \theta(B)\varepsilon_t$

ここで、B はバックシフト演算子（ $Bx_t = x_{t-1}$ ）を表し、 $\phi(B)$ は自己回帰多項式、 $\theta(B)$ は移動平均多項式、d は分数階差分パラメータ、そして $\varepsilon_t$ は白色ノイズ過程です。自己回帰多項式と移動平均多項式は以下のように表されます

$\phi(B) = 1 - \phi_1B - \phi_2B^2 - ... - \phi_pB^p$

$\theta(B) = 1 + \theta_1B + \theta_2B^2 + ... + \theta_qB^q$

分数階差分演算子は二項展開により以下のように表現されます

$(1-B)^d = \sum_{k=0}^{\infty} \binom{d}{k}(-B)^k = \sum_{k=0}^{\infty} \frac{\Gamma(k-d)}{\Gamma(-d)\Gamma(k+1)}B^k$

ここで、 $\Gamma(\cdot)$ はガンマ関数です。長期記憶性は差分パラメータ d とHurst指数 H の間に以下の関係があります

$H = d + 0.5$

差分パラメータ d の値は時系列の性質を特徴づけます。d が -0.5 から 0 の間にある場合、時系列は反持続的過程となり、負の長期記憶性を示します。d が 0 の場合は短期記憶過程となり、0 から 0.5 の間では長期記憶過程となって正の長期記憶性を示します。0.5 以上では非定常過程となります。
　　
スペクトル密度関数の観点からは、ARFIMAプロセスは低周波数領域で以下のような漸近的な振る舞いを示します

$f(\lambda) \sim c|\lambda|^{-2d} \quad \text{as} \quad \lambda \to 0$

ここで λ は周波数、c は定数です。ARFIMAモデルのパラメータ推定には、最尤推定法が用いられ、その対数尤度関数は以下のように与えられます

$L(\phi,d,\theta|\mathbf{x}) = -\frac{n}{2}\log(2\pi) - \frac{1}{2}\log|\Sigma| - \frac{1}{2}\mathbf{x}'\Sigma^{-1}\mathbf{x}$

ここで Σ は共分散行列です。実際の推定では、この尤度関数を数値的に最大化することで、パラメータの推定値を得ます。予測においては、h期先予測値は条件付き期待値として以下のように表現されます

$\hat{x}{t+h|t} = \mathbb{E}[x{t+h}|x_t,x_{t-1},...]$

最後に

長期記憶性の確認方法からそれから派生して、時系列に関するあれこれの解説も含みましたがどうだったでしょうか？個人的には、長期記憶性の確認方法自体は、単一の結果にとどまることなく、複数の指標を基に評価すべきであると感じています。今回の記事がどなたかの助けとなれば幸いです。最後までお読みいただきありがとうございました！