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

2023/01/10に公開約11,000字

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変量時系列 y_n ホワイトノイズ v_n
AR次数 m MA次数 l
AR係数 a_j MA係数 b_j
  • この時のホワイトノイズの性質は以下の通り。

    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モデルの特殊な場合

  1. l=0の時は、ARモデルとなる。AR(m)と書く。自己回帰モデル

    y_n = \sum_{j=1}^m{a_jy_{n-j}} + v_n

  2. m=0の時は、MAモデルとなる。MA(l)と書く。移動平均モデル

    y_n = v_n - \sum_{j=1}^l{b_jv_{n-j}}

  3. m=l=0の時は、ホワイトノイズを表す。

    y_n = v_n

  4. m=1,l=0の時、マルコフ過程となる。

    • マルコフ過程:1期前の値だけに依存する確率過程のこと。
      y_n = a_1y_{n-1} + v_n
  5. m=1,l=0,a_1=1の時、ランダムウォークという。

    • ランダムウォーク:確率pポイント、確率1-p-1ポイントのゲームを繰り返すような確率モデルのこと。
      y_n = y_{n-1} + v_n

ARIMAモデル

  • ARIMAモデルとは、平均に関して非定常なデータに対して、階差をとることで定常データに近づけて、モデリングする手法である。
    • y_nを(平均)非定常時系列とする。
    • その時、1時点前との差分を\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を1周期の長さとし、\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(Back-shift operator):ラグをとるという意味である。数式で表すと以下の通り。L(Lag operator)と表記されることもある。

    • By_n \equiv y_{n-1}:原系列データy_nBを施すと、y_{n-1}という過去データになる。
    • B^ky_n = y_{n-k}k階差分行った場合の表し方。
  • 以上を用いて、ARMAモデルの数式当てはめると以下の式が導出される。(元の式の右辺第1項を移行している。)

    \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}
  • ここでいう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_ny_{n-j}がどれくらい関係性があるかを見るにあたって、重要な指標である。

  • 共分散でも測れるのでは?となるが、実際に計算していくと0という結果が出てこないとに気づく。0であれば、関係ないと言えるが、それが出ないとなると指標として成立しない。

  • ARモデルの回帰式y_n = \sum_{j=1}^m{a_j^my_{n-j}} + v_nがある時、a_j^mは次数mのARモデルのj番目の係数とみなすことができる。

  • その中でも、次数mm番目の係数については、偏自己相関係数b_m と呼ばれる。

    b_m \equiv a_m^m

  • 偏自己相関係数については、直交化されているので数学的に便利な性質を持っているらしい。

  • また、利用される場面として、音声分析、音声合成、画像認識が挙げられる。

偏自己相関係数b_mとAR係数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^1a_2^2がわかっていれば、求めることができる。

  • もう一つの関係として、ある次数における係数がすべてわかっていたとすると、その1つ下の次数の係数もすべて導出できるということがある。式で書くと以下の通り。

a_i^{m-1} = \frac{a_i^m + a_m^ma_{m-i}^{m-1}}{1-(a_m^m)^2} ~~~~~ (i=1,\cdots,m-1)
  • 例えば、次数4の時の係数が4つともわかっていた時、次数3の時の係数が3つともわかるということである。

パワースペクトル


  • 本記事シリーズの(2)の時に記載した通り、パワースペクトルは、自己共分散C_kをフーリエ変換したものになる。

  • 本編は数式だらけなため、細かい内容は省略するが、やっていることとしては、パワースペクトルp(f)をインパルス応答関数g_jで表すということである。

\begin{align}p(f) &= \sum_{k=-\infin}^\infin{C_k\exp{(-2\pi{ikf})}}\\ &= \sigma^2|\sum_{j=0}^{\infin}{g_j\exp{(-2\pi{ijf})}}|^2\\ \sum_{j=0}^{\infin}{g_j\exp{(-2\pi{ijf})}} &= \bigg({1 - \sum_{j=1}^m{a_j\exp{(-2\pi{ijf})}}}\bigg)^{-1} \bigg({1 - \sum_{j=1}^l{b_j\exp{(-2\pi{ijf})}}}\bigg)\\ p(f) &= \sigma^2\frac{|{1 - \sum_{j=1}^l{b_j\exp{(-2\pi{ijf})}}}|^2}{|{1 - \sum_{j=1}^m{a_j\exp{(-2\pi{ijf})}}}|^2} \end{align}
  • 結局、スペクトルが表現したいのは、波であるが、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個の山が必要なときは、ARの次数は2k以上必要となる。谷の場合も同様で、MAの次数は2k以上必要となる。


特性方程式


ARオペレータa(B)とMAオペレータ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がARモデルの式で表現できるということ。
        \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}
  • ここで、今後、便利な方程式として利用したいため、g(B)Bの部分(Back-shift Operator)を、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_nl変量正規ホワイトノイズ
      • E[v_n] = \begin{bmatrix}0 \\ \vdots \\ 0s \end{bmatrix}:平均は0(正規分布)
      • 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≠mv_nは自分自身と独立(無相関)
      • E[v_ny_m^T] = 0 ~~~~~ n>mv_ny_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が多変量ARモデルに従うときは、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に近づくというもの。
coh_{jk}(f) = \frac{a_{jk}(f)^2}{p_{jj}(f)p_{kk}(f)} ~~~~~ コヒーレンシー
名称 数式 中身
振幅スペクトル a_{jk}(f) =\sqrt{c_{jk}(f)^2 + d_{jk}(f)^2}
位相スペクトル \phi_{jk}(f) =arctan(\frac{d_{jk}(f)}{c_{jk}(f)})
クロススペクトル p_{jk}(f) = a_{jk}(f)\exp(i\phi_{jk}(f))

フィードバックシステム・パワー寄与率


フィードバックシステム

  • 多変量の時系列データの場合、ある変数Aに対して、ある変数Bが影響を及ぼすことがある。逆も然りで、変数Bに対して、変数Aが影響を及ぼすこともある。

  • これをフィードバックシステムという。

    • 例1:農業生産物の生産量が上がると価格が下がるが、その後生産量下げることで価格が上がる。
    • 例2:石油価格が上がるとエコカーの販売量が上がり、石油価格が下がるとエコカーの販売量が下がる。
  • これは、つまり、インパルス応答関数の形が変わることを意味する。

    • 1変量の場合:開ループ。ある周波数に対して、変数Aを通した時のその後の周波数の変化を見るが、それは一方通行である。
    • 多変量の場合:閉ループ。変数Aを通すまでは同様だが、その出力を変数Bへ渡す。さらに、変数Bの出力を変数Aへ返すという風に、ループが発生する。

2変量の場合のフィードバックシステムを考えてみる

  • まず、y_n,x_nをそれぞれフィードバックシステムに組み込んだ時の数式として表現する。
    • 外乱:外部からの干渉というようなイメージ。具体的には、変数Aを通してでた出力から変数Bの入力までの間に、外乱u_nが加わるということ。ホワイトノイズとして、処理したいが、このままだとx_nとの間に相関がないとは言えないため、その処理は難しい。
\begin{align}y_n &= \sum_{j=1}^{\infin}{a_jx_{n-j}+u_n}\\ x_n &= \sum_{j=1}^{\infin}{b_jy_{n-j}+v_n}\end{align}
名称 数式
ある変数Aのインパルス応答関数 a_j
ある変数Bのインパルス応答関数 b_j
互いに無相関な外乱 u_n,v_n
  • x_n,y_nの観測値を用いて、インパルス応答関数a_j,b_jと外乱u_n,v_nの特性を推定していく。

  • これを解いていくのであれば、一見、通常の最小二乗法でいい気もする。ただ、そのためには、誤差項の部分(今回でいえば外乱)が独立である必要があった。

  • 現状、x_nu_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について代入をし、式を1つにする。

\begin{align}y_n - \sum_{j=1}^{\infin}{a_jx_{n-j}} &= \sum_{i=1}^{\infin}{c_i(y_{n-i} - \sum_{j=1}^{\infin}{a_jx_{n-i-j}}) + \varepsilon_n}\\ y_n &= \sum_{i=1}^{\infin}{c_iy_{n-i}} + \sum_{j=1}^{\infin}{a_jx_{n-j}} - \sum_{i=1}^{\infin}\sum_{j=1}^{\infin}{c_ia_jx_{n-i-j}} + \varepsilon_n \end{align}
  • さらに、この式を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

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