[R] [Stan] でベクトル ARMA (VARMA) を推定

概要

書きたいネタは他にもいくつかあるが, 下準備に時間の掛かりそうなものばかりなので, とりあえず簡単なネタで投稿. Stan は最近使い始めたばかりなので非効率なところがあるかも*1.

時系列データのモデリングでよく使われる, 自己回帰移動平均 (ARMA) 過程のベクトル版であるベクトルARMA (VARMA) モデルの簡単な説明をする
VARMA モデルを stan で推定する場合のサンプルコードも紹介アンド解説 (RStan を利用)

VARMA モデルとは

ARMA モデルについては, 専門書では Hayashi (2000), 沖本 (2010) などで取り上げられている*2し, より簡易な解説もネット上に多く存在する. 数式を用いてある程度厳密に説明した記事だと, tjo.hatenablog.com
tjo.hatenablog.com

とかウィキペ自己回帰移動平均モデル - Wikipedia とかでも概要がわかる. 手前味噌だが ill-identified.hatenablog.com でもよい. ARMA が ARIMA に拡張できるように VARIMA も存在するが, 今回は VARMA の話をする.

VAR から VARMA へ

ここからが読み飛ばすとわけが分からなくなる説明*3. まずは MA の要素のない VAR の場合から. $y_{1},\,\cdots,y_{n}$ という $n$ 個の時系列があって, それぞれが自己回帰モデルに従うとする.

$\begin{align} y_{1,t}= & \phi_{1,1}y_{1,t-1}+\phi_{1,2}y_{t-2}+\cdots+\varepsilon_{1,t}\\ y_{2,t}= & \phi_{2,1}y_{2,t-1}+\phi_{2,2}y_{t-2}+\cdots+\varepsilon_{1,t}\\ \vdots\end{align}$

複数の時系列をこのような連立方程式で表せるのだが, 上記の式ではそれぞれの時系列は他の時系列と全く別個に動いていることになる. 現実には, 複数の時系列は互いに何らかの相関関係があることが多い. それは以前の記事
ill-identified.hatenablog.com
で紹介した CAR モデルのような空間的なつながりによるものかもしれないし. あるいはもっと簡単に, 気温の推移とアイスクリームやおでんの売り上げの推移のような関係かもしれない. とにかく, それぞれの時系列どうしで相関が全くないことを頭から否定できる場合は少ない*4. そこで VAR では, 各時系列は, 自己の過去の値だけでなく, 他の全ての時系列の過去の値によって決まるように,

$\begin{align} \begin{bmatrix}y_{1,t}\\ y_{2,t}\\ \vdots\\ y_{n,t} \end{bmatrix}= & \begin{bmatrix}\psi_{1,1} & \cdots & \psi_{1,n}\\ \vdots & \ddots & \vdots\\ \psi_{n,1} & \cdots & \psi_{n,n} \end{bmatrix}\begin{bmatrix}y_{1,t-1}\\ y_{2,t-1}\\ \vdots\\ y_{n,t-1} \end{bmatrix}+\cdots+\begin{bmatrix}\varepsilon_{1,t}\\ \vdots\\ \varepsilon_{n,t} \end{bmatrix}\tag{1}\label{eq:VAR}\end{align}$

のようにベクトルで表現する. 1期前の部分だけを書き直すと,

$\begin{align} y_{1,t}= & \psi_{1,1}y_{1,t-1}+\psi_{1,2}y_{2,t-1}+\cdots+\psi_{1,n}y_{n,t-1}+\varepsilon_{1,t}\\ y_{2,t}= & \psi_{2,1}y_{2,t-1}+\psi_{2,2}y_{2,t-2}+\cdots+\psi_{2,n}y_{n,t-1}+\varepsilon_{1,t}\\ & \vdots\end{align}$

となる (つまり VAR(1) ). さらに 2期, 3期前と書き足せるが, 見づらくなるので, ベクトルを使って表現する.

$\begin{align} \boldsymbol{y}_{t}= & \begin{bmatrix}y_{1,t}\\ \vdots\\ y_{n,t} \end{bmatrix},\\ \boldsymbol{\varepsilon}_{t}= & \begin{bmatrix}\varepsilon_{1,t}\\ \vdots\\ \varepsilon_{n,t} \end{bmatrix}\end{align}$

というふうに, 同じ $t$ 時点の時系列どうしを1つの列ベクトルにまとめると,

$\begin{align} \boldsymbol{y}_{t}= & \boldsymbol{\Psi}_{1}\boldsymbol{y}_{t-1}+\boldsymbol{\Psi}_{2}\boldsymbol{y}_{t-2}+\cdots+\boldsymbol{\varepsilon}_{t}\tag{\ref{eq:VAR}'}\label{eq:VAR2}\end{align}$

となり, あたかもベクトルでない単なる AR のように書ける. この(\ref{eq:VAR2}) が VAR モデルである.

次に, 移動平均 (MA) 過程は

$\begin{align} y_{t}= & \theta_{0}+\theta_{1}\varepsilon_{t-1}+\theta_{2}\varepsilon_{t-2}+\cdots\end{align}$

のように過去の誤差項の値で決定されていたから, これも VAR と同様に, $n$ 個の時系列 $y_{1,t},\,\cdots,\,y_{n,t}$ は自身の過去の誤差項だけでなく, 他の時系列の過去の誤差項を加味して決まると考えると, ベクトル移動平均 (VMA) 過程は

$\begin{align} \begin{bmatrix}y_{1,t}\\ y_{2,t}\\ \vdots\\ y_{n,t} \end{bmatrix}=\begin{bmatrix}\theta_{0,1}\\ \vdots\\ \theta_{0,n} \end{bmatrix}+ & \begin{bmatrix}\theta_{1,1} & \cdots & \theta_{1,n}\\ \vdots & \ddots & \vdots\\ \theta_{n,1} & \cdots & \theta_{n,n} \end{bmatrix}\begin{bmatrix}\varepsilon_{1,t-1}\\ \varepsilon_{2,t-1}\\ \vdots\\ \varepsilon_{n,t-1} \end{bmatrix}+\cdots\tag{2}\label{eq:VMA}\end{align}$

と書ける. よって, こちらもベクトル $\boldsymbol{y}_{t}$ と $\boldsymbol{\varepsilon}_{t}$ を使って,

$\begin{align} \boldsymbol{y}_{t}=&\boldsymbol{\theta}_{0}+\boldsymbol{\Theta}_{1}\boldsymbol{\varepsilon}_{t-1}+\boldsymbol{\Theta}_{2}\boldsymbol{\varepsilon}_{t-2}+\cdots+\boldsymbol{\Theta}_{q}\boldsymbol{\varepsilon}_{t-q}+\boldsymbol{\varepsilon}_{t}\tag{2'}\label{eq:VMA2} \end{align}$

と表現できる. 以上 (\ref{eq:VAR2}), (\ref{eq:VMA2}) から, VARMA は一般に

$\begin{align} \boldsymbol{y}_{t}=& \boldsymbol{\Psi}_{1}\boldsymbol{y}_{t-1}+\cdots+\boldsymbol{\Psi}_{p}\boldsymbol{y}_{t-p}\\ &+\boldsymbol{\theta}_{0}+\boldsymbol{\Theta}_{1}\boldsymbol{\varepsilon}_{t-1}+\cdots+\boldsymbol{\Theta}_{q}\boldsymbol{\varepsilon}_{t-q}+\boldsymbol{\varepsilon}_{t}\tag{3}\label{eq:VARMA}\end{align}$

と書くことができる. なお, 複数時系列があるので, $p,\,q$ はそれぞれ, 全ての時系列のうちで一番長くラグをとる時系列に対応する. 例えば $1$ 番目が ARMA $(2,\,2)$ , 2番目が ARMA $(3,\,1)$ , 残りが全て AR (1) ならば, VARMA ( $3,\,2$ ) になる. ラグの短いものは対応する係数がゼロになる. スカラをベクトルに変えただけなので, 差分も同じように定義できるから, VARIMA について考えることもできるが, 今回は VARMA だけを扱う.

VARMA の反転条件

VARMA の式が (\ref{eq:VARMA}) で与えられたので, 次に ARMA と同じように定常的 (安定的) な VARMA かどうかを判定する条件が存在するかどうかを考える必要がある. Hayashi (2000) の畳み込みを使った説明を借りると簡単になる (かもしれない) . まず, 一旦ベクトルでない ARMA の話に戻る. まず, 前提として MA $(\infty)$ は定常的であるとする*5. ARMA は, $y_{t}$ のラグ項を左に移動して,

$\begin{align} y_{t}-\psi_{1}y_{t-1}-\cdots-\psi_{p}y_{t-1}=\theta_{0}+ & \theta_{1}\varepsilon_{t-1}+\cdots+\theta_{q}\varepsilon_{t-q}+\varepsilon_{t}\tag{4}\label{eq:arma2}\end{align}$

となる. ここでラグ多項式を導入する. ラグ多項式とは,

$\begin{align} \psi(\mathrm{L}):= & 1-\psi_{1}\mathrm{L}-\psi_{2}\mathrm{L}^{2}-\cdots-\psi_{p}\mathrm{L}^{p}\tag{5}\label{eq:lagpoly}\end{align}$

のように係数とラグ演算子のべき乗による多項式である (ラグ演算子を入れた場合は厳密にはフィルタと呼ぶ). ラグ演算子 $\mathrm{L}$ は, $\mathrm{L}y_{t}=y_{t-1}$ , $\mathrm{L}^{2}y_{t}=y_{t-2}$ , ..., $\mathrm{L}^{p}y_{t}=y_{t-p}$ というふうに扱われるので, ラグ多項式と $y_{t}$ をかけた $\psi(\mathrm{L})y_{t}$ は (\ref{eq:arma2}) の左辺と等しい. 同様に, 右辺の MA についても,

$\begin{align} \theta(\mathrm{L})= & 1+\theta_{1}\varepsilon_{t-1}+\cdots+\theta_{q}\varepsilon_{t-q}\tag{6}\label{eq:lagpolyma}\end{align}$

と定義すると, (\ref{eq:arma2}) は

$\begin{align} \psi(\mathrm{L})y_{t}= & \theta_{0}+\theta(\mathrm{L})\varepsilon_{t}\end{align}$

となる. ここで, $\psi(1)=:1-\psi_{1}-\psi_{2}-\cdots\psi_{p}$ とするなら, $\mu:=\theta_{0}/\psi(1)$ とおくと, $\theta_{0}$ を移動して,

$\begin{align} \psi(\mathrm{L})y_{t}-\psi(1)\mu= & \theta(\mathrm{L})\varepsilon_{t}\end{align}$

となる. $\mu$ は定数なので, ラグ演算子をかけても変わらないから, $\psi(1)\mu=\psi(\mathrm{L})\mu$ である. よって,

$\begin{align} \psi(\mathrm{L})(y_{t}-\mu)= & \theta(\mathrm{L})\varepsilon_{t}\end{align}$

となる. $\psi(\mathrm{L})$ の逆数が存在するなら,

$\begin{align} y_{t}-\mu= & \psi(\mathrm{L})^{-1}\theta(\mathrm{L})\varepsilon_{t}\end{align}$

となる. $\psi(\mathrm{L})^{-1}\theta(\mathrm{L})$ は新たなラグ多項式となるので, この式は $y_{t}$ が MA で表現できること, ひいては定常的であることを示している. AR が MA に変換されるだけでなく, 逆に $\theta(\mathrm{L})$ の逆があるならば、MA を AR に変換することもできる。ここから, これらのラグ多項式の逆が存在する条件を, 「反転性条件」という. 特に AR から MA へ変換できる場合, 定常性も同時に証明できるため, AR から MA への反転性条件は時系列が定常的である条件でもある. そして, この定常性条件, つまり $\psi(\mathrm{L})$ の逆数が存在する条件が,「特性方程式

$\begin{align} \psi(z)= & 1-\psi_{1}z-\psi_{2}z^{2}-\cdots-\psi_{p}z^{p}=0\end{align}$

の $z$ の全て解が1より大きい」である*6. これが AR(1) なら, もっと簡単に $|1-\psi_{1}z|=0$ となり, $|\psi_{1}|<1$ が定常性の条件になる $\psi_{1}$ は過去の値に掛ける係数だから, 1より小さければ過去の異常値は時間経過で減衰していくということで, 定常性が直感的にイメージできる. また, 特性方程式は逆に, 「

$\begin{align} z^{p}-\psi_{1}z^{p-1}-\cdots\psi_{p-1}z-\phi_{p}= & 0\end{align}$

の全ての解が1より小さい」とも表現できる. とはいえ, 実際に推定する際には安定性を確認するために特性方程式を解くことはしない*7. ここでは単に, 「定常であるかどうかには明確な条件がある」とだけ分かればよい.

この話を VARMA の場合にも当てはめる. VMA ( $\infty$ ) もまた定常的なので, 同様に, 行列化したラグ多項式 $\boldsymbol{\Psi}(\mathrm{L})$ と $\boldsymbol{\theta}(\mathrm{L})$ を使うと,

$\begin{align} \boldsymbol{\Psi}(\mathrm{L})(\boldsymbol{y}_{t}-\boldsymbol{\mu})= & \boldsymbol{\Theta}(\mathrm{L})\boldsymbol{\varepsilon}_{t}\end{align}$

となる. よって逆行列 $\boldsymbol{\Psi}(\mathrm{L})^{-1}$ が存在すれば, VARMA が VMAと等しいことになるので, ARMA の場合と同様に, 「特性方程式

$\begin{align} \left|\mathbf{I}_{n}-\boldsymbol{\Psi}_{1}z^{p-1}-\boldsymbol{\Psi}_{2}z^{p-2}-\cdots-\boldsymbol{\Psi}_{p}\boldsymbol{z}^{p}\right|= & 0\end{align}$

の全ての解の絶対値が1より大きいかどうか」が定常性条件になる.

推定方法

VAR は過去の値に影響して決まるモデルなので, SUR の一種とみなすこともできる*8. しかし, VARMA の場合は MA のせいで系列相関が発生するため, システム GMM や最尤法を利用する. そこで今回は, 分布の変更も容易な stan で最尤法をベースにコーディングしてみた. といっても, Stan のマニュアルには MA(2), ARMA (1, 1) のコーディング例があるのでそれをもとに作ることができる. しかも stan はベクトルでの記述ができるため, VARMA への拡張はかなり簡単である.ただし, ARMA の例では誤差項を model ブロックのローカル変数に代入してるだけなので, このままでは予測値の計算に使えない. 予測値を求めるには過去の誤差項の列をパラメータとして扱う必要があるので, MA(2) の例をベースに作ることになる.

ARMA は過去の値に依存するため, 尤度が過去の値が所与の条件付き密度関数の積になる (これも厳密な説明は教科書を参照). つまり, $T$ 期間あるデータに対して, 尤度は

$\begin{align} \mathcal{L}(\psi,\,\theta)= & f(y_{T}|y_{T-1},\psi,\,\theta)\times f(y_{T-1}|y_{T-2}, \cdots ) \times \cdots \times f(y_{2}|y_{1},\cdots ) \times f(y_{1}) \end{align}$

というふうに表される. そのため, stan のコードとしては

$\begin{align} y_{t}\sim & \mathcal{N}(\psi y_{t-1} + \theta \varepsilon_{t-1} + \varepsilon_{t},\sigma^{2})\end{align}$

のような再帰的な式を for ループで繰り返すことで尤度を表現できる.

ARMA(1,1)

Phi は stan では累積分布を計算する関数名として予約されているのでパラメータの名前には使えない (phi や PHI は可能). VARMA (1,1) の場合, プログラムは以下のようにした.

gist8148a97d7b5f6e1656d5

事前分布については, ARMA で平均 0, 分散 2 の正規分布を $\psi$ , $\theta$ の事前分布としていたので, VARMA においても, 互いに独立な $\mathcal{N}(0,\,2)$ とした. 平均値も互いに独立な $\mathcal{N}(0,\,10)$ とした. 一方で, ARMA では $\varepsilon_{t}$ の分散をコーシー分布にしているのに対し, 共分散行列は対応する分布がないため, 共分散行列の事前分布によく使われる逆ウィシャート分布で代用した (32行目). いずれにせよ, 事前分布は 29~32 行目の記述を変更することで好きなように設定できる. VARMA バージョンは, 時系列の数が変数 N で与えられているため, N=1 とすれば ARMA にも対応できる. 42行目以降は $y_{t}$ の予測値を求めるプログラムである. 予測期間は T_forecast という変数として入力データに含める.

そして以下が R 上で stan を実行するコードである. forecast パッケージの auto.arima() と arima() の結果とも比較できるようにした. なお, MCMC は次元が大きくなると計算に時間がかかることに留意.

gist9e9eb0ff9a1dd53533bf

やってみるとわかるが, auto.arima() においても次数が正しく識別されることは少ない.

f:id:ill-identified:20160214193647p:plain

予測値のプロットは, 現状は少々面倒で, stan で生成したサンプルを取り出して R 側で予測区間などを計算してグラフにする必要がある. グラフ作成用プログラムは
http://heartruptcy.blog.fc2.com/blog-entry-141.html を参考にした. forecast パッケージのデフォルト値にあわせて, 80 % と 90% の区間を表示している.

2次元のVARMA(1,1)

forecast::arima.sim() は 1変量の時系列しか扱えないため, 次の2次元の VARMA (1,1) は mvrnorm() 関数を使って乱数列を作成している. 一方で, 推定については, R では vars というパッケージがあって, これは VAR の計算ができるが, VARMA は計算できないため, 今度は stan の結果だけを示す. これは, 実際には計算の簡単さから VARMA より VAR がよく好まれるということである*9. そもそも, 単なるガウシアン (V)ARMA ならば MCMC よりニュートン法のほうがはるかに高速なので, 今回のプログラムをそのまま使う実用上のメリットはない. MCMC である意味が発生するのは, もっと複雑な分布を使いたいときだろう.

gist9ab2f579189ff347cc95

f:id:ill-identified:20160214193912p:plain