Time Series Analysis：ノート5章最尤推定（その1）

以下の本の5章を読みます。キャラクターは架空のものです。私の誤りは私に帰属します。お気付きの点がありましたらコメントでご指摘いただけますと幸いです。

作者:Hamilton, James D.
発売日: 1994/01/11
メディア: ハードカバー

f:id:cookie-box:20190101155733p:plain:w60

5章はそのタイトルも「最尤推定」ですね。これまでは ARMA 過程

$Y_t = c + \phi_1 Y_{t-1} + \cdots + \phi_p Y_{t-p} + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \cdots + \theta_q \varepsilon_{t-q}$ のパラメータ $\theta = ( c, \, \phi_1, \, \cdots, \, \phi_p, \, \theta_1, \, \cdots, \, \theta_q, \, \sigma^2 )^\top$ が所与の下で、 $Y_t$ がどんな性質をもつかとか、あるステップまで観測した下で先のステップの値がどのように予測できるかを扱ったのですよね。対してこの章ではパラメータを未知として、 $Y_t$ の観測値 $(y_1, \cdots, y_T)$ からそれを推定しにいくということです。

これまで： $\theta$ が所与の下で $Y_t$ の期待値や自己共分散、条件付き期待値を調べた。
ここでは： $\theta$ が未知の下で $Y_t$ の観測値 $(y_1, \cdots, y_T)$ から $\theta$ を推定する。

そして $\theta$ の推定のやり方としては、 $(Y_1, \cdots, Y_T)$ の確率密度関数 $f_{Y_T, \, \cdots, \, Y_1}(y_T, \, \cdots, \, y_1 \, ; \, \theta)$ の観測値における値を最大にするような $\hat{\theta}$ を推定値とするということです。なんと、そのような考え方があったとは…。

f:id:cookie-box:20190101160814p:plain:w60

それが「最尤推定」だからね。その確率密度関数は、パラメータ $\theta$ の関数と考えると $\theta$ の尤度関数というね。尤度関数の最大点 $\hat{\theta}$ は $\theta$ の最尤推定量だね。最尤推定量は観測値の関数 $\hat{\theta}(y_1, \cdots, y_T)$ だね。

それで、「最尤推定をするためにはホワイトノイズ $\varepsilon_t$ の分布の形を特定しなければならない」？なぜですか？

例えば AR(1) 過程 $Y_t = c + \phi Y_{t-1} + \varepsilon_t$ のある1ステップだけの観測値 $y_1$ が手に入ったとき、尤度関数 $L(c, \phi, \sigma^2) = f_{Y_1}(y_1 \, ; c, \phi, \sigma^2)$ ってどうなる？過程は定常とするよ。

$y_1$ だけ観測したとき、ですか。パラメータがある $(c, \phi, \sigma^2)$ だったときにその $y_1$ がどれくらい観測されやすいのかを考えなければなりませんね。うーん…そもそも定常 AR(1) 過程には期待値があったはずです。パラメータが $(c, \phi, \sigma^2)$ であるような定常 AR(1) 過程の期待値 $E(Y_t)$ を計算して、 $y_1$ がその期待値に近ければその $(c, \phi, \sigma^2)$ の尤度は大きいでしょうし、期待値から遠ければその $(c, \phi, \sigma^2)$ の尤度は小さいのではないかと思います。定常 AR(1) 過程の期待値を復習します。AR(1) 過程の式の両辺の期待値をとります。

$\begin{split} E(Y_t) &= E(c + \phi Y_{t-1} + \varepsilon_t) \\ \Longleftrightarrow \; E(Y_t) &= c + \phi E(Y_{t-1}) \end{split}$ いま過程が定常であることを仮定していますから、期待値は時刻によらず $E(Y_t) = E(Y_{t-1})$ のはずです。これより、 $E(Y_t) = c/ (1 - \phi)$ を得ます。…AR(1) 過程の期待値はわかりましたが、この期待値に照らし合わせて $y_1$ がどれくらい観測されやすいかは、AR(1) 過程の分散にもよりますね。分散が大きければ期待値から離れた値も観測されやすいですから。分散も先ほどと同様に求められます。 $\begin{split} V(Y_t) &= V(c + \phi Y_{t-1} + \varepsilon_t) \\ \Longleftrightarrow \; V(Y_t) &= V(\phi Y_{t-1}) + V(\varepsilon_t)+ 2{\rm Cov}(\phi Y_{t-1}, \varepsilon_t) \\ \Longleftrightarrow \; V(Y_t) &= \phi^2 V(Y_{t-1}) + \sigma^2 \end{split}$ 分散もまた時刻によらず $V(Y_t) = V(Y_{t-1})$ ですから、 $V(Y_t) = \sigma^2/ (1 - \phi^2)$ を得ます。過程の期待値と分散がわかったので、これに照らし合わせて $y_1$ がどれくらい観測しやすいかを求めればいいです。

期待値と分散だけで分布の形って決まる？

！き、決まりませんでした…。で、では、歪度や尖度も求めれば…。

きりがないな…というか、定常（弱定常）って歪度や尖度まで一定とはいっていないからね。そもそもそのやり方で求まらないよ。

なんと…詰みました…。

だから扱いやすいホワイトノイズを選べばいいんだよ。正規ホワイトノイズをね。それなら任意の $t$ で $Y_t$ の分布は正規分布だよ。正規分布は再生性をもつから、正規分布にしたがうノイズを足していっても正規分布にしたがうからね。なら、期待値と分散だけで分布の形が決まるね。

確かに正規分布であれば期待値 $\mu$ と分散 $\sigma^2$ だけで確率密度関数が特定できます。

$\displaystyle f(x) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{(x-\mu)^2}{2\sigma^2}\right)$ でしたら $Y_t$ の分布は先ほど求めた期待値と分散を上式に当てはめて、以下になるはずです。
$\displaystyle f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2/ (1 - \phi^2)}} \exp \left( - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)$ これが尤度関数 $L(c, \phi, \sigma^2)$ に他なりません。観測値 $y_1$ が手に入ったとき、 $L(c, \phi, \sigma^2)$ を最も大きくする $(c, \phi, \sigma^2)$ が $y_1$ を生成した AR(1) 過程のパラメータだったのだろうと推定すればいいんですね。正規分布のお陰で詰みを回避できました…しかし、ホワイトノイズが正規分布でなかったらどうなるんです？確率密度関数が特定できなくてやはり詰みませんか？

この章では扱われていないと思う。5.2～5.6 節の全てのタイトルに明示的に Gaussian って付いてるし。正規分布じゃなかったら尤度関数をかき下すことはできないからね。じゃあ正規分布じゃなかったらどうするのかは、maximum likelihood estimates of non-Gaussian ARMA processes とかで調べたら色々文献が出てくるんじゃないかな。

ところで、いまはある1ステップだけの観測値 $y_1$ が手に入ったときの尤度関数を計算したけど、普通は連続する観測値 $(y_1, \cdots, y_T)$ からパラメータを推定するよね。じゃあ、連続する2ステップ $(y_1, y_2)$ を観測したときの尤度関数はどうなる？

2点観測したのであれば、その確率密度関数は、1点観測したときの確率密度関数を2乗すればいいのでは？

$Y_1$ と $Y_2$ は独立じゃないからね… $Y_2$ が $Y_1$ に依存しているからこそ AR（自己回帰）過程なんだけど…。

！す、すみません…AR(1) 過程は $Y_t = c + \phi Y_{t-1} + \varepsilon_t$ ですから、 $Y_2 = c + \phi Y_1 + \varepsilon_2$ ですね。 $Y_2$ はこれにしたがって生成されなければなりません。これは、 $Y_1 = y_1$ が観測されている下では、 $c + \phi y_1$ にホワイトノイズ $\varepsilon_2$ を足したものですから、

$\displaystyle f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{ \bigl( y_2 - c - \phi y_1 \bigr)^2}{2\sigma^2}\right)$ です。でしたら $(y_1, y_2)$ の同時分布は、 $f_{Y_2, Y_1}(y_2, y_1 \, ; c, \phi, \sigma^2) =f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) f_{Y_1}(y_1 \, ; c, \phi, \sigma^2)$ より、以下ですね。これが $(y_1, y_2)$ を観測したときの尤度関数です。 $\displaystyle f_{Y_2, Y_1}(y_2, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^2} \exp \left( - \frac{ \bigl( y_2 - c - \phi y_1 \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)$

うん。さらに連続する3ステップ $(y_1, y_2, y_3)$ を観測したときはこうなるね。

$\begin{split} f_{Y_3, Y_2, Y_1}(y_3, y_2, y_1 \, ; c, \phi, \sigma^2) &= f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi, \sigma^2) f_{Y_2, Y_1}(y_2, y_1 \, ; c, \phi, \sigma^2) \\ &= \; f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) \\ &= \; f_{Y_3 | Y_2}(y_3 \, | y_2 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_2 | Y_1}(y_2 \, | y_1 \, ; c, \phi, \sigma^2) \\ & \quad \cdot f_{Y_1}(y_1 \, ; c, \phi, \sigma^2) \end{split}$ 最後の式への式変形は、AR(1) 過程が1ステップ前までにしか依存しないことによったよ。具体的に確率密度関数をかき下すとこうだね。 $\displaystyle f_{Y_3, Y_2, Y_1}(y_3, y_2, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^3} \exp \left( - \sum_{t=2}^3 \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)$ これはこのまま一般化できて、 $(y_1, \cdots, y_T)$ を観測したときの確率密度関数は以下になるね。
$\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_1 \, ; c, \phi, \sigma^2) = \frac{\sqrt{1 - \phi^2}}{\left( \sqrt{2 \pi \sigma^2} \right)^T} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} - \frac{ \bigl( y_1 - c/ (1 - \phi) \bigr)^2}{2\sigma^2/ (1 - \phi^2)}\right)$ これが尤度関数 $L(c, \phi, \sigma^2)$ であり、これを最大化する $(c, \phi, \sigma^2)$ を探すことになるけど、通常はこれの対数をとった対数尤度関数 $\mathcal{L}(c, \phi, \sigma^2)$ を最大化するから、テキストではこちらで表記しているね。それで、上の式は（対数を取れば）予測誤差 $y_t - c - \phi y_{t-1}$ の2乗を足していっているようにみえる。これを尤度関数の予測誤差分解（prediction-error decomposition）というみたいだね。

おお、AR(1) 過程の尤度関数はこれで求まりましたね。…ん？ 119ページからは AR(1) 過程の尤度関数を求める別の方法が紹介されていますね？ $y \equiv (y_1, \cdots, y_T)^{\top}$ と定義するようです。このようなベクトルを定義してどうするんでしょうか…。

$y$ が多変量正規分布にしたがうと考えるんだね。期待値ベクトル $\mu$ と分散共分散行列 $\Omega$ は以下になるね。

$\displaystyle \mu = \left(\begin{array}{c} E(Y_1) \\ E(Y_2) \\ \vdots \\ E(Y_T) \end{array}\right)= \left(\begin{array}{c} c/ (1 - \phi) \\ c/ (1 - \phi) \\ \vdots \\ c/ (1 - \phi) \end{array}\right)$ $\displaystyle \Omega = \left( \begin{array}{c} V(Y_1) & {\rm Cov}(Y_1, Y_2) & \cdots & {\rm Cov}(Y_1, Y_T) \\ {\rm Cov}(Y_1, Y_2) & V(Y_2) & \cdots & {\rm Cov}(Y_2, Y_T) \\ \vdots & \vdots & \ldots & \vdots \\ {\rm Cov}(Y_1, Y_T) & {\rm Cov}(Y_2, Y_T) & \cdots & V(Y_T) \end{array} \right)= \frac{\sigma^2}{1-\phi^2}\left(\begin{array}{c} 1 & \phi & \cdots & \phi^{T-1} \\ \phi & 1 & \cdots & \phi^{T-2} \\ \vdots & \vdots & \ldots & \vdots \\ \phi^{T-1} & \phi^{T-2} & \cdots & 1 \end{array}\right)$ ちなみに AR(1) 過程の自己共分散は1時点ずれるごとに分散に $\phi$ が1つかかってくるからね。 ${\rm Cov}(Y_t, Y_{t-1}) = {\rm Cov}(c + \phi Y_{t-1} + \varepsilon_t, Y_{t-1}) = \phi {\rm Cov}(Y_{t-1}, Y_{t-1}) = \phi V(Y_{t-1})$ そうすると直接 $y$ の確率密度関数が出せる。 $\displaystyle f_Y(y \, ; c, \phi, \sigma^2) = \frac{1}{ \sqrt{ (2 \pi)^T | \Omega | } } \exp \left( - \frac{1}{2} (y - \mu)^\top \Omega^{-1} (y - \mu)\right)$ これはさっき求めた確率密度関数と一致する。計算には $\Omega = \sigma^2 V$ とおいたときに $L^{\top} L = V^{-1}$ を満たす $L$ を利用しているね。

いやいやいや、 $L$ が降って湧いていますよね？ $L^{\top} L = V^{-1}$ が成り立つことは検算して確かめよというようにありますが、この $L$ はどうやって考えついたものなんです？

うーん、 $V^{-1}$ から $L$ を求めるのは LU 分解ともちょっと違うかな？でも、 $L$ にはちゃんと意味があるよ。だって、 $\Omega^{-1} = \sigma^{-2} L\top L$ を確率密度関数に代入したら、

$\displaystyle f_Y(y \, ; c, \phi, \sigma^2) = \frac{1}{ \sqrt{ (2 \pi)^T | \Omega | } } \exp \left( - \frac{1}{2 \sigma^2} (y - \mu)^\top L^\top L (y - \mu)\right)$ となるから、 $\hat{y} = L(y-\mu)$ とおくと $\hat{y}$ は各成分が独立で各成分の分散が $\sigma^2$ になるわけだよね。 $y$ の各成分に対して上から順に、「前の成分と独立にした上で、平均を $0$ にして、分散を $\sigma^2$ にする」という操作をやっていけば $\hat{y}$ になるよ。 $L$ はその変換を行列表示したものに他ならない。AR(1) 過程は1ステップ前に依存するだけだからそんな変換は求まるし、そんな変換を求めておけば、分散共分散行列の逆行列を計算する必要もないよね。

なるほど？それで、122ページは実際に対数尤度関数の最大点を求める方法ですね。あれ？対数尤度関数そのものではなくて、 $y_1$ が所与の場合の対数尤度関数を最大化しているんですね？

その場合は最小2乗法をかき下せるからね。長い系列を観測している場合、最初の観測値の影響は無視できるし。

つづいたらつづく