時系列解析初心者が始める時系列解析(7)

2023/01/12に公開約15,300字

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}
名称 数式 備考
定常時系列 y_n -
ホワイトノイズ v_n v_n \backsim N(0,\sigma^2)
E[v_nv_k]=0 ~~~~~ n≠k
E[v_ny_{n-j}]=0 ~~~~~ j>0
自己回帰の次数 m -
自己回帰係数 a_j -
  • また、その時のスペクトルは下記のようになる。(分母にフーリエ変換が施され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 と表す。

  • では、1期先の予測 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期以上の予測(長期予測)の場合、どういった形で表せるかを見てみる。

  • まず、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+2|n} = a_1y_{n+1|n} + a_2y_{n|n} + \cdots + a_my_{n+2-m|n} + v_{n+2|n}
  • この時、以下の性質が成り立つ。

    • y_{n|n} = y_n
    • y_{n+2-m|n} = y_{n+2-m}
    • v_{n+2|n} = 0
  • 以上の性質が成り立つと、少し式を整理することができる。ポイントは第1項の部分で、この部分だけ、時刻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が与えられているのは、第1項、第2項まで。これが、期を進めるごとに増えていく。

    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の場合ではどうなるのか。(2)と(3)の式を用いて、式展開してみると以下のようになる。

    \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}に置き換えられるため、0になる。)

    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つの手法を以降で紹介していく。

    1. Yule-Walker法
    2. 最尤法
    3. 最小二乗法
    4. PARCOR法(3種類)

Yule-Walker法とLevinson Algorithm


Yule-Walker法で解くというのは?

  • 最もオーソドックスな方法。数式で書くと以下の通り。(C_kは自己共分散関数→同シリーズの(2)の記事参照)

    \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が求まるので、それを用いて、推定対象であるAR係数a_jと分散\sigma^2を推定していく。

  • 手順としては、まず、(9)の式を解いて、a_jを求める。

    • 第1項にあたる\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について、式展開し、それぞれの項に、(9)の解を代入し、\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の計算量が必要になる。それを、Levinsonアルゴリズムを用いることで、2M^2までに抑えることが可能。

    M=100の時

    \frac{2M^2}{M^4/12} = \frac{24}{M^2} = \frac{24}{10000} = 0.0024

  • 係数、分散、AICの計算については、m=0m=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の時(偏自己相関PARCORを用いる - 式(12))

    \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のみ。式(18)では、y_2も同様に切り出すわけだが、その時の条件は\thetay_1となる。これを繰り返し行なっていくことで、式(20)にたどり着く。
      \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}
  • 式(20)の形にしたことで、正規分布の確率密度関数を用いる時の分散は\sigma^2が使える。(y_nを対象にするため、冒頭から記述しているARモデルの1次式となるため)

    • ただし、下記のように表現できるのは、条件として、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

  • 手順1として、\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

  • 手順1で得られた式について、\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)を作成し、整理する。

    • ポイントとして、第2項については、元々\sigma^2だったところに、上の式を代入して、整理している。
    • 第1項については、代入すると式がゴチャゴチャするので、そのままにしていると想定。(特に言及がなかった)
      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法である。この方法では、共分散行列を使用せずにデータ行列から直接推定するため、トータルの計算量が抑えられるというメリットがある。

  • それ以外にも 「精度が2倍に向上」 する点や、「モデルの自由度(変数選択、次数選択)」 が高い点、「計算上の便利さ(データが増えた場合でも柔軟に対応ができる)」 という点がある。


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次の前向きのARモデルとm-1次の後向きのARモデルの予測誤差の相関係数を計算して求めるということである。

  • 数式的に考えるとまず、後向きARモデルを定義する必要がある。どこがポイントかというと、n-mの部分である。この点は前向きであれば、nなので、この点が後向きARモデルを指していると考えることができる。

    y_{n-m} = \sum_{i=1}^{m-1}a_j^{m-1}y_{n-m+j} + w_{n-m}^{m-1}

  • 実際に式変形をしていき、自己共分散関数Cではなく、前向きARモデルと後向きARモデルの相関係数で表現してみる。

    • 式(24)については、自己共分散関数を期待値の形の表現に戻していると捉えると良い。
    • 式(25)については、\bigg(y_n - \sum_{i=1}^{m-1}{a_j^{m-1}y_{n-j}}\bigg)がノイズを表しているので、前向きARモデルのノイズv_n^{m-1}に置き換わっている。
    • 式(26)が少々複雑。シンプルに見るとy_{n-m}が後向きARモデルのノイズw_{n-m}^{m-1}に置き換わっている。これは、先ほど定義した後向きARモデルの式の性質があるからこそ成り立つ置き換えである。それは、E[v_n^{m-1}y_{n-m-j}] = 0 ~~~~ (j= 1,\dots,m-1)という性質である。後向きARモデルの第1項を見ると、y_{n-m-j}があり、式(26)に後向きARモデルの右辺を代入すると、まさにその第1項が0になる。それが成り立つということは、y_{n-m} = w_{n-m}^{m-1}となるので、置き換えることができる。
\begin{align}C_m - \sum_{j=1}^{m-1}{a_j^{m-1}C_{m-j}} &= E\bigg\{\bigg(y_n - \sum_{i=1}^{m-1}{a_j^{m-1}y_{n-j}}\bigg) y_{n-m}\bigg\}\\ &= E(v_n^{m-1}y_{n-m})\\ &= E(v_n^{m-1} w_{n-m}^{m-1})\\ &\cong \frac{1}{N-m}\sum_{n=m+1}^{1}{(v_n^{m-1}w_{n-m}^{m-1})} \end{align}

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)が後向きARモデルのノイズを表しているので、w_{n-m}^{m-1}に置き換えている。
    • 式(30)は、式(26)での説明の通りのため、それを行うと、2乗の形になるということである。
      \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が成り立つらしくそれを用いると最終の式は3つの形に変形できる。

    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つの手法を各5項目で評価

  • それぞれの手法を紹介したわけだが、ここで実際にどういった特徴があるのかを計5項目で評価してみる。
定常性 数値精度 モデル精度 速度 自由度
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]を用いて、Yule-Walker方程式を作成する。

    \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_jVを置き換えられるようにする。

    • ここで行っていることは、元の式のB_0y_nを左辺に移行して、それをy_nでまとめて出てきた係数(I-B_0)を両辺で割っているということである。
    • Iは単位行列で、行列が出てくる式での1の役割。
      y_n = (I-B_0)^{-1}\sum_{i=1}^m{B_jy_{n-j}} + (I-B_0)^{-1}w_n
  • 式変形した結果を用いて、A_jVを置き換える。

    • 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

ログインするとコメントできます