時系列解析　ことはじめ ③

概要

本ポストでは, 定常時系列モデルのパラメータの推定方法を紹介する. 最近のRやPythonの時系列解析パッケージは, 難しい推定方法を気にしなくてよいほど高度な機能を手軽に試せる.
が, パラメータ推定法にはいろいろと学びや応用がきく内容があるので, あえてゼロから考えていくこととする.

一変量ARモデルの推定方法

ここではAR(m)モデルを例に考えていく. ARMAモデルを扱わない理由は, ARMAモデルの推定には状態空間モデルがとても便利であるため, 状態空間モデル（後日紹介）を説明した際にARMAモデルの推定を行うとしよう.また扱うモデルは1変量時系列モデルに限定している.

改めて,AR(m)モデルの式を提示する.AR(m)モデルでは, $n$時点の時系列データ$y_n$を, $m$時点前からのデータ$y_{n-m}, y_{n-m+1},\dots,y_{n-1}$と, 白色雑音$\nu_n$の線形和で表した以下の式を利用する.

$$y_n = \sum_{j=1}^{m} a_j y_{n-j} + \nu_n, \nu_n \sim N(0, \sigma^2)$$

このAR(m)モデルでのパラメータ推定対象は, m個のAR係数$a_j$と分散$\sigma^2$である($m+1$個), また, 過去どの時点までのデータを考慮するのが最良か, という次数$m$の選択もある.

ARモデルでは複数の推定方法があり, 本ポストでは, 『時系列解析』の講座で取り上げられた4つの方法, 1)Yule-Walker法, 2)最尤法, 3)最小二乗法, 4)PARCOR法を解説する.

Yule-Walker法

Yule-Walker法では, 以下の自己共分散関数$C_k$と,AR係数$a_j$の一次方程式を解くことで推定値を得るというものである.

$$\begin{align} C_0 &= \sum_{j=1}^{m}a_j C_j + \sigma^2 \\ C_k &= \sum_{j=1}^{m} a_j C_{j-k}, (k = 1,2,...,m-1) \\ \end{align}$$

趣味の人向けの導出

AR(m)モデルの両辺に, $n-k$時点の時系列データ$y_{n-k}$をかけて期待値をとると,

$$\begin{align} E(y_n y_{n-k}) &= \sum_{j=1}^{m} a_j E(y_{n-j}y_{n-k}) + E(\nu_{n} y_{n-k}) \notag\\ \Leftrightarrow C_{k} &= \sum_{j=1}^{m} a_j C_{j-k} + E(\nu_{n} y_{n-k}) \notag\\ \end{align}$$

となる. ここでARモデルの定義を思い出すと, $n$時点の時系列データ$y_n$と, それよりも以後の白色雑音$\nu_{n+j} (j = 1,2,...)$は無相間($E(\nu_{n+j}y_n ) = 0$)であることや, どう時点の場合は$E(\nu_{n}y_n ) = \sigma^2$であることを使うと,

$$\begin{align} C_0 &= \sum_{j=1}^{m} a_j C_{j} + \sigma^2 \\ C_{k} &= \sum_{j=1}^{m} a_j C_{j-k} (k = 1,2,...) \\ \end{align}$$

となることにより導出された.

特に式(2)に着目すると.

$$\begin{align} C_1 &= a_1 C_0 + a_2 C_1 + \dots + a_m C_{m-1} \notag\\ C_2 &= a_1 C_{-(1)} + a_2 C_0 + \dots + a_m C_{m-2} \notag\\ \vdots &= \vdots \notag\\ C_{m-1} &= a_1 C_{-(m-2)} + a_2 C_{-(m-1)} + \dots + a_m C_{0} \notag\\ \end{align}$$

一変量時系列では, 自己共分散関数は偶関数($C_k = C_{-k}$)より,以下の一次方程式に書き直せる.

$$\begin{align} C_1 &= a_1 C_0 + a_2 C_1 + \dots + a_m C_{m-1} \notag\\ C_2 &= a_1 C_{1} + a_2 C_0 + \dots + a_m C_{m-2} \notag\\ \vdots &= \vdots \notag\\ C_{m-1} &= a_1 C_{(m-2)} + a_2 C_{(m-1)} + \dots + a_m C_{0} \notag\\ \end{align}$$

いま, 自己共分散関数$C_k$を観測データからの推定値$\hat{C}_k$に置き換え, かつ上式を行列表現に変えると,

$$\begin{pmatrix} \hat{C}_1 \\ \hat{C}_2 \\ \vdots \\ \hat{C}_m \end{pmatrix} = \begin{pmatrix} \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{pmatrix} \begin{pmatrix} a_1 \\ a_2 \\ \vdots \\ a_m \end{pmatrix}$$

と書き直せる. つまり, 自己共分散関数の推定値$\hat{C}_k$が手元にあれば, 上式のm次の一次方程式を解くことでAR係数の推定量が得られる. 分散$\sigma^2$については, 新たに推定されたAR係数と自己共分散関数の推定値が出そろえば, 最後に求めることができる(式(1)に代入するだけ）

さてAR係数$a_j$を推定するための式が出てきたので, あとはm次の一次方程式を解くだけだが, ここで計算量について考えると, 最良のARモデルを探すために次数$m$も探索する場合には, 計算量は$O(m^4)$となってくる. 次数の高いARモデルになると計算量がボトルネックとなる.

そこで, 計算量を$O(m^2)$程度まで抑えるLevinson's algorithmを使えると便利

Levinson's algorithmによる効率的な解法

さきほどの通常のYule-Walker法では, 次数$m$ごとに自己共分散関数$C_k$とAR係数$a_j$の$m$元一次方程式を毎回解いている.

Levinson's algorithmでは, 毎回$m$元一次方程式を最初から解くようなことはせず, 次数$m-1$のときのAR係数$a_{j}^{m-1}$を活用し, 次数$m$でのARモデルのAR係数$a_{j}^{m}$を逐次的に求めるという方法をとる.

具体的には, $m = 0$からスタートして, $m = M$までパラメータの推定量を逐次的に求めていく.

1. $m = 0$のときの初期値を用意

$$\begin{align} \sigma^2 &= C_0 \\ AIC_0 &= N(\log{2\pi \hat{\sigma}_{0}^2 + 1}) + 2 \\ \end{align}$$

2.$m = 1,2,\dots,M$のとき以下の反復式を解く

$$\begin{align} a_{m}^{m} &= (C_m - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{m-j})(\sigma_{m-1}^2)^{-1} \\ a_{j}^{m} &= a_{j}^{m-1} - a_{m}^{m} a_{m-j}^{m-1} \\ \sigma_{m}^2 &= \sigma_{m-1}^2 (1 - (a_{m}^{m})^2) \\ AIC_m &= N(\log{2\pi \hat{\sigma}_{m}^2 + 1}) + 2(m+1) \\ \end{align}$$

m=2ぐらい求めてみる

Levinson's algorithmのイメージを掴むために, $m=0,1,2$ぐらいまで列挙してみる.

$m=0$のとき

$$\begin{align} \sigma^2 &= C_0 \notag\\ AIC_0 &= N(\log{2\pi \hat{\sigma}_{0}^2 + 1}) + 2 \notag\\ \end{align}$$

$m=1$のとき

$$\begin{align} a_{1}^{1} &= (C_1 - 0)(\hat{\sigma}_{0}^{2})^{-1} \notag\\ \sigma_{1}^2 &= \sigma_{0}^2 (1-(a_{1}^{1})^2) \notag\\ AIC_1 &= N(\log{2\pi \hat{\sigma}_{1}^2 + 1}) + 2(1+1) \notag\\ \end{align}$$

$m=2$のとき

$$\begin{align} a_{2}^{2} &= (C_2 - a_{1}^{1}C_1)(\hat{\sigma}_{1}^{2})^{-1} \notag\\ a_{1}^{2} &= a_{1}^{1} - a_{2}^{2} a_{1}^{1} \notag\\ \sigma_{2}^2 &= \sigma_{1}^2 (1-(a_{2}^{2})^2) \notag\\ AIC_2 &= N(\log{2\pi \hat{\sigma}_{2}^2 + 1}) + 2(2+1) \notag\\ \end{align}$$

...
のように, 上から一つずつ計算していけば良い.

最尤法

さて次は最尤法でのパラメータ推定方法を説明する. 最尤法では,パラメータ$\theta = (a_1,a_2,\dots,a_m,\sigma^2)$を直接推定する. そして, 時系列モデルの残差部分が正規分布$N(0, \sigma^2)$に従うことを仮定している. つまり,

$$p(\nu_n) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp \Big({-\frac{1}{2\sigma^2}\nu_{n}^{2}}\Big)$$

とかける. ここで, 残差$\nu_n$自体は, 時系列データ$y_n$およびAR係数$a_j$から求まるため,上の確率モデルを, 時系列データ$y_n$およびAR係数$a_j$を使って書き直した以下の式を使う.

$$\begin{align} p(y_n | \theta, y_1, y_2, \dots, y_N) &= \frac{1}{ \sqrt{2\pi \sigma^2}} \exp \Big\{ - \frac{1}{2\sigma^2} \Big( y_n - \sum_{j=1}^{m} a_j y_{n-j} \Big)^{2} \Big\} \\ \log p(y_n | \theta, y_1, y_2, \dots, y_N) &= -\frac{1}{2} \log (2\pi \sigma^2) - \frac{1}{2\sigma^2} \Big( y_n - \sum_{j=1}^{m} a_j y_{n-j} \Big)^{2} \\ \end{align}$$

次に, ARモデルから,時系列データ$y_n$に対する尤度および対数尤度を求める.

$$\begin{align} L(\theta) &= p(y_1, y_2, \dots, y_N | \theta) \\ &= p(y_1 | \theta)p(y_2 | \theta, y_1) \cdots p(y_N | \theta, y_1, y_2, \dots, y_{N-1}) \\ &= \prod_{n=1}^{N} p(y_n | \theta, y_1, y_2,\dots, y_{n-1}) \\ l(\theta) = \log L(\theta) &= \sum_{n=1}^{N} \log p(y_n | \theta, y_1, y_2,\dots, y_{n-1}) \\ &= \sum_{n=1}^{N} \Big\{ -\frac{1}{2} \log (2\pi \sigma^2) - \frac{1}{2 \sigma^2} \Big( y_n - \sum_{j=1}^{m} a_j y_{n-j} \Big)^{2} \Big\} \\ &= -\frac{N}{2}\log(2\pi \sigma^2) - \frac{1}{2\sigma^2}\sum_{n=1}^{N}\Big( y_n - \sum_{j=1}^{m} a_j y_{n-j}\Big)^2 \end{align}$$

この対数尤度をパラメータ$\theta = (a_1,a_2,\dots,a_m,\sigma^2)$に対し最大化することを考える.

例えば, 分散$\hat{\sigma}^2$を求めるには, 対数尤度を分散$\sigma^2$に関して微分し,

$$\begin{align} \frac{\partial l(\theta)}{\partial \sigma^2} &= -\frac{N}{2\sigma^2} + \frac{1}{2(\sigma^2)^2} \sum_{n=1}^{N} \Big( y_n - \sum_{i=1}^{m} a_j y_{n-j} \Big)^2 = 0 \\ \rightarrow \hat{\sigma}^2 &= \frac{1}{N}\sum_{n=1}^{N} \Big( y_n - \sum_{i=1}^{m} a_j y_{n-j} \Big)^2 \end{align}$$

と求まる. 同様に, AR係数$a_j$で対数尤度を微分した式を解くことで, AR係数$a_j$の推定量が得られる. 実際には解析的に求まらないような複雑な尤度関数となることも多く, その場合は数理最適化手法を用いることとなる.

しかし, もう一点注意が必要ななのは, $n \leq m$の場合である.
いくつか対処法があり,

1. $y_1, \dots, y_m$だけは, 同時分布$p(y_1, y_2, \dots, y_m | \theta)$で正直に解くという方法
2. より低次のAR係数を使って, 同じく条件付き分布を解くという方法がある

『時系列解析』の講座で北川先生がさらっとおっしゃっていたので完全に理解できているわけではないが, 次の理由からだろう.

$m$次のARモデルのAR係数$a_{1}^{m},a_{2}^{m},\dots, a_{m}^{m}$がもとまっていたら,$m-1$次のARモデルのAR係数$a_{1}^{m-1},a_{2}^{m-1},\dots, a_{m-1}^{m-1}$が逆算できる.これを使うと,

$$y_m = \sum_{j=1}^{m-1}a_{j}^{m-1} y_{m-j}$$

という推定式が作れるので,

$$p(y_m | \theta, y_1, y_2, \dots, y_{m-1}) = -\frac{1}{2} \log (2\pi \sigma^2) - \frac{1}{2\sigma^2} \Big( y_m - \sum_{j=1}^{m-1}a_{j}^{m-1} y_{m-j} \Big)^{2}$$

と, $n \> m$のときと同じように条件付き分布を使って計算できる.
同じく, $p(y_{m-1} | \theta, y_1, y_2, \dots, y_{m-2})$を出すなら, $m-1$次のARモデルのAR係数から, $m-2$次のARモデルでのAR係数を求める...ように再帰的に求めたら解決する.

しかし, とても面倒くさいので厳密な最尤法かつ楽な方法として状態空間モデルが便利, だそう.

最小二乗法

真正面から最小二乗法を考える前に,先ほどの尤度法の式を一旦近似することを考える.
さきほどは, $n \leq m$の関係のとき, 条件付き分布がややこしいと述べた.

なので, $n \leq m$の場合をまったく無視してみると,

$$\begin{align} p(y_n | y_1, y_2, \dots, y_{n-1}, \theta) &= p(y_n | y_{n-m}, \dots, y_{n-1}, \theta) \notag\\ &= \frac{1}{ \sqrt{2\pi \sigma^2}} \exp \Big\{ - \frac{1}{2\sigma^2} \Big( y_n - \sum_{j=1}^{m} a_j y_{n-j} \Big)^{2} \Big\} \notag\\ l(\theta) &= \sum_{n=1}^{N}\log p(y_n | y_1, y_2, \dots, y_{n-1}, \theta) \notag\\ &\simeq \sum_{n=m+1}^{N}\log p(y_n | y_1, y_2, \dots, y_{n-1}, \theta) \notag\\ &= -\frac{N-m}{2}\log(2\pi \sigma^2) - \frac{1}{2\sigma^2}\sum_{n=m+1}^{N}\Big( y_n - \sum_{j=1}^{m} a_j y_{n-j}\Big)^2 \notag \end{align}$$

と書き直せる.この近似尤度をパラメータ$a_j, \sigma^2$で微分して両辺を0とすると

$$\begin{align} \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} \Big( y_n - \sum_{i=1}^{m} a_j y_{n-j} \Big)^2 = 0 \\ \rightarrow \hat{\sigma}^2 &= \frac{1}{N-m}\sum_{n=m+}^{N} \Big( y_n - \sum_{i=1}^{m} a_j y_{n-j} \Big)^2 \end{align}$$

推定した$\hat{\sigma}^2$を尤度関数に代入すると,

$$l(\sigma^2, a_j) = -\frac{N-m}{2}\log (2\pi \hat{\sigma}^2) - \frac{N-m}{2}$$

となる. 次に, この尤度関数をAR係数$a_j$で微分して, 両辺を0とおいた式から推定量$\hat{a_j}$が得られる. しかし, よく式をみると, この尤度関数をAR係数$a_j$で最小化するといいうことは, 推定量$\hat{\sigma}^2$をAR係数$a_j$で最小化することと一致する. さらに, 推定量$\hat{\sigma}^2$の式に戻ると,

$$\max_{a_j}l'(\hat{\sigma}^2,a_j) \Leftrightarrow \min_{a_j}\hat{\sigma}^2 \Leftrightarrow \min_{a_j}\sum_{n=m+1}^{N}(y_n - \sum_{j=1}^{m} a_j y_{n-j})^2$$

一番右の項は, ARモデルの誤差項となっていることから, 誤差項の二乗の最小化を求めていることになる. どうしたことか, 最尤法を少し緩和すると, みなさんご存知最小二乗法が出てきた.

さて, あとはいつものように行列式をもちだして

$$y = \begin{pmatrix} y_{m+1} \\ \vdots \\ y_N \end{pmatrix}, Z = \begin{pmatrix} y_{m} & \cdots & y_1 \\ \vdots & & \vdots \\ y_{N-1} & \cdots & y_{N-m} \\ \end{pmatrix}, a = \begin{pmatrix} a_{1} \\ \vdots \\ a_m \end{pmatrix}, \nu = \begin{pmatrix} \nu_{m+1} \\ \vdots \\ \nu_N \end{pmatrix},$$

$y = Za + \nu$と書き直せ, $\parallel \nu \parallel^2 = \parallel y - Za \parallel^2$ を最小化することで, $a = (Z^{\top}Z)^{-1}Z^{\top}y$と求まる.

Householder法

この最小二乗法を効率よく方法としてHouseholder法を説明する.

いま任意の直交変換$U$を考える.$U$は直交変換なので, $y$および$Za$に変換を加えた$Uy$および$UZa$の内積は保存される. つまり, $\parallel \nu \parallel^2 = \parallel y - Za \parallel^2 = \parallel Uy - UZa \parallel^{2}$

$X$を$[z | y]$と行列（ベクトル）を埋め込んだ$N \times m + 1$次元の行列とおき, 任意の直行変換$U$がこの$X$を上三角行列に変換するようなものを探す. つまり,

$$UX = S = \begin{pmatrix} s_{11} & \cdots & s_{1,m} & s_{1, m+1} \\ & \ddots & \cdots & \cdots \\ & & s_{m,m} & s_{m, m+1} \\ & & & s_{m+1, m+1} \\ & 0 & & \\ \end{pmatrix}$$

となる変換$U$が見つかったとする.

すると,

$$\begin{align} \parallel Uy - UZa \parallel_{N}^{2} &= \begin{pmatrix} s_{1, m+1} \\ \vdots \\ s_{m, m+1} \\ s_{m+1, m+1} \\ 0 \end{pmatrix} - \begin{pmatrix} s_{1, 1} & \cdots & s_{1,m} \\ & \ddots & \vdots \\ & & s_{m, m} \\ & 0 & 0 \\ \end{pmatrix} \begin{pmatrix} a_1 \\ \vdots \\ a_m \\ \end{pmatrix} \parallel_{N}^{2} \\ &= \parallel \begin{pmatrix} s_{1, m+1} \\ \vdots \\ s_{m, m+1} \\ \end{pmatrix} - \begin{pmatrix} s_{1, 1} & \cdots & s_{1,m} \\ & \ddots & \vdots \\ 0 & & s_{m, m} \\ \end{pmatrix} \begin{pmatrix} a_1 \\ \vdots \\ a_m \end{pmatrix} \parallel_{m}^{2} + s_{m+1, m+1}^2 \end{align}$$

後ろの$s_{m+1, m+1}^2$は十分に小さいので無視すると, 結局次の式

$$\begin{pmatrix} s_{1, 1} & \cdots & s_{1,m} \\ & \ddots & \vdots \\ 0 & & s_{m, m} \\ \end{pmatrix} \begin{pmatrix} a_1 \\ \vdots \\ a_m \end{pmatrix} = \parallel \begin{pmatrix} s_{1, m+1} \\ \vdots \\ s_{m, m+1} \\ \end{pmatrix}$$

を解くことで, 推定量$\hat{a_j}$が得られる. そして, 左辺に上三角行列があることから, 推定量$\hat{a_j}$は次数$m$から後ろ向きに解くことができる.

$$\begin{align} \hat{a}_m &= \frac{s_{m,m+1}}{s_{m,m}} \\ \hat{a}_{i} &= \frac{s_{i, m+1} - s_{i, i+1}\hat{a}_{i+1} - \dots s_{i,m}\hat{a}_{m} }{s_{i,i}}, i = m-1, ..., 1 \\ \hat{\sigma^2} &= \frac{s_{m+1,m+1}^2}{n} \end{align}$$

求めた$\hat{a}_{i}$をつかって, $\hat{a}_{i-1}$を解くというような反復処理と推定量がすべて求まるになっている. ベーシックな逆行列を使った解法と異なり, 逆行列を求めないことから計算コストが良いのが魅力だ.

実際にHouserholder法を使えば, 簡単な手続きで推定量が得られることがわかった.ただ問題としては, この$X = [Z | y]$を上三角行列にする直交変換$U$をどう作るかだ

実はそのような直交変換$U$の作り方は決まっている. 例えば, ある単位ベクトル$w$があるとき, $w$に直交する平面$M$を考える. 次に, ベクトル$a$と, ベクトル$b$を用意する. ここで, ベクトル$a$の平面$M$の鏡映変換$U$が, $Ua = b$という関係になっているとする. つまり, ベクトル$a$とベクトル$b$が平面$M$を境に鏡写しになっている.

そのような$U$を $U = I - 2ww^{\top}$とおくと, $U$は直交変換になっている.
実際

$$\begin{align} UU^{\top} &= (I - 2ww^{\top})(I - 2ww^{\top})^{\top} \\ &= I - 4ww^{\top} + 4ww^{\top}ww^{\top} \\ &= I - 4ww^{\top} + 4ww^{\top} = I\\ \end{align}$$

となっている. 式の移り変わりで, $w$が単位ベクトルであることから, $w^{\top}w = I$となることを使った.

たしかに, 単位ベクトル$w$を使って, $U = I - 2ww^{\top}$と書けば, $U$は直交変換となっている. では, 次の問題として,$Ua = b$となるような直交変換$U$を存在するだろうか, 存在するとすればこの直交変換$U$をどう作るか, もっと具体的には, どのような単位ベクトル$w$　を定義したらよいだろうか？

ここで, ベクトル$a$とベクトル$b$が同じ長さであれば, $Ua = b$となるような直交変換$U$を存在する. 実際, $w = \frac{a-b}{\parallel a - b \parallel}$とおけば,

$$\begin{align} Ua &= (I - 2ww^{\top})a \\ &= a - \frac{(a-b)(a-b)^{\top} }{\parallel a - b \parallel^2} 2a \\ &= a - \frac{(a-b)(a-b)^{\top} }{\parallel a - b \parallel^2} {(a-b) + (a +b)} \\ &= a - (a-b) - \frac{(a-b)(\parallel a \parallel^2 \parallel b \parallel^2)^{\top} }{\parallel a - b \parallel^2} \\ &= b \end{align}$$

と確かに変換されている.

再び『時系列解析』の講座では, このHouseholder変換をつかって, 行列を三角行列にする方法を説明されていた. しかし, 北川先生が解説だけでは私自身まだ理解が及ばず, 私なりに理解した三角化法を説明する.

まず, さきほどのHouseholder変換の良さは, $Ua = b$となるような直交変換$U$は, ベクトル$a$とベクトル$b$が同じ長さであれば存在するという主張がされていることだ. つまり, ベクトル$a$とベクトル$b$の長さが同じであれば, 各ベクトル成分の配分に特段指定はないのがミソだと理解した.^[1]

このことを理解し, $m \times m$の行列$X$を三角化することを考える. いま, 行列$X$の一番目の列ベクトルに着目し, 残りの要素を$X_{-a}$とおく.つまり,

$$X = \begin{pmatrix} x_{1,1} & x_{1,2} & \cdots & x_{1,m} \\ x_{2,1} & x_{2,2} & \cdots & x_{2,m} \\ \vdots & \vdots & & \vdots \\ x_{n,1} & x_{n,2} & \cdots & x_{n,m} \\ \end{pmatrix} = [a | X_{-a}]$$

で, 列ベクトル$a_1$ は, $(x_{1,1}, x_{2,1},...,x_{n,1})^{\top}$となる.
ここで, $a_1$と同じ大きさの列ベクトル$b_1$を$(x'_{1,1}, 0,...,0)^{\top}$とおく.ベクトル$a_1$とベクトル$b_1$は同じ大きさなので, $x'_{1,1} = \mp \parallel a_1 \parallel$と置き換えられる. そして, ベクトル$a_1$からベクトル$b_1$に変換する直交変換$U_1$は, $w_1 = \frac{a_1-b_1}{\parallel a_1 - b_1\parallel}$と定めて, $U = I - 2w_1w_{1}^{\top}$で求まる.

では, この$U$を使って$X$を変換すると,

$$U_1X = [U_1a_1 | U_1X_{-a} ] = \begin{pmatrix} x'_{1,1} & x'_{1,2} & \cdots & x'_{1,m} \\ 0 & x'_{2,2} & \cdots & x'_{2,m} \\ \vdots & \vdots & & \vdots \\ 0 & x'_{n,2} & \cdots & x'_{n,m} \\ \end{pmatrix}$$

となり, 一番左の列が, $(1,1)$成分を除いてすべて 0に置き換わった.^[2]

同じく,次は２列目のベクトルを$a_2 =(x'_{1,2}, x'_{2,2},...,x'_{n,2})^{\top}$とおき, 変換後のベクトルを$b_2 =(x'_{1,2}, x''_{2,2},0, ...,0)^{\top}$とするような直交行列$U_2$を作る.そして, 今度は$U_1X$に左からかけると

$$U_2U_1X = [U_2b_1 | U_2a_2 | U_2X_{-b}] = \begin{pmatrix} x'_{1,1} & x'_{1,2} & \cdots & x''_{1,m} \\ 0 & x''_{2,2} & \cdots & x''_{2,m} \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & x'_{n,m} \\ \end{pmatrix}$$

とさらに三角化行列に近づいた...

という変換を$m$列ベクトルまで続けると, 最終的には$U = U_mU_{m-1}...U_{1}$が, 元の行列$X$を上三角化行列にする変換行列として得られる.

気になる人向け

説明を省いたが, $U_2b_1 = b_1$となっていることを確認してみよう.
まず, $b_1$は, 第一成分だけがある列ベクトルである($(x'_{1,1}, 0,...,0)^{\top}$)

次, ベクトル$a_2$からベクトル$b_2$へ変換する$U_2$をみてみる.
$w_2 = \frac{a_2 - b_2}{\parallel a_2 - b_2 \parallel}$であるが, ここで, ベクトル$a_2$からベクトル$b_2$の第一成分は同じ値なので, 単位ベクトル$w_2$の第一成分は0となる.ここでは,$w_2 = (0, w_{2},w_{3},...,w_{m})^{\top})$とおいておく.

つぎに, この$w_2$を使って変換行列$U_2$を求めたい. ここで, $w_2 w_{2}^{\top}$は次のような形となっている.

$$w_2 w_{2}^{\top} = \begin{pmatrix} 0 & 0 & \cdots & 0 \\ 0 & w_{2,2} & \cdots & w_{2,m} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & w_{m,2} & \cdots & w_{m,m} \\ \end{pmatrix}$$

さて, $U_2 = I - 2w_2 w_{2}^{\top}$をもとめると,

$$U_2 = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 - 2w_{2,2} & \cdots & -2w_{2,m} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & -2w_{m,2} & \cdots & 1- 2w_{m,m} \\ \end{pmatrix}$$

$U_2 b_2$の計算では, $U_2$の一行目(=$(1,0,...,0)$)との掛け算のみ値が反映され, 他の$U_2$の行とは積が0になるので, 結局, ベクトル$b_1$の第一成分しか残らない.
よって, $U_2b_1 = b_1$となるのでした. おわり

PARCOR法

最後にPARCOR法について説明する.途中まではYule-Walker法と同じで,

$$\begin{align} a_{j}^{m} &= a_{j}^{m-1} - a_{m}^{m} a_{m-j}^{m-1} \notag\\ a_{m}^{m} &= (C_m - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{m-j})(\sigma_{m-1}^2)^{-1} \notag\\ \end{align}$$

を解くことを考える. Yule-Walker法では, Levinson's algorithmを使い, 時系列データ$y_n$から事前に推定された自己共分散$\hat{C}_{k}$を用いて推定量$\hat{a}_k$を求めていた.

PARCOR法では, 自己共分散$\hat{C}_{k}$を用いず, データから直接推定することを考える.^[3]

$$\begin{align} C_m - \sum_{j=1}^{m-1}a_{j}^{m-1}C_{m-j} &= E\{(y_n - \sum_{j=1}^{m-1}a_j^{m-1}y_{n-j})y_{n-m} \} \notag\\ &= E(\nu_{n}^{m-1}y_{n-m})\notag\\ &= E(\nu_n^{m-1}w_{n-m}^{m-1})\notag\\ &\simeq \frac{1}{N-m} \sum_{n=m+1}^{N}(\nu_n^{m-1}w_{n-m}^{m-1}) \end{align}$$

と, 推定量$a_{m}^{m}$を得るための一部の式を**自己共分散$\hat{C}_{k}$を用いず出せた.
続けて, 推定量$\hat{\sigma}^2$については,

$$\begin{align} \hat{\sigma}_{m-1}^2 = C_0 - \sum_{j=1}^{m-1} a_j^{m-1}C_j &= E\{(y_{n-m}- \sum_{j=1}^{m-1}a_j^{m-1}y_{n-m + j})y_{n-m} \} \notag\\ &= E(w_{n-m}^{m-1}y_{n-m}) \notag\\ &= E(w_{n-m}^{m-1})^2 \notag\\ &\begin{equation} \simeq \begin{cases} \frac{1}{N-m} \sum_{m+1}^{N} (w_{n-m}^{m-1})^2 \notag\\ \frac{1}{N-m} \{ \sum_{m+1}^{N} (w_{n-m}^{m-1})^2 \sum_{m+1}^{N} (\nu_{n-m}^{m-1})^2 \}^{1/2} \notag\\ \frac{1}{2(N-m)} \{ \sum_{m+1}^{N} (w_{n-m}^{m-1})^2 + \sum_{m+1}^{N} (\nu_{n-m}^{m-1})^2\} \notag\\ \end{cases} \end{equation} \end{align}$$

と表すことができる. 最後の式で３パターンにわかれたのは, $E(w_{n-m}^{m-1})^2$の近似方法がいろいろあるためである. 一番上は, Partial Regressionとよばれ, 標本平均を求めている.真ん中の近似では, 相乗平均, 一番下では, 相加平均となっている. それぞれ, PARCOR法, Burg法(MEM) と呼ばれている. これらのパターンを考慮して, さきほど求めた推定量$a_m^m$の式に代入すると,

$$\begin{equation} a_m^m = \begin{cases} \sum_{n=m+1}^{N} (\nu_n^{m-1} w_{n-m}^{m-1}) \Big\{ \sum_{m+1}^{N} (w_{n-m}^{m-1})^2 \Big\}^{-1} \\ \sum_{n=m+1}^{N} (\nu_n^{m-1} w_{n-m}^{m-1}) \Big\{ \sum_{m+1}^{N} (w_{n-m}^{m-1})^2 \sum_{m+1}^{N} (\nu_{n-m}^{m-1})^2 \Big\}^{-1/2} \\ \sum_{n=m+1}^{N} (\nu_n^{m-1} w_{n-m}^{m-1}) \Big\{ \sum_{m+1}^{N} (w_{n-m}^{m-1})^2 + \sum_{m+1}^{N} (\nu_{n-m}^{m-1})^2 \Big\}^{-1} \\ \end{cases} \end{equation}$$

と,$a_m^m$を**自己共分散$\hat{C}_{k}$を用いず推定する式ができた. あとは, データから直接$a_m^mを求めることができる.

まとめ

本ポストでは, 定常時系列モデルの推定方法を紹介した. 特に一変量ARモデルでは, 1)Yule-Walker法, 2)最尤法, 3)最小二乗法, 4)PARCOR法があった.

ちなみに, 『時系列解析』の講座では, 各種法の特徴を次の軸で示していた.

	定常性	数値精度	モデル精度	速度	自由度
Yule-Walker法	○	△	△	○	△
最尤法	◎	○	◎	×	○
最小二乗法	×	◎	○	◎	◎
PARCOR法	◎	○	○	◎	△

最小二乗法についている「定常性×」は, 最小二乗法で推定されたARモデルは定常ではないかもしれないということらしい（？）. ただ, 非定常な時系列データには非定常なモデルとして推定するため, 良い面もある. その他に, Householder法での解法がとても便利なので, 注目しても良いかと思う.

次回は, ようやくだが, 時系列モデルやプログラム実装を試していこう.

雑記

数式ばかり書いていたらとんでもない文字量になったので, 一度ポストを切り分けた💦
北川先生の解説はわかりやすいが, たまに理解が追いつけなくなる部分もある（それは私自身の知識不足が原因だが,,,). 改めて解説を自己解釈したり, 数式を書いてみることも大事だとおもった.

参考

https://ocw.u-tokyo.ac.jp/course_11477/

いつもお世話になっている時系列講座. 本ポストはだいたい4章,6章辺りの内容を解説している.

https://ocw.kyoto-u.ac.jp/wp-content/uploads/2021/03/2004_gazoushoriron_06.pdf

直交変換の勉強のために参照した

脚注

1. もちろん, $a = b$となるのは無し

2. $U_1X_{-a}$については, どのような変換結果になっていようが（現段階では）関心がない. まずは,左からひとつずつ三角行列に近づけていっている..と考えよう.

3. Levinson's algorithmでは推定量の自己共分散$\hat{C}_{k}$を使って,$\hat{a}_k$を求めていたので, 推定精度に問題があるかもしれない,...というリスクを回避しようと試みといえる