2020/05/22

時系列分析 (3) - モデルの推定と診断

いよいよ、時系列データからモデルを作ることを考えます。モデルを推定し、そのモデルの善し悪しを診断する方法を紹介します。

モデルの推定

時系列モデルの推定とは、時系列データがどのような種類の確率過程に従うか、そしてそのパラメータは何かを推定することです。モデルが推定できれば、未来の時点の値の予測に用いることができます。
モデルを推定する方法としては、最小二乗法 (Ordinary least squares/OLS) と、最尤推定 (Maximum likelihood estimation/MLE) が良く知られています。最小二乗法はシンプルかつ良い推定量を導出できるのでよく用いられますが、複雑なモデルでは適用できないことがあり、その場合は最尤推定を用いることになります。

最小二乗法

最小二乗法は、モデルによって説明できない項 (誤差項) を最小にするようにパラメータを推定する方法です。
例えば、ある過程$y_t$を$AR(1)$モデル$y_t = c + \phi y_{t-1} + \epsilon_t, \epsilon_t \sim WN(\sigma^2)$に当てはめることを考えましょう。$c'$と$\phi'$がそれぞれ$c$と$\phi$の推定値だとすると、このモデルで説明できない値$e_t$は$e_t = y_t - c' - \phi'y_{t-1}$となります。これを残差 (residuals) と呼びます。この残差の絶対値を小さくすれば、良いモデルに近づくことになります。最小二乗法では、残差平方和 (Sum of squared residuals/SSR) を最小にするようなパラメータを選択します。このモデルでは、残差平方和は$SSR = \sum_{t=2}^T(y_t - c' - \phi'y_{t-1})^2$となります。この$SSR$を最小化するような$c', \phi'$ (これを OLS 推定量と呼び、ここでは$c'', \phi''$ で表すことにする) を求めれば良いことになります。そのため、$c', \phi'$それぞれで偏微分します。
$$
\begin{eqnarray*}
\frac{\partial SSR}{\partial c'} &=& -2\sum_{t=2}^T (y_t - c' - \phi'y_{t-1}) \\
\frac{\partial SSR}{\partial \phi'} &=& -2\sum_{t=2}^T y_{t-1}(y_t - c' - \phi'y_{t-1})
\end{eqnarray*}
$$
この2式が$0$になるような$c', \phi'$が$c'', \phi''$となるので、$c'', \phi''$は以下の2式を満たします。
$$
\begin{eqnarray*}
\sum_{t=2}^T (y_t - c'' - \phi''y_{t-1}) = 0 \\
\sum_{t=2}^T y_{t-1}(y_t - c'' - \phi''y_{t-1}) = 0
\end{eqnarray*}
$$
これを正規方程式と呼び、解くと以下のようになります。ただし、$\bar{y}_{a,b} = \frac{1}{b-a+1}\sum_{t=a}^b y_t$ ($y_a \sim y_b$の平均)です。
$$
\begin{eqnarray*}
c'' &=& \bar{y}_{2,T} - \phi''\bar{y}_{1,T-1} \\
φ'' &=& \frac{\sum_{t=2}^T (y_t - \bar{y}_{2,T})(y_{t-1} - \bar{y}_{1,T-1})}{\sum_{t=2}^T (y_{t-1} - \bar{y}_{1,T-1})^2}
\end{eqnarray*}
$$
誤差項の分散$\sigma^2$の OLS 推定量$s^2$は、OLS 残差の不偏標本分散なので、以下のようになります。
$$
s^2 = \frac{1}{T - 2}\sum_{t=2}^T (y_t - c'' - \phi''y_{t-1})^2
$$
OLS 推定量は、不偏推定量 (平均的に真の値をとらえる)、一致推定量 (標本数が大きくなれば真の値に限りなく近づく) の性質を持ちます。また、線形不偏推定量の中で最小の分散共分散行列を持ち、$\epsilon_t \sim iidN(\sigma^2)$であれば任意の不偏推定量の中で最小の分散共分散行列を持ちます。つまり、OLS はこういった条件の下で最適性を持っています。OLS 推定量を基準化すると、正規分布に漸近するため、仮説検定ではこの性質も活用します。

最尤法

最尤法は、モデルが観測値を実現しやすくなるようにパラメータを選択する方法です。もう少し具体的に言うと、観測値が実現される確率を尤度関数 (または単に尤度) というパラメータの関数として見て、この尤度関数を最大化するという方法です。尤度関数を最大化するパラメータを最尤推定値と呼びます。最尤推定値は解析的に求めることができる場合もありますが、解析解が無く数値的に求めなければならない場合もあります。解析的に求める場合は尤度関数を微分することになりますが、尤度関数の最大化と尤度関数の対数を取ったもの (対数尤度と言う) の最大化は同じ結果となるので、微分しやすい対数尤度を用いることが多いです。

例1 簡単な確率の推定

例えば簡単な例として、「サイコロを10回振ったら、1が2回出た」という観測値があったとします。ここから、1が出る確率$p$を最尤法で推定してみましょう。

尤度関数$L(p)$は、「1が出る確率が$p$のときに、10回中2回1が出る確率」ですから、$L(p) = {}_{10}\rm{C}_2p^2(1-p)^8$となります。このままでは微分しづらいので対数をとると、$\log L(p) = \log({}_{10}\rm{C}_2p^2(1-p)^8) = \log {}_{10}\rm{C}_2 + 2 \log p + 8 \log(1-p)$となります。
微分すると、$$\frac{d}{dp} \log L(p) = \frac{2}{p} - \frac{8}{1-p} = \frac{2 - 10p}{p(1-p)}$$なので、これが$0$となるような$p$が最尤推定値となります。$2 - 10p = 0$すなわち$p = 0.2$となります。

例2 AR(1) のパラメータ推定

最小二乗法で例示したものと同様の$AR(1)$のパラメータ推定を最尤法で行ってみましょう。
今回は$y_t = c + \phi y_{t-1} + \epsilon_t, \epsilon_t \sim iidN(0, \sigma^2)$とします。最小二乗法の時は$\epsilon_t \sim WN(\sigma^2)$としていましたが、$\epsilon_t$の分布が分からないと尤度関数が定まらないので、一般の$WN$ではなく、$iidN$としています。

尤度関数は、同時密度を条件付き確率に分解する$P(X, Y) = \frac{P(X | Y)}{P(Y)}$を用いると、確率密度関数$f()$を使って$$L(c, \phi, \sigma^2) = f_{Y_T, Y_{T-1}, \ldots, Y_2|Y_1}(y_T, y_{T-1}, \ldots, y_2|y_1; c, \phi, \sigma^2) = \prod_{t=2}^T f_{Y_t|Y_{t-1}}(y_t|y_{t-1}; c, \phi, \sigma^2)$$と書けます。対数をとると、$$\log L(c, \phi, \sigma^2) = \sum_{t=2}^T \log f_{Y_t|Y_{t-1}}(y_t|y_{t-1}; c, \phi, \sigma^2)$$となります。ここで$\epsilon_t \sim iidN(0, \sigma^2)$を用いると、$y_t|y_{t-1} \sim N(c + \phi y_{t-1},\sigma^2)$より$$
f_{Y_t|Y_{t-1}}(y_t|y_{t-1}; c, \phi, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}exp\left(-\frac{(y_t-(c+\phi y_{t-1}))^2}{2\sigma^2}\right)
$$なので、対数尤度は$$
\log L(c, \phi, \sigma^2) = -\frac{T}{2}log 2\pi - \frac{T}{2}\log\sigma^2 - \sum_{t=2}^T\left(\frac{(y_t-c-\phi y_{t-1})^2}{2\sigma^2}\right)
$$となります。
この対数尤度の最大化は、$$\sum_{t=2}^T (y_t - c - \phi y_{t-1})^2$$の最小化になりますが、これは最小二乗法で見た残差平方和そのものです。つまり、最尤推定量は OLS 推定量に一致することになります。

なお、最尤推定量は不偏推定量とは限りません。例えば、$\sigma^2$の最尤推定量は、$\sigma^2 = \frac{1}{T}\sum_{t=2}^T(y_t - c - \phi y_{t-1})^2$ですが、これは OLS 推定量とは一致しません。OLS 推定量の方が不偏推定量で、最尤推定量は不偏推定量ではないのです。最尤推定量は一致推定量ではあります。また基準化すると正規分布に漸近します。最尤推定量は、一致推定量の中で漸近的に最小の分散共分散行列を持ちます。

偏自己相関

パラメータ推定以前の問題として、そもろもモデルとして MA, AR, ARMA のどれが適切であるかを選択する必要があります。これらのモデルの違いは、自己相関の組み込み方の違いですから、選択にあたってはまずは標本自己相関関数を見て判断することになります。ある次数$q$から先が$0$に近いようであれば、$MA(q)$で当てはめられる可能性が高いわけです。しかしそうではない場合、AR 過程とみるか、ARMA 過程とみるかは、標本自己相関関数からの判断は難しくなってきます。この判断にあたっては、偏自己相関を用いるのが適切です。
$k$次-自己相関は、$y_t$と$y_{t-k}$の相関でしたが、$k$次-偏自己相関は、$y_t$と$y_{t-k}$の相関からその間にある$y_{t-1}, \ldots, y_{t-k+1}$の影響を除去したものです。自己相関関数と同様、$k$の関数として見て、偏自己相関関数とも言います。

自己相関関数は、$AR(p), ARMA(p, q)$では減衰、$MA(q)$では$q$次を越えると$0$になる性質がありましたが、これに対して偏自己相関関数は、自己回帰項の影響が取り除かれることによって、$AR(p)$では$p$次を越えると$0$になる性質があります。$MA(q), ARMA(p, q)$では減衰します (MA は$AR(\infty)$で書き直せることからわかります)。

従って、標本自己相関関数と標本偏自己相関関数を見ることによって、MA, AR, ARMA のいずれが適しているかを見定めることができます。MA, AR であれば切断箇所から次数も見積もることができます。ARMA の次数は自己相関、偏自己相関からは分かりません。またもちろん、標本自己相関関数、標本偏自己相関関数には真の値との間に誤差があることも忘れてはなりません。

情報量規準

あるモデル1つを仮定すれば、そのパラメータの推定量を最尤法などで求めることはできますが、そもそもどのモデルが最適なのかを決める術が必要です。その道具となるのが情報量規準です。
情報量規準は、モデルの当てはまりの良さと、モデルの複雑さに対するペナルティから構成されています。情報量規準$IC$は以下のような一般形で表されます。$$
IC = -2\log L(\Theta) + p(T)k
$$
$\Theta$は最尤推定値ベクトルで、$L(\Theta)$は最大尤度、$k$はパラメータ数、$p(T)$はペナルティ関数です。つまり第1項がモデルの当てはまりの良さを表現しており、第2項が複雑さに対するペナルティを表しています。$IC$の小さい方が、良いモデルであると考えます。

AICとSIC

ペナルティ関数に選び方によって$IC$の値は変わってきます。良く使われるのが赤池情報量規準 ($AIC$) と Schwarts 情報量規準 ($SIC$) です。$SIC$はベイズ情報量規準 ($BIC$) とも呼ばれています。$AIC$は$p(T) = 2$、$SIC$は$p(T) = \log T$としたもので、全体としては以下のようになります。$$
\begin{eqnarray*}
AIC &=& -2\log L(\Theta) + 2k \\
SIC &=& -2\log L(\Theta) + k\log T
\end{eqnarray*}
$$
$T \ge 8$では$SIC > AIC$です。通常は$T \ge 8$でしょうから、$SIC$は$AIC$よりもペナルティが大きいことになります。すなわち、$SIC$に基づく方がより小さいモデル (パラメータの少ないモデル) を選択する傾向があるということになります。また、$AIC$は$T$を大きくしても、パラメータの多すぎるモデルを選択してしまう確率が$0$になりませんが、$SIC$は$T$を大きくするほど最適なモデルを選択する確率が上がっていきます。ただし、$AIC$が過大なパラメータ数を持つモデルを選択しても、多すぎる部分のパラメータの最尤推定量は$0$に漸近するので、$AIC$も十分に有用です。

モデルの診断

最適とされるモデルが選択されたら (あるいは候補が絞り込まれたら)、そのモデルが適切かどうかの診断を行います。この診断で適切と判断されなければ、これまで候補に挙がったモデルの中に適切なモデルが無かったということになるので、おそらくはモデルの候補の範囲を広げて推定からやり直すことになるでしょう。

例えば真のモデルが MA, AR, ARMA 等だとすると、$\epsilon_t \sim WN(\sigma^2)$ですから、$\epsilon_t$は自己相関を持たないはずです。そこで、最適とされたモデルにおける残差が自己相関を持たないことをかばん検定を用いて検定するという方法が使えます。ただし、パラメータ数の分だけ自由度が下がることを考慮し、$Q(m)$と$\chi^2(m-k)$ ($k$はパラメータ数) とで検定を行います。