時系列解析初心者が始める時系列解析(6)
ARMAモデルによる時系列の解析
参考：https://elf-c.he.u-tokyo.ac.jp/courses/384
1変量ARMAモデル
時系列モデルを考えるにあたって
- 時系列モデルでは、ホワイトノイズが重要な役割を占める。
- 過去データである
x_{n-1} y_n - その表現しきれない部分をホワイトノイズとする。
自己回帰移動平均(ARMA)モデル
- ARMAモデルでは、目的変数である
y_n y_{n-j} v_n v_{n-j} - 数式にすると以下の通りである。
y_n = \sum_{j=1}^m{a_jy_{n-j}} + v_n - \sum_{j=1}^l{b_jv_{n-j}}
| 名称 | 数式 | 名称 | 数式 |
|---|---|---|---|
| 1変量時系列 | ホワイトノイズ | ||
| AR次数 | MA次数 | ||
| AR係数 | MA係数 |
-
この時のホワイトノイズの性質は以下の通り。
v_n \backsim N(0,\sigma^2)~~~正規分布
E[v_nv_k]=0 ~~~ n≠k ~~~ 無相関
E[v_ny_{n-j}]=0 ~~~ j>0 ~~~ 過去の時系列と無相関 -
以上のように、定義されるモデルのことをARMA(AutoRegressive Moving Average)モデルという。表現として、ARMA(
l,m
ARMAモデルの特殊な場合
-
l=0 m
y_n = \sum_{j=1}^m{a_jy_{n-j}} + v_n -
m=0 l
y_n = v_n - \sum_{j=1}^l{b_jv_{n-j}} -
m=l=0
y_n = v_n -
m=1,l=0 - マルコフ過程：１期前の値だけに依存する確率過程のこと。
y_n = a_1y_{n-1} + v_n
- マルコフ過程：１期前の値だけに依存する確率過程のこと。
-
m=1,l=0,a_1=1 - ランダムウォーク：確率
p １ 1-p -1
y_n = y_{n-1} + v_n
- ランダムウォーク：確率
ARIMAモデル
- ARIMAモデルとは、平均に関して非定常なデータに対して、階差をとることで定常データに近づけて、モデリングする手法である。
-
y_n - その時、１時点前との差分を
\Delta{y}_n = y_n - y_{n-1} - それを任意の
k z_n = \Delta^k{y_n} - 以上を踏まえると、以下のように式変形できる。
z_n = \sum_{j=1}^m{a_jz_{n-j}} + v_n - \sum_{j=1}^l{b_jv_{n-j}}
-
SARIMAモデル
- ARIMAモデルに、季節周期性の成分を加えたもの。
-
p \Delta_p = 1- B^p - それを用いて、周期性を表した
z_n z_n = \Delta{_p}y_n = y_n - y_{n-p} y_n \Delta_p - 以上を踏まえると、以下のように式変形できるが、見た目は、ARIMAモデルと全く同じである。
z_n = \sum_{j=1}^m{a_jz_{n-j}} + v_n - \sum_{j=1}^l{b_jv_{n-j}}
-
インパルス応答関数
Back-shift operatorによる数式の表現
-
B L -
By_n \equiv y_{n-1} y_n B y_{n-1} -
B^ky_n = y_{n-k} k
-
-
以上を用いて、ARMAモデルの数式当てはめると以下の式が導出される。（元の式の右辺第１項を移行している。）
\bigg({1 - \sum_{j=1}^m{a_jB^j}}\bigg)y_n = \bigg({1 - \sum_{j=1}^l{b_jB^j}}\bigg)v_n -
この式より、ARオペレータとMAオペレータを定義する。
a(B) = {1 - \sum_{j=1}^m{a_jB^j}} ~~~ ARオペレータ
b(B) = {1 - \sum_{j=1}^l{b_jB^j}} ~~~ MAオペレータ -
すると、より簡単に表すことができる。
a(B)y_n = b(B)v_n
インパルス応答関数 g_j
-
a(B)y_n = b(B)v_n y_n=
y_n = a(B)^{-1}b(B)v_n -
ここで、
a(B)^{-1}b(B) g
g(B) \equiv a(B)^{-1}b(B) = \sum_{j=0}^\infin{g_jB^j},~~~~~g_1 = 1 -
よって、
g(B) - (2)→(3)のところでは、
B^jv_n
\begin{align}y_n &= g(B)v_n\\ &= \sum_{j=0}^\infin{g_jB^jv_n}\\ &= \sum_{j=0}^\infin{g_jv_{n-j}} \end{align}
- (2)→(3)のところでは、
-
ここでいう
g_j g_j:j=0,1,... -
インパルス応答関数が何を表しているのかというと、ある事象がデータのどういった影響を与えたか、である。イメージとしては、最初に大きな衝撃(インパルス)を与えて、その後放っておくとどのように変化し落ち着くのかをみるというようなものである。
自己共分散関数
自己共分散関数をインパルス応答関数を用いて表す。
-
自己共分散関数
C_k g_j
\begin{align}C_0 &= \sum_{j=1}^m{a_jC_j} + \sigma^2(1 - \sum_{j=1}^l{b_jg_j})\\ C_k &= \sum_{j=1}^m{a_jC_{k-j}} - \sigma^2\sum_{j=1}^l{b_jg_{j-k}} \end{align} -
この式は覚えるしかないものだが、この式からわかることとしては、パラメータである
m,l,a_j,b_j,\sigma^2 g_j g_j C_0,C_1,\cdots
偏自己相関係数とAR係数の関係
偏自己相関係数(PARCOR)
-
偏自己相関係数：時系列データにおいて、
y_n y_{n-j} -
共分散でも測れるのでは？となるが、実際に計算していくと0という結果が出てこないとに気づく。0であれば、関係ないと言えるが、それが出ないとなると指標として成立しない。
-
ARモデルの回帰式
y_n = \sum_{j=1}^m{a_j^my_{n-j}} + v_n a_j^m m j -
その中でも、次数
m m b_m
b_m \equiv a_m^m -
偏自己相関係数については、直交化されているので数学的に便利な性質を持っているらしい。
-
また、利用される場面として、音声分析、音声合成、画像認識が挙げられる。
偏自己相関係数b_m a
-
どんな関係があるのかというと、ある偏自己相関係数がわかると、それ以外の係数がわかるということである。式で書くと以下の通り。
a_i^m = a_i^{m-1} - a_m^m{a_{m-i}^{m-1}} ~~~ (i=1,\cdots,m-1) -
例えば、
a_1^2 a_1^1 a_2^2 -
もう一つの関係として、ある次数における係数がすべてわかっていたとすると、その1つ下の次数の係数もすべて導出できるということがある。式で書くと以下の通り。
- 例えば、次数4の時の係数が4つともわかっていた時、次数3の時の係数が3つともわかるということである。
パワースペクトル
-
本記事シリーズの(2)の時に記載した通り、パワースペクトルは、自己共分散
C_k -
本編は数式だらけなため、細かい内容は省略するが、やっていることとしては、パワースペクトル
p(f) g_j
-
結局、スペクトルが表現したいのは、波であるが、AR,MA,ARMAでそれぞれ表す波には特徴がある。
- AR：山を表現するモデル
- MA：谷を表現するモデル
- ARMA：山と谷を表現するモデル
-
上記の性質を数式で表すには、(9)の式の対数を取るとわかる。
\log{p(f)} = \log{\sigma^2} -2\log{\biggm\vert{1 - \sum_{j=1}^m{a_j\exp{(-2\pi{ijf})}}}\biggm\vert} + 2\log{\biggm\vert{1 - \sum_{j=1}^l{b_j\exp{(-2\pi{ijf})}}}\biggm\vert} -
この時、スペクトルの山は
\biggm\vert{1 - \sum_{j=1}^m{a_j\exp{(-2\pi{ijf})}}}\biggm\vert \biggm\vert{1 - \sum_{j=1}^l{b_j\exp{(-2\pi{ijf})}}}\biggm\vert -
また、スペクトルにおいて、
k 2k 2k
特性方程式
ARオペレータa(B) b(B)
-
まず、それぞれの特性方程式を定義する。
a(B) = 1 - \sum_{j=1}^{m}{a_jB^j} = 0
b(B) = 1 - \sum_{j=1}^{l}{b_jB^{j}} = 0 -
この時、それぞれに性質が存在する。
-
a(B) = 0 -
b(B) = 0 - 反転可能：b(B)はMAオペレータであるように、MAモデルで表されているが、これの逆数
b(B)^-1
\begin{align}b(B)^{-1} &= \sum_{i=0}^{\infin}{g_iB^i}\\ y_n &= \sum_{i=1}^{\infin}{g_iy_{n-i}} + v_n \end{align}
- 反転可能：b(B)はMAオペレータであるように、MAモデルで表されているが、これの逆数
-
-
ここで、今後、便利な方程式として利用したいため、
g(B) B B^{-1} \lambda = B^{-1} -
すると、
a(B),b(B)
a(\lambda) = \lambda^m - \sum_{j=1}^{m}{a_j\lambda^{m-j}} = 0
b(\lambda) = \lambda^m - \sum_{j=1}^{l}{b_j\lambda^{l-j}} = 0 -
また、こうすることで性質についても変化が起こる。
-
a(\lambda) = 0 -
b(\lambda) = 0
-
多変量ARモデル
多変量ARモデルを定義する
-
そもそもなぜARMAモデルではないのかというところがあるが、これは、ARモデルと同様に定義ができるという点、ARモデルで色々なことが表現できる点、ARMAモデルでモデリングする場合は、様々な面で負荷が高くなるという理由から、ARモデルで定義することになる。
-
まず、多変量ということなので、ベクトル、行列表現がそれぞれ出てくる。
y_n = \sum_{j=1}^{m}{A_jy_{n-j} + v_n} ~~~~~ 多変量ARモデル -
それぞれの項についての説明は下記の通り。
-
y_n = (y_n(1),\cdots,y_n(l))^T -
A_j = \begin{bmatrix}a_j(1,1) & \cdots & a_j(1,l)\\ \vdots & \ddots & \vdots\\ a_j(l,1) & \cdots & a_j(l,l) \end{bmatrix} -
v_n l -
E[v_n] = \begin{bmatrix}0 \\ \vdots \\ 0s \end{bmatrix} E[v_nv_n^T] = \begin{bmatrix}\sigma_{11} & \cdots & \sigma_{1l}\\ \vdots & \ddots & \vdots\\ \sigma_{l1} & \cdots & \sigma_{ll}\end{bmatrix} = W ~~~~~~ \sigma_{ij} = \sigma_{ji}\rightarrow 対称行列 -
E[v_nv_m^T] = 0 ~~~~~ n≠m v_n -
E[v_ny_m^T] = 0 ~~~~~ n>m v_n y_n
-
-
相互共分散関数
-
相互共分散関数
C_k
C_k = E[y_ny_{n-k}^T]
C_k = \begin{bmatrix}C_k(1,1) & \cdots & C_k(1,l)\\ \vdots & \ddots & \vdots\\ C_k(l,1) & \cdots & C_k(l,l)\end{bmatrix} -
この時、Yule-Walker方程式を用いて、
C_0,C_k
\begin{align}C_0 &= \sum_{j=1}^{m}{A_jC_{-j}+W}\\ C_k &= \sum_{j=1}^{m}{A_jC_{k-j}}~~~~~(k=1,2,\cdots)\end{align}
クロススペクトル
-
相互相関関数をフーリエ変換したもの。そのため、
C_k
P(f) = \sum_{k=-\infin}^{\infin}{C_k\exp(-2\pi{ikf})}
P(f) = \begin{bmatrix}p_{11}(f) & \cdots & p_{1m}(f)\\ \vdots & \ddots & \vdots\\ p_{m1}(f) & \cdots & p_{mm}(f)\end{bmatrix} -
また、
y_n P(f) A_{jk}(f) = I - \sum_{j=1}^{m}{a_j(j,k)\exp{-2\pi{ijf}}} -
()^*
P(f) = A(f)^{-1}W(A(f)^{-1})^*
-
多変量ARの場合、
p(f) c(f) id(f)
p_{jk}(f) = c_{jk}(f) + i{d_{jk}(f)} -
この後、コヒーレンシーを求めることになるため、それを求めるために必要なパーツをいくつか用意する。
- コヒーレンシー：2つの時系列間の相関を周波数領域で見る指標。2つの時系列の相関があれば、1に近づくというもの。
| 名称 | 数式 | 中身 |
|---|---|---|
| 振幅スペクトル | ||
| 位相スペクトル | ||
| クロススペクトル |
フィードバックシステム・パワー寄与率
フィードバックシステム
-
多変量の時系列データの場合、ある変数Aに対して、ある変数Bが影響を及ぼすことがある。逆も然りで、変数Bに対して、変数Aが影響を及ぼすこともある。
-
これをフィードバックシステムという。
- 例１：農業生産物の生産量が上がると価格が下がるが、その後生産量下げることで価格が上がる。
- 例２：石油価格が上がるとエコカーの販売量が上がり、石油価格が下がるとエコカーの販売量が下がる。
-
これは、つまり、インパルス応答関数の形が変わることを意味する。
- 1変量の場合：開ループ。ある周波数に対して、変数Aを通した時のその後の周波数の変化を見るが、それは一方通行である。
- 多変量の場合：閉ループ。変数Aを通すまでは同様だが、その出力を変数Bへ渡す。さらに、変数Bの出力を変数Aへ返すという風に、ループが発生する。
2変量の場合のフィードバックシステムを考えてみる
- まず、
y_n,x_n -
外乱：外部からの干渉というようなイメージ。具体的には、変数Aを通してでた出力から変数Bの入力までの間に、外乱
u_n x_n
-
外乱：外部からの干渉というようなイメージ。具体的には、変数Aを通してでた出力から変数Bの入力までの間に、外乱
| 名称 | 数式 |
|---|---|
| ある変数Aのインパルス応答関数 | |
| ある変数Bのインパルス応答関数 | |
| 互いに無相関な外乱 |
-
x_n,y_n a_j,b_j u_n,v_n -
これを解いていくのであれば、一見、通常の最小二乗法でいい気もする。ただ、そのためには、誤差項の部分(今回でいえば外乱)が独立である必要があった。
-
現状、
x_n u_n -
そこでとるアプローチが外乱
u_n,v_n -
ARモデルを用いてホワイトノイズ化する。数式で書くと以下の通り。
\begin{align}u_n &= \sum_{j=1}^{\infin}{c_ju_{n-j}+\varepsilon_n}\\ v_n &= \sum_{j=1}^{\infin}{d_jv_{n-j}+\delta_n}\end{align} -
ここで、式(14)(16)にて、
u_n
-
さらに、この式を
y_{n-j},x_{n-i} c_i
\begin{align}a_j(1,1) &= c_j\\ a_j(1,2) &= a_j - \sum_{i=1}^{j-1}{c_ia_{j-i}}\end{align} -
式(20),(21)を用いて、式(19)を変換する。
y_n = \sum_{j=1}^{\infin}{a_j(1,1)y_{n-j}} + \sum_{j=1}^{\infin}{a_j(1,2)x_{n-i}} + \varepsilon_n -
また、式(15)についても式(17)を用いてることで、同様な変形ができる。(流れは一緒)
x_n = \sum_{j=1}^{\infin}{a_j(2,1)y_{n-j}} + \sum_{j=1}^{\infin}{a_j(2,2)x_{n-i}} + \delta_n -
以上の式から、「2変量のARモデル」、「ノイズのモデル」、「変数A,Bのインパルス応答関数」の式を定義でき、それを用いて、推定対象である係数
a_j,b_j u_n,v_n
パワー寄与率
- 数式で表すとかなり複雑であった。（動画には記載あり。）
- ニュアンスとしては、ある時系列A,B,Cがあった時に、そのパワースペクトルをノイズごとに分解して、その時系列がどのノイズに起因するかを表すといったもの。
- 主成分分析でもあったように、どの要素（今回でいえばノイズ）がそのデータを説明できているのかというふうに解釈できる。
- この点の解釈は現時点では、不明確なので、正式なものとは捉えないようにしてほしい。
Discussion