時系列解析初心者が始める時系列解析(7)
ARモデルの推定
参考：https://elf-c.he.u-tokyo.ac.jp/courses/385
★：著書「経済・ファイナンスデータの計量時系列分析」に記載がある内容を表している。
ARモデルによる予測★
自己回帰モデル（ARモデル）
- ARモデルのおさらいとして、数式とそれぞれの文字が表す内容について、整理する。
\begin{align} y_n = \sum_{j=1}^{m}{a_jy_{n-j}} + v_n \end{align}
| 名称 | 数式 | 備考 |
|---|---|---|
| 定常時系列 | - | |
| ホワイトノイズ |
|
|
| 自己回帰の次数 | - | |
| 自己回帰係数 | - |
-
また、その時のスペクトルは下記のようになる。（分母にフーリエ変換が施されrている。）
p(f) = \frac{\sigma^2}{|1-\sum_{j=1}^{m}{a_j\exp({-2\pi{ijf}}})|^2}~, ~~~~~ -\frac{1}{2} ≤ f ≤ \frac{1}{2} -
この時、次数によって、推定結果は異なる。そのため、複数の次数を仮定して、検証する必要がある。ちなみに、次数が増えると複雑化していく。
-
複数の次数を用いて、実行していくわけだが、どの次数が適しているのかは、対数尤度を見るだけだとわからない（次数を増やすと常に最小値を取り、更新し続けていく性質があるため）。
-
そこで判断基準として、情報量基準（今回はAIC）を用いて、モデルの適切な次数を判断していくわけである。
ARモデルにおける予測
-
(1)
y_n = a_1y_{n-1} + \cdots + a_my_{n-m} + v_n -
n n-1 \backsim n-m m -
その時の予測値を
y_{n|n-1} \varepsilon_n = v_n -
では、１期先の予測
y_{n+1} n n+1
y_{n+1} = a_1y_{n} + \cdots + a_my_{n+1-m} + v_{n+1} -
また、時刻
n y_{n+1} y_{n+1|n}
y_{n+1|n} \equiv a_1y_{n} + \cdots + a_my_{n+1-m} -
その際の予測誤差については以下の通り。
\varepsilon_{n+1} = y_{n+1} - y_{n+1|n} = v_{n+1} -
y_{n+1} y_{n+1|n} - 平均：
E[\varepsilon_{n+1}] = E[v_{n+1}] = 0 - 分散：
E[\varepsilon^2_{n+1}] = E[v^2_{n+1}] = \sigma^2
- 平均：
長期予測の誤差
ARモデルによる長期予測
-
２期以上の予測（長期予測）の場合、どういった形で表せるかを見てみる。
-
まず、2期先
n+2 n+1
\begin{align} y_{n+2} = a_1y_{n+1} + a_2y_n + \cdots + a_my_{n+2-m} + v_{n+2} \end{align} -
その時の時刻
n y_{n+2} n |n
-
この時、以下の性質が成り立つ。
y_{n|n} = y_n y_{n+2-m|n} = y_{n+2-m} v_{n+2|n} = 0
-
以上の性質が成り立つと、少し式を整理することができる。ポイントは第１項の部分で、この部分だけ、時刻
n
\begin{align} y_{n+2|n} = a_1y_{n+1|n} + a_2y_{n} + \cdots + a_my_{n+2-m} \end{align} -
また、3期先の場合も同様に式を構築することができ、時刻
n
y_{n+3|n} = a_1y_{n+2|n} + a_2y_{n+1|n} + \cdots + a_my_{n+3-m}
長期予測の誤差
-
前出の通り、1期先
n+1
\varepsilon_{n+1} = y_{n+1} - y_{n+1|n} = v_{n+1} -
では、2期先
n+2
\begin{align}\varepsilon_{n+2} &= y_{n+2} - y_{n+2|n}\\ &= a_1(y_{n+1} - y_{n+1|n}) + v_{n+2}\\ &= a_1\varepsilon_{n+1} + v_{n+2} \end{align} -
この時の期待値の性質は、式展開すると、0であり、これは一定。（第1項については、
v_{n+1}
E[\varepsilon_{n+2}] = a_1E[\varepsilon_{n+1}] + E[v_{n+2}] = 0 -
分散の性質については、1期先予測時とは少々異なる。より先を予測するということになるため、イメージ通り、誤差（ばらつき）が大きくなる。
E[\varepsilon^2_{n+2}] = a^2_1E[\varepsilon^2_{n+1}] + 2a_1E[\varepsilon_{n+1}v_{n+2}] + E[v^2_{n+2}] = (1+a^2_1)\sigma^2
1変量ARモデルの推定
何故、一般性のあるARMAモデルではなく、ARモデルで推定を行なっていくのか
-
一般に、ARMAモデルの方が汎用性が高いが、その分推定する難易度も上がってしまうという特徴がある。状態空間モデルまで学習したのち、学習していく。
-
ARモデルのメリットとして、ARMAモデルの推定（最尤法等）に比べて、ARモデルの推定（最小二乗法等）の方が簡単であるという点がある。そういった性質から、実用性も高いと評価されている。
ARモデルの同定問題
(関数)同定問題：数式空間を探索する回帰分析の一つで、与えられたデータセットに対して、正確かつ単純な最も相応しいモデル(関数)を見つける問題のこと。
-
ARモデルの数式およびホワイトノイズの仮定。
\begin{align} y_n = \sum_{j=1}^{m}{a_jy_{n-j}} + v_n \end{align}, ~~~~~ v_n \backsim N(0,\sigma^2) -
これらの中で、推定の対象となるのは、AR係数である
a_j~~(j=1,\dots,m) \sigma^2 -
また、次数
m -
これらの推定対象をどのように解いていくかというところで、計4つの手法を以降で紹介していく。
- Yule-Walker法
- 最尤法
- 最小二乗法
- PARCOR法(３種類)
Yule-Walker法とLevinson Algorithm★
Yule-Walker法で解くというのは？
-
最もオーソドックスな方法。数式で書くと以下の通り。（
C_k
\begin{align}C_0 &= \sum_{j=1}^m{a_jC_j} + \sigma^2\\ C_k &= \sum_{j=1}^m{a_jC_{j-k}}~~~~~~~(k = 1,2,\dots) \end{align} -
この時、データが与えられると、
C_k a_j \sigma^2 -
手順としては、まず、(9)の式を解いて、
a_j - 第１項にあたる
\hat{C}_k
\begin{bmatrix}\hat{C}_0 & \hat{C}_1 & \cdots & \hat{C}_{m-1}\\ \hat{C}_1 & \hat{C}_0 & \cdots & \hat{C}_{m-2}\\ \vdots & \vdots & \ddots & \vdots \\ \hat{C}_{m-1} & \hat{C}_{m-2} & \cdots & \hat{C}_{0}\end{bmatrix} \begin{bmatrix}a_1\\ a_2\\ \vdots\\ a_m\end{bmatrix} = \begin{bmatrix}\hat{C}_1\\ \hat{C}_2\\ \vdots\\ \hat{C}_m\end{bmatrix}
- 第１項にあたる
-
その後、(8)の式を分散
\sigma^2 \sigma^2
\hat{\sigma}^2 = \hat{C}_0 - \sum_{j=1}^m{\hat{a}_j\hat{C}_j}
Yule-Walker法では何をやっているのか？
-
Yule-Walker法では、「予測誤差分散の期待値を最小にするように係数
a_j -
そもそも、ARモデルの式の分散
\sigma^2
\sigma^2 = E[v^2_n] -
その時、
v_n n y_n n-j
E[v^2_n] = E[(y_n - \sum_{j=1}^m{a_jy_{n-j}})^2] -
この式を式展開すると以下のようになる。
E[(y_n - \sum_{j=1}^m{a_jy_{n-j}})^2] = E[y_ny_n] - 2\sum_{j=1}^m{a_jE[y_ny_{n-j}]} + \sum_{j=1}^m\sum_{i=1}^m{a_ja_i E[y_{n-j}y_{n-i}]} -
それぞれの項に含まれる期待値については、つまりは自己共分散関数であるので、さらに以下のように整理する。
E[y_ny_n] - 2\sum_{j=1}^m{a_jE[y_ny_{n-j}]} + \sum_{j=1}^m\sum_{i=1}^m{a_ja_i E[y_{n-j}y_{n-i}]}\\ = C_0 - 2\sum_{j=1}^m{a_jC_j} + \sum_{j=1}^m\sum_{i=1}^m{a_ja_i C_{i-j}} -
整理した式について、係数
a_j 0 a_j
\frac{\partial{\sigma^2}}{\partial{a_j}} = -2C_j + 2\sum_{i=1}^m{a_iC_{i-j}} = 0, ~~~~~(j = 1,\dots,m)
Levinson Algorithm
-
Yule-Walker方程式を特にあたって、次元数
m m -
Levinson(レヴィンソン)のアルゴリズムは、その問題を解決する、計算量を抑えてYule-Walker方程式ためのアルゴリズムである。
-
実際、素直に
M M^4/12 2M^2 M=100 \frac{2M^2}{M^4/12} = \frac{24}{M^2} = \frac{24}{10000} = 0.0024 -
係数、分散、
AIC m=0 m=1,\dots,M -
m=0
\begin{align}\sigma^2_0 &= C_0\\ AIC_0 &= N(\log{2\pi \hat{\sigma}^2_0} + 1) + 2 \end{align} -
m=1,\dots,M
\begin{align}(a) ~~~~~ a_m^m &= (C_m - \sum_{j=1}^{m-1}{a_j^{m-1}C_{m-j}})(\sigma_{m-1}^2)^{-1}\\ (b) ~~~~~ a_j^m &= a_j^{m-1} - a_m^m a_{m-j}^{m-1}, ~~~~~ (j=1,\dots,m-1)\\ (c) ~~~~~ \sigma^2_m &= \sigma^2_{m-1}\{{1 - (a_m^m)^2}\}\\ (d) ~~~~~ AIC_m &= N(\log{2\pi\hat\sigma^2_m + 1}) + 2(m+1)\end{align} -
-
Levinson Algorithm の導出については、割愛させてください。
最尤法★
定義通りに解くならば
-
パラメータ
\theta=(a_1,\dots,a_m,\sigma^2)^T -
y \sim N(\mu,C) y \mu C - ちなみに、平均
\mu 0 C - 式展開した先は複雑に見えるが、結局は正規分布の確率密度関数の公式に平均と分散を当てはめているだけ。
- このあとは、対数尤度として計算し、最尤推定値を求めていく。
L(\theta) = p(y_1,\dots,y_N|\theta) = (2\pi)^{-\frac{N}{2}} |C|^{-\frac{1}{2}}\exp{\bigg\{-\frac{1}{2}(x-\mu)^T C^{-1} (x-\mu) \bigg\}}
- ちなみに、平均
各時刻の条件付き分布に分解する方法
-
L(\theta) = p(y_1,\dots,y_N|\theta) \sigma^2 -
どのように分解していくかは以下の通り。
- 式(17)のように、
y_1 \theta y_2 \theta y_1
\begin{align} L(\theta) &= p(y_1,\dots,y_N|\theta)\\ &= p(y_1|\theta)p(y_2,\dots,y_N|y_1,\theta)\\ &= p(y_1|\theta) p(y_2|y_1,\theta) p(y_3,\dots,y_N|y_1,y_2,\theta)\\ &=\cdots\\ &= \Pi_{n=1}^N{p(y_n|y_1,\cdots,y_{n-1},\theta)}\end{align}
- 式(17)のように、
-
式(20)の形にしたことで、正規分布の確率密度関数を用いる時の分散は
\sigma^2 y_n - ただし、下記のように表現できるのは、条件として、
n>m -
n≤m n≤m
{p(y_n|y_1,\cdots,y_{n-1},\theta)} = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\bigg\{-\frac{1}{2\sigma^2} \bigg(y_n - \sum_{j=1}^m{a_jy_{n-j}}\bigg)^2\bigg\} - ただし、下記のように表現できるのは、条件として、
Householder変換に基づいた最小二乗法
最小二乗法で、尤度の近似を行う
-
\sum_{n=1}^N 1 \sim N 1 \sim m -
数式で表すと以下の通り。
\begin{align}l(\theta) &= \sum_{n=1}^N{\log{p(y_n|y_1,\dots,y_{n-1})}}\\ &\cong \sum_{n=m+1}^N{\log{p(y_n|y_{n-m},\dots,y_{n-1})}}\\ &= -\frac{N-m}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{n=m+1}^N (y_n -\sum_{j=1}^m{a_jy_{n-j}})^2 \end{align} -
注意点として、情報量基準でモデル比較をする際は、無視したデータ点の数は全て均一にする必要がある。ここがブレると情報量基準が正しい軸にならなくなってしまう。
近似最尤法〜最小二乗法
-
式(23)については、尤度関数を近似したものを式展開したというふうに捉えることができる。この式をもとに、最小二乗法を用いて、パラメータ
\theta=(\sigma^2,a_j)
l'(\sigma^2,a_j)=-\frac{N-m}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2} \sum_{n=m+1}^N (y_n -\sum_{j=1}^m{a_jy_{n-j}})^2 -
手順１として、
\sigma^2 =0
\frac{\partial{l'(\sigma^2,a_j)}}{\partial{\sigma^2}}=-\frac{N-m}{2\sigma^2} + \frac{1}{2(\sigma^2)^2} \sum_{n=m+1}^N (y_n -\sum_{j=1}^m{a_jy_{n-j}})^2 = 0 -
手順１で得られた式について、
\hat\sigma^2=
\hat\sigma^2 = \frac{1}{N-m}\sum_{n=m+1}^{N}{(y_n -\sum_{j=1}^m{a_jy_{n-j}})^2} -
これをもとに、
l'(\hat\sigma^2, a_j) - ポイントとして、第２項については、元々
\sigma^2 - 第１項については、代入すると式がゴチャゴチャするので、そのままにしていると想定。（特に言及がなかった）
l'(\hat\sigma^2, a_j) = -\frac{N-m}{2} \log{(2\pi{\hat\sigma^2})} - \frac{N-m}{2}
- ポイントとして、第２項については、元々
-
以上より、
l'(\hat\sigma^2, a_j) a_j \hat\sigma^2
\max_{a_j}{l'(\hat\sigma^2, a_j)} \lrArr \min_{a_j}{\hat\sigma^2} \lrArr \min_{a_j}{\sum_{n=m+1}^N{(y_n - \sum_{j=1}^m{a_jy_{n-j}})}}
Householder法に置き換える
-
Householder法の詳細については、下記URLを参照してください。
https://zenn.dev/udai_sharelearn/articles/time_series_estimation_choice#モデルの推定・選択 -
最小二乗法を素直に解くとするならば、
y=Za+v v a -
y,a,v Z
||v||^2 = ||y-Za||^2
\hat{a} = (Z^TZ)^{-1}Z^Ty \rightarrow これを解く
-
-
ただ、この解き方よりも優れた解き方というのがあり、それがHouseholder法である。この方法では、共分散行列を使用せずにデータ行列から直接推定するため、トータルの計算量が抑えられるというメリットがある。
-
それ以外にも 「精度が２倍に向上」 する点や、「モデルの自由度（変数選択、次数選択）」 が高い点、「計算上の便利さ（データが増えた場合でも柔軟に対応ができる）」 という点がある。
PARCOR法
PARCOR法
-
Yule-Walker法の方程式をLevinsonアルゴリズムで解く際、以下の式にあるように、自己共分散関数を用いる必要があった。
\hat a_m^m = (\hat C_m - \sum_{j=1}^{m-1}{\hat a_j^{m-1}\hat C_{m-j}})(\hat\sigma_{m-1}^2)^{-1} -
PARCOR法では、
\hat a_m^m -
考え方としては、
m-1 m-1 -
数式的に考えるとまず、後向きARモデルを定義する必要がある。どこがポイントかというと、
n-m n
y_{n-m} = \sum_{i=1}^{m-1}a_j^{m-1}y_{n-m+j} + w_{n-m}^{m-1} -
実際に式変形をしていき、自己共分散関数
C - 式(24)については、自己共分散関数を期待値の形の表現に戻していると捉えると良い。
- 式(25)については、
\bigg(y_n - \sum_{i=1}^{m-1}{a_j^{m-1}y_{n-j}}\bigg) v_n^{m-1} - 式(26)が少々複雑。シンプルに見ると
y_{n-m} w_{n-m}^{m-1} E[v_n^{m-1}y_{n-m-j}] = 0 ~~~~ (j= 1,\dots,m-1) y_{n-m-j} 0 y_{n-m} = w_{n-m}^{m-1}
PARCOR法での予測誤差分散
-
PARCOR法を適応した際の予測誤差分散についても定義をしていくが、考え方は先ほどと同様で、後向きモデルのノイズ
w -
まず、予測誤差分散の式について確認する。今回の式変形の対象である。これは、Yule-Walker方程式の際に一度出てきている。
\hat\sigma_{m-1}^2 = C_0 - \sum_{j=1}^{m-1}{a_j^{m-1}C_j} -
では、この式の右辺を先ほどと似た要領で変形していく。
- 式(28)については、式(24)と同様で、自己共分散関数を期待値の形に変形している。
- 式(29)については、
\bigg(y_{n-m} - \sum_{i=1}^{m-1}{a_j^{m-1}y_{n-m+j}}\bigg) w_{n-m}^{m-1} - 式(30)は、式(26)での説明の通りのため、それを行うと、２乗の形になるということである。
\begin{align}C_0 - \sum_{j=1}^{m-1}{a_j^{m-1}C_j} &= E\bigg\{\bigg(y_{n-m} - \sum_{i=1}^{m-1}{a_j^{m-1}y_{n-m+j}}\bigg)y_{n-m}\bigg\}\\ &= E(w_{n-m}^{m-1}y_{n-m})\\ &= E(w_{n-m}^{m-1})^2 \end{align}
-
また、
E(w_{n-m}^{m-1})^2 = E(v_n^{m-1})^2
E(w_{n-m}^{m-1})^2 \cong \begin{cases}\frac{1}{N-m}\sum_{n=m+1}^N(w_{n-m}^{m-1})^2\\~\\ \frac{1}{N-m}\bigg\{\sum_{n=m+1}^N(w_{n-m}^{m-1})^2 \sum_{n=m+1}^N(v_{n}^{m-1})^2\bigg\}^{\frac{1}{2}}\\~\\ \frac{1}{2(N-m)}\bigg\{\sum_{n=m+1}^N(w_{n-m}^{m-1})^2 + \sum_{n=m+1}^N(v_{n}^{m-1})^2\bigg\} \end{cases}
各推定法の特徴と数値例
4つの手法を各５項目で評価
- それぞれの手法を紹介したわけだが、ここで実際にどういった特徴があるのかを計５項目で評価してみる。
| 定常性 | 数値精度 | モデル精度 | 速度 | 自由度 | |
|---|---|---|---|---|---|
| Yule-Walker法 | ◯ | △ | △ | ◯ | ◯ |
| 最尤法 | ◎ | ◯ | ◎ | × | ◯ |
| 最小二乗法 | × | ◎ | ◯ | ◎ | ◎ |
| PARCOR法 | ◎ | ◯ | ◯ | ◎ | △ |
- どの面でも優れているという手法はない。実際やってみての判断になる。Rでもそれがパラメータとして、用意されている。（おそらくPythonでも。未確認。）
多変量ARモデルの推定(Levinson-Durvin法、Householder法)
Yule-Walker法
-
多変量ARモデルの場合はどうなるのかというところを説明していく。
-
多変量ARモデルの数式については下記の通り。
y_n = \sum_{i=1}^m{A^m_jy_{n-j}} + v_n, ~~~~~~~ v_n \sim N(0,V_m) -
この時自己共分散関数
C_k = E[y_ny_{n-k}^T]
\begin{align} C_0 &= \sum_{i=1}^m{A_jC_{-j}} + V\\ C_k &= \sum_{i=1}^m{A_jC_{k-j}} ~~~~~~ (k=1,2,\dots) \end{align} -
式(32)について、行列を用いて表現する。
- この時、自己共分散関数の行列はテンプリッツ行列ではあるものの、対称ではないことに注意が必要。
\begin{bmatrix}{C}_0 & {C}_{-1} & \cdots & C_{1-m}\\ C_1 & C_0 & \cdots & C_{2-m}\\ \vdots & \vdots & \ddots & \vdots \\ C_{m-1} & C_{m-2} & \cdots & C_{0}\end{bmatrix} \begin{bmatrix}A_1\\ A_2\\ \vdots\\ A_m\end{bmatrix} = \begin{bmatrix}C_1\\ C_2\\ \vdots\\ C_m \end{bmatrix}
- この時、自己共分散関数の行列はテンプリッツ行列ではあるものの、対称ではないことに注意が必要。
-
この方程式の解き方だが、1変量の時と同様、Levinsonアルゴリズムを用いることができる。
-
ただ、1変量と異なる点として、前向きのARモデルと後向きのARモデルが異なっている。これまでは同一という性質を使って、式変形等を行ってきたので、この違いは大きい。推定するパラメータも増えてしまう。
- 前向きモデル：
y_n = \sum_{i=1}^m{A^m_jy_{n-j} + v_n}, ~~~~~~~ v_n \sim N(0,V_n) - 後向きモデル：
y_n = \sum_{i=1}^m{B^m_jy_{n+j} + u_n}, ~~~~~~~ u_n \sim N(0,U_n) - 推定するパラメータ：
\{ ((A_j^m, B_j^m), ~~ j=1,\dots,m), V_m, U_m \}, ~~~~~ m=1,\dots,M
- 前向きモデル：
-
それぞれの計算式については、割愛。気になる人は動画を確認してしてください。
最小二乗法
-
最小二乗法については、1変量の時のような解き方をしてしまうと、多変量である分、計算量が膨大になってしまう。
y_n = \sum_{i=1}^m{A_jy_{n-j}} + v_n, ~~~~~~~ v_n \sim N(0,V) -
そういった問題をクリアするために、工夫として、同時応答モデルというのを作成し、計算量を抑える、
y_n = B_0y_n + \sum_{i=1}^m{B_jy_{n-j}} + w_n, ~~~~~~~ w_n \sim N(0,W) -
この式の中身である、
B_0 W
B_0 = \begin{bmatrix} 0 & 0 & \cdots & 0 \\~\\ b_0(2,1) & 0 & \cdots & 0\\ ~ \\ \vdots & \ddots & \ddots & \vdots\\~\\ b_0(k,1) & \cdots & b_0(k,k-1) & 0 \end{bmatrix},~~~~~ W = \begin{bmatrix} \sigma_1^2 & & & 0 \\~\\ & \sigma_2^2 & & \\ ~ \\ & & \ddots & \\~\\ 0 & & & \sigma_k^2 \end{bmatrix} -
先ほど定義した式について、式変形を行い、
A_j V - ここで行っていることは、元の式の
B_0y_n y_n (I-B_0) -
I
y_n = (I-B_0)^{-1}\sum_{i=1}^m{B_jy_{n-j}} + (I-B_0)^{-1}w_n
- ここで行っていることは、元の式の
-
式変形した結果を用いて、
A_j V A_j = (I-B_0)^{-1}B_j V = (I-B_0)^{-1}W(I-B_0)^{-T}
-
ここで変換したことで得られるメリットは、
W -
また、単純に計算量が減少する点や、モデリング時、変数（説明変数、目的変数）ごとに異なるが決められる点がメリットとしてある。様々な条件下で推定を行う上でこのメリットは大きいと言える。
変数選択の方法と選択例
変数選択の方法
-
重回帰モデルでは、あるデータに対してのモデル構築する際、不要な変数については、取り除くことが可能であった。
-
しかし、多変量時系列モデルの場合はそうはいかない。目的変数＝説明変数という関係があるために、説明変数を切り落としてしまうと、対応する目的変数も切り落としてしまう。すると、各条件下で構築されたモデルを比較して評価することができなくなってしまう。
-
そこで、対策として、2つ以上のグループで切り分けてモデルを構築し、その時のAICを足し算することで、比較して評価するという方法をとるらしい。
-
2つ以上のグループに切り分ける際に、Aグループで不要だった説明変数を、Bグループとして独立してモデリングするということである。
-
言葉での説明は難しいので、動画での具体例を見るとわかりやすいかと思う。
Discussion