時系列解析をゆるっと解説(1)

動機

時系列データ分析を絶賛勉強中
時系列解析の全体像がまとまっている資料を作りたかった

導入

参考書「時系列解析入門(北川源四郎著)」
構造化を目的としているので個々の手法等の詳細については説明していません。あくまでも概要です。

時系列解析

参考書の内容を構造化してまとめた

動機
導入
時系列解析

時系列解析の目的

記述(description)
時系列を図示したり、標本自己共分散関数、標本自己相関関数、ピリオドグラムなどの木ジョン的な記述統計量を用いて時系列の特徴を完結に表現すること、グラフによって表現すること。
モデリング(modeling)
与えられた時系列に対し、その変動の仕方を表現する時系列モデルを構成し、時系列の確率的構造を解析すること。解析の対象や目的に応じて適切な時系列モデルを選択し、そのモデルに当てはまるパラメータを推定することが必要。
予測(prediction)
時系列が互いに相関を持つことを利用して、現在までに得られた情報から今後の変動を予測すること。
信号抽出(signal extraction)
時系列から目的に応じて必要な信号や情報を取り出すこと。対象の特徴や目的に応じて適切なモデリングを行うことが必要。

前処理

変数変換
あたいの増加に伴って変動の大きさも増大するという特徴を持っていたり、分布が正規分布とは著しく異なっていたりした場合に $z_n = h(y_n)$ のような変換をすることで扱いやすいデータにする。 ex) ロジット変換, Box-Cox変換
差分
顕著なトレンドがわかる場合は差分系列を求めて解析を行う。複数回繰り返することでn次成分を除去することができる。
前期比, 前年同期比
経済時系列などでは、比を求めて成長率を検出することがある。
移動平均
変動の激しい時系列を滑らかにする簡便な方法。

時系列データの記述方法

共分散関数
-自己共分散関数(標本共分散関数)
時系列に定常性が仮定できると平均値関数は時刻に依存しない一定の値となる。自己共分散関数は $C_k = E[(y_n - \mu)(y_{n-k} - \mu)$ ]としてかけ、kをラグと呼ぶ。 $k=0$ 以外の値が0になる時を特別ホワイトノイズ(white noise)と呼ぶ。
-自己相関関数(標本自己相関関数)
上の $y_n$ とのん相関係数をラグkの関数とみなしたものを自己相関関数と呼ぶ。
-相互共分散関数, 相互相関関数
上記の例は1変量系列の場合だったが、多変量解析の場合は"相互"の関係を記述するようになる。この時関数は $l*l$ 行列(l: 変数の数)で表記される。
スペクトル
-スペクトル
ラグkが大きくなるとき自己共分散関数 $C_k$ が急速に減衰し
$$ \sum_{k=-\infty}^{\infty} C_k \leq \infty $$
となる時フーリエ変換が以下の式で定義できる。パワースペクトルまたはスペクトルと呼ばれる。

$$ p(f) = C_0 + 2 \sum_{k=^\infty}^\infty C_k cos 2 \pi k f $$

ピリオドグラム
スペクトルの式の $C_k$ の代わりに標本自己共分散関数を代入して得られるものをピリオドグラムと呼ぶ。また、定義域を連続区域に拡張したものを標本スペクトルと呼ぶ。
$$ p_j = \hat{C_0} + 2 \sum_{k=1}^{N-1} C_k cos 2 \pi k f $$
ピリオドグラムの平均と平滑化
標本スペクトルではデータ数をいくら大きくしても真のスペクトルに近づかない問題点があった。そこで $f_j = \frac{i} {2 L}$ とすることで定数Lで定められるデータ量の平均を取るようにしてデータ数の増加に伴って真のスペクトルに近づくようにした。
しかし、これだけでは厳密には正しくなく、さらにスペクトルウィンドウと呼ばれる関数を用いて平滑化する。ウィンドウにはハミングとハニングがある。
$$ \hat{p_j} = \sum_{i=-m}^{m} W_i p_j-i $$

モデリングの評価基準とパラメータ推定・モデル選択

KL情報量
カルバック・ライブラー情報量は以下の式で定義される。
$$ I(g;f) = E_y \log {\frac{g(Y)} {f(Y)}} $$ KL情報量が小さいほど確率分布 $f$ は真の確率分布 $g$ に近いと考えることができる。 $B(g;f)=-I(g;f)$ は一般化されたエントロピーと呼ばれ、KL情報量に基づくモデル構築はエントロピー最大化原理である。しかし、KL情報量が実際の統計的モデルの評価に直接用いられることはない。それは統計解析が行われる場面では、真の分布が未知でKL情報量を計算することができないからである。真の分布 $g(y)$ の代わりに $g(y)$ から独立に観測されたデータが観測できているのが現実的な状況である。この場合は、(途中式省略)対数尤度 $\sum_{n=1}^N \log f(y_n)$ を最大にするようなモデルを選べば、近似的にエントロピーを最大にすることができる。対数尤度を最大化する方法として最小二乗法(ハウスホルダー法)やニュートン法などの数値的最適化のアルゴリズムがある。
AIC(赤池情報量規準)
異なるモデル間でも尤度比較できるようにしたのがこの情報量規準。以下の式で与えられる。
$$ AIC = -2 l(\hat{\theta}) + 2k = -2 (最大対数尤度) + 2 (パラメータ数) $$
ちなみに、Box-Cox変換などの変数変換を施した際のAICは以下のように更新される。 $$ z_n = h(y_n) $$ $$ g(h) = | \frac{dh}{dy} | f(h(y)) $$ $$ AIC_{z}' = AIC_z - 2 \sum \log | \frac{dh}{dy} | $$

ARMAモデル(時系列モデリング手法)

ARMAモデル
自己回帰移動平均モデル(AutoRegressive Moving Average model, ARMA)は自己回帰の次数( $m$ )、自己回帰係数( $a_i$ )、移動平均次数( $l$ )と移動平均係数( $b_i$ )を用いるモデル。
$$ y_n = \sum_{i=i}^{m} a_i y_n-i + v_n $$

$$ -\sum_{i=1}^lb_i v_n-i $$

ARMAモデルのインパルス応答関数
$B y_n \equiv y_{n-1}$ によって定義される時間シフトオペレータを用いることで、(途中式省略)ARMAモデルは現在までに加わっ白色雑音 $v_n$ だけを用いて
$$ y_n = g(B) v_n $$
と無限次の移動平均モデルで表現できる。時刻 $n=0$ に加わったノイズが $i$ 期後に時系列に及ぼす影響の大きさを表す $g_i$ を用いてARMAモデルのインパルス応答関数が求められる。

$$ y_n = \sum_{i=0}^{\infty} g_i v_n-i $$

* ARMAモデルの自己共分散関数
インパルス応答関数が与えられると自己共分散関数を求めることができる。特に $l=0$ のARモデルの場合の式はユールウォーカー方程式と呼ばれる。

AR係数とPARCORの関係
m-1次のARモデルの係数 $a_i^{m-1}$ とm次のARモデルの係数 $a_i^{m}$ の間には次のような関係が成り立つ。
$$ a_i^{m} = a_i^{m-1} - a_m^{m} a_{m-i}^{m-1} $$ $a_m^{m}$ はm次のPARCOR(偏自己相関係数)と呼ばれる。m次までのPARCORが与えられると以上の関係を繰り返し適用することんいより2次~M次のARモデルの係数を全て求めることができる。逆の場合も求めることができる。以上からm次のARモデルの係数 $a_1^{m},\cdots, a_m^{m}$ を推定することと、m次までのPARCOR $a_1^{1},\cdots, a_m^{m}$ を推定することは同等である。
パワースペクトル
(詳細略)ARMAモデルのパワースペクトルを求めることができる。AR次数やMA次数とスペクトルの山谷の間には密接な関係があり、k個の山を表現するためには2k以上のAR次数が、k個の谷を表現するためには2k以上のMA次数が必要になる。また、スペクトルの山谷の場所や高さは特性方程式の複素数根の偏角と絶対値に関連する。
ARモデルの推定
与えられたデータにARモデルを当てはめるためには、次数mを決定し、自己回帰係数 $a_m,\cdots,a_m$ と分散 $\sigma^{2}$ を推定することが必要となる。そのためのアプローチとしては、
1.数値的方法で最尤推定パラメータ $\hat{\theta}$ (尤度最大化)→AIC最小化法で次数mの評価
2.ユールウォーカー法(+レビンソンのアルゴリズム)
3.最小二乗法
4.PARCOR法

局所定常ARモデル(時系列モデリング手法)

局所定常ARモデル
現実の時系列は非定常なものが多い。非定常時系列のモデリングとして局所定常ARモデルがあり、それは時間空間を適当な小区間に分割し、各小区間では定常と仮定してARモデルを当てはめることで非定常時系列を近似的に表現するモデル。
任意個の区間への自動分割
(1)定常性を仮定する小区間の長さLと各区間で当てはめるARモデルの最高次数mを決める。
(2) $i$ 番目の区間にm次までのARモデルを当てはめ、 $AIC_0$ を求める。
(3) $i+1$ 番目の区間にm次までのARモデルを当てはめ、 $AIC_1$ を求め、 $AIC_D = AIC_0 + AIC_1$ とする。
(4) $i, i+1$ の区間を統合してm次までのARモデルを当てはめて $AIC_P$ を求める。
(5)AICの値を比べて(3),(4)どちらのモデルが良いか判断し、次の区間に進み、同様のステップを踏む。
以上のステップにより分割ができる。データの追加のたびに計算する必要があるが、ハウスホルダー変換の性質を使うと効率よく計算することができる。
変化時点の精密な推定
前項では区間ごとにARモデルを当てはめて分割を行なったが、区間内で変化していることまではわからない。よってより精密に変化時点を知るためには $[1,n - 1$ , [n,N]]の二つにさらに分割しそれぞれにARモデルを当てはめて変化時点の決定をする。 $n$ を変化させてそれぞれのAICを比較するので大量の計算をようすることになるがハウスホルダー変換のデータ追加の方法を利用すると効率よく計算できる。