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

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

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

前回： Time Series Analysis：ノート5章（その1） - クッキーの日記

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

前回の復習ですが、ノイズが正規ホワイトノイズであるような定常 AR(1) 過程 $Y_t = c + \phi Y_{t-1} + \varepsilon_t$ の実現値 $(y_1, \cdots, y_T)$ の同時分布は以下になりました。この分布は「最初に出てくる $y_1$ はどんな分布だろう」「 $y_1$ を観測した上で次の $y_2$ はどんな分布だろう」と積み上げていけば導くことができます。

定常ガウシアン AR(1) 過程の実現値の同時分布

$\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)$ 実現値 $(y_1, \cdots, y_T)$ からパラメータ $(c, \phi, \sigma^2)$ を推定するには上の同時分布が最大になるような $(c, \phi, \sigma^2)$ を選べばいいわけです（実際にはこれの対数を取って最大化しますが）。これを最尤推定といいます。

119~122ページでは上の同時分布を求める別解として、「 $y \equiv (y_1, \cdots, y_T)^{\top}$ はどんな多変量正規分布にしたがうだろうか？」というアプローチで同じ同時分布を導きました。

122~123ページでは「 $y_1$ が所与」という考えやすい場合で実際に最尤推定値を求めますが、前回は流してしまったので改めて追ってみます。 $y_1$ が所与なら上の同時分布は以下のように化けますね。

定常ガウシアン AR(1) 過程の実現値の同時分布（ $y_1$ が所与の場合）

$\displaystyle f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_2 \, ; y_1, c, \phi, \sigma^2) = \frac{1}{\left( \sqrt{2 \pi \sigma^2} \right)^{T-1}} \exp \left( - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2} \right)$ 対数を取るとこうですね。 $\displaystyle \mathcal{L}(c, \phi, \sigma^2) \equiv \log f_{Y_T, \cdots, Y_1}(y_T, \cdots, y_2 \, ; y_1, c, \phi, \sigma^2) = - \frac{T-1}{2} \log \left( 2 \pi \sigma^2 \right) - \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^2}$ これを $(c, \phi, \sigma^2)$ について最大化したいわけですから、これらで偏微分してみましょう。 $\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial c} = \sum_{t=2}^T \frac{ 2 \bigl( y_t - c - \phi y_{t-1} \bigr)}{2\sigma^2}$ $\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial \phi} = \sum_{t=2}^T \frac{ 2 y_{t-1} \bigl( y_t - c - \phi y_{t-1} \bigr)}{2\sigma^2}$ $\displaystyle \frac{\partial \mathcal{L}(c, \phi, \sigma^2)}{\partial \sigma^2} = - \frac{T-1}{2\sigma^2} + \sum_{t=2}^T \frac{ \bigl( y_t - c - \phi y_{t-1} \bigr)^2}{2\sigma^4}$ これらの式をイコールゼロとしたものを満たす $(c, \phi, \sigma^2)$ が最尤推定値 $(\hat{c}, \hat{\phi}, \hat{\sigma}^2)$ になりますが、まず上の2式より $(\hat{c}, \hat{\phi})$ が導かれ、次に最後の式より $\hat{\sigma}^2$ が導かれますね。

$y_1$ が所与の定常ガウシアン AR(1) 過程の最尤推定値 $(\hat{c}, \hat{\phi}, \hat{\sigma}^2)$

$\displaystyle \left(\begin{array}{cc} T-1 & \sum_{t=2}^T y_{t-1} \\ \sum_{t=2}^T y_{t-1} & \sum_{t=2}^T y_{t-1}^2 \end{array}\right) \left(\begin{array}{c} \hat{c} \\ \hat{\phi} \end{array}\right) = \left(\begin{array}{c} \sum_{t=2}^T y_{t-1} \\ \sum_{t=2}^T y_{t-1}y_t \end{array}\right)$ $\displaystyle \hat{\sigma}^2 = \frac{1}{T-1} \sum_{t=2}^T \bigl( y_t - \hat{c} - \hat{\phi} y_{t-1} \bigr)^2$ $y_1$ が所与なら最尤推定値を求めることは容易ですね。ただ、「言い換えると、これは最小2乗回帰（OLS）である」とありますが、どういうことなんでしょう？

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

いま、以下の前者の問からその最尤推定値にたどり着いたよね。でもそれって、以下の後者の問から出発したのと同じになってるんだよね。

$y_1$ が所与の定常ガウシアン AR(1) 過程の、実現値 $(y_2, \cdots, y_T)$ に対する最尤推定値は何だろうか？
$y_1$ が所与の AR(1) 過程の、実現値 $(y_2, \cdots, y_T)$ に対する最小2乗推定値は何だろうか？

だって $(\hat{c}, \hat{\phi})$ を求めるのに $\bigl( y_t - c - \phi y_{t-1} \bigr)^2$ を最小化しているけど、これは各 $y_t$ を $c + \phi y_{t-1}$ で予測したときの2乗誤差に他ならないし、 $\hat{\sigma}^2$ は残差の分散になっているしね。ちなみに、後者は定常性とガウス性を仮定していない。「誤差2乗和を最小化し、残差の分散を求める」という作業は定常性もガウス性も要求しないからね。

あっと、最小2乗推定と等価ということはわかりました。しかし、なら最初から最小2乗推定を紹介していただければよかったのではありませんか？定常性とガウス性を仮定して頑張って尤度関数を求めて、「実は単に最小2乗推定するのと同じだったんだけどね」では私が報われなくありませんか？

そんなことないよ。確かに最小2乗法はいつでもできるよ。最小2乗推定値自体はいつでも求まると思う。でも、それがまともな推定値になっていることを支持する材料ってないからね。定常性とガウス性をもつ場合には最小2乗推定値は最尤推定値になっていることがいまわかったけど、逆にいうとそうでない場合は変な値になっているかもしれないってことだからね。まあそう言い出すと「なら最尤推定値ならまともなのか」って話になるけど、単なる最小2乗推定値よりはまともだと信じられる場面が多いってことだと思ってる。

な、なるほど…？それで、123ページからは AR(p) 過程の場合ですね。AR(p) 過程の同時分布は…

まず「最初の $p$ 点 $(y_1, \cdots, y_p)^{\top}$ はどんな多変量正規分布にしたがうだろうか？」を考えた上で、
その後「 $(y_1, \cdots, y_p)$ を観測した上で次の $y_{p+1}$ はどんな分布だろう」を積み上げていっています。

AR(1) 過程でみた2つアプローチの合わせ技をやっていますね。素朴なアプローチではできないんでしょうか？

…5.3.8 式に AR(2) の場合が具体的に書き下してあるから、AR(2) で求めて検算してみればいいんじゃないかな…と思ったけど、素朴にはできない気がするかな。まず、何も観測していない下での最初の1点の分布 $f_{Y_1}(y_1)$ は AR(1) のときと一緒だよね。期待値と分散は AR(2) のそれにしなければならないけど。あと、直近の $p$ 点 $(y_1, y_2)$ を観測した下での $y_3$ の分布 $f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1)$ も AR(1) のときと同じ要領だね。

$\displaystyle f_{Y_1}(y_1 \, ; c, \phi_1, \phi_2, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2/ (1 - \phi_1^2 - \phi_2^2)}} \exp \left( - \frac{ \bigl( y_1 - c/ (1 - \phi_1 - \phi_2) \bigr)^2}{2\sigma^2/ (1 - \phi_1^2 - \phi_2^2)}\right)$ $\displaystyle f_{Y_3 | Y_2, Y_1}(y_3 \, | y_2, y_1 \, ; c, \phi_1, \phi_2, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( - \frac{ \bigl( y_3 - c - \phi_1 y_2 - \phi_2 y_1 \bigr)^2}{2\sigma^2}\right)$ でもこれらの中間の $f_{Y_2 | Y_1}(y_2 \, | y_1)$ がこんな感じじゃ求められない気がするんだよね。だって、「何も観測していない下での」でもないし、「直近の $p$ 点を観測した下での」でもないから、どっちのアプローチも使えない。結局 AR(2) の共分散構造を満たすように $(y_1, y_2)$ の同時分布を真面目に求めないといけない気がする。その分散共分散行列が125ページの一番上の式だと思う。

素朴にはいかなさそうなんですね。…あれ、126ページには「ガウス性をもたない場合の最尤推定」とありますね。5章ではガウス性をもたない場合は扱わないはずでは？

…ここにかいてあるのは、ガウス性をもたない過程にガウス性をもつと仮定した最尤推定をすることを「準最尤推定」ということ。生データがガウス性をもたない場合でもシンプルな変換でガウス性をもつようにできることがあること。とかかな。その変換の例が 126ページの下の方にある Box-Cox 変換だね。

この変換は $\lambda$ というパラメータを与えて用いるんですね。どのように決定すればいいんでしょうか。

無論ガウス分布に近くなるように決めるんだけど、scipy.stats.boxcox は $\lambda$ を与えなければ勝手に決めてくれるって。

scipy.stats.boxcox — SciPy v1.5.4 Reference Guide

文明の利器…127ページの一番上には非ガウスの場合の ARMA 過程の推定に関する研究が紹介されていますね。この本自体が17年前の本なのでそれを考慮して参考にしなければなりませんが…。5.4~5.6節は MA 過程、ARMA 過程の話ですが一旦とばします。5.6 節の最後には、ガウシアン ARMA 過程を最尤推定するのに最もシンプルなアプローチはカルマンフィルタであるとありますね。また、他の参考文献も紹介されています。

それで、尤度 $\mathcal{L}(\theta)$ を出す式が特定できた下で、それを最大化する $\hat{\theta}$ を如何に求めるかという話が残っていますが（ $y_1$ が所与の場合の最小2乗法はやりましたが）、5.7 節は数値的に最大化する方法ですね。紹介されているのは以下でしょうか。

グリッドサーチ
最急降下法
Newton-Raphson 法
Davidon-Fletcher-Powell 法（準 Newton 法）

それで、5.8 節が最尤推定を用いた統計的検定で、5.9 節の「不等式制約」というのは、パラメータ探索は「過程が定常である」「分散は正である」のような制約の下で行わなければならないので、そのために使えるテクニックですかね。

つづいたらつづく