時系列モデルのパラメーター推定

時系列モデルを考えた時、推定しないと値が分からないパラメーターがいくつかあります。今回の記事では、パラメーターの推定方法として、最小二乗法と最尤推定を紹介します。
その後、python で計算させてみて、正しく推定できるか確認します。

一応3冊が参考文献です。

最小二乗法

ARモデルのパラメーターを推定します。一つの手法として、回帰分析でもおなじみの最小二乗法があります。
AR(p)モデルは以下の式です。
$$\begin{eqnarray}
y_t&=& c+ \sum_i ^p\phi _i y_{t-i} +\epsilon _t \\
\epsilon _t &\sim & { \rm W.N. } (\sigma ^2 )
\end{eqnarray}$$
\( y_t \)を 特徴量\(y_{t-1} , \cdots , y_{t-p} \)から予測すると思うと、上式は重回帰分析の式と同じ形です。回帰分析の記事があるので見てみてください。
https://masamunetogetoge.com/single-regression

重回帰分析

変数が沢山ある場合の重回帰分析の解説をします。行列とベクトルを用いて式を書く事で、シンプルに結論までたどり着くことが出来ます。また、多重共線性という概念が自然に出てくることを見ます。

masamunetogetoge.com

2019.06.13

\(AR(1) \)の場合で、\( \phi , c , \sigma \)を最小二乗法で予測しましょう。\(\epsilon _t \)を誤差項だと思い、最小化します。
$$\begin{eqnarray}
\epsilon _t &=& y_t- c- \phi y_{t-1} \\
E_t &=& ( \epsilon _t )^2 \\
E&=& \sum _{t=1}^T E_t
\end{eqnarray}$$
\(E \)を\(c , \phi \)の関数だと思って偏微分して0とおき、\( c, \phi \)について解く事でパラメーターの推定値が求められます。
$$\begin{eqnarray}
\frac{ \partial E}{\partial c} &=& \sum (y_t -c – \phi y_{t-1} )=0 \\
\frac{ \partial E}{\partial \phi } &=& \sum y_{t-1} (y_t -c – \phi y_{t-1} )=0
\end{eqnarray}$$
この2つの式と、平均値を使った2乗の式¹から、以下のように推定値が計算出来ます。
$$\begin{eqnarray}
\hat{c} &=& \bar{y_t} – \overline{y_{t-1}} \hat{\phi} \\
\hat{\phi} &=& \frac{ \sum (y_t – \bar{y_t} )(y_{t-1} – \overline{y_{t-1}} ) }{\sum ( y_{t-1} – \overline{y_{t-1}} )^2}
\end{eqnarray}$$
ただし、\(\bar{y_t} =\sum_{t=1}^T y_t /T, \overline{y_t} =\sum_{t=1}^T y_{t-1}/T \)です。

\(AR(p) \)モデルについても、重回帰分析と全く同様に計算出来ます。

最尤法

最尤法でARモデルのパラメーターを推定します。
その為に、 ホワイトノイズは正規分布に従っているとします。つまり、全ての\(t \)について
$$\begin{eqnarray}
\epsilon _t \sim \mathcal{N} (0, \sigma ^2 )
\end{eqnarray}$$
とします。もちろん、独立性を仮定します。こうすると、パラメーターは\(c, \phi , \sigma ^2 \)になります。最尤法を使うにあたって、初期値\(y_0 \)が得られているとします。ベイズの定理を何度も使う事で、
$$\begin{eqnarray}
p(y_T, y_{T-1}, \cdots , y_0 ) = \prod_{i=1}^{T} p(y_{t_i}|y_{t_{i-1}})
\end{eqnarray}$$
が成り立ちます。²
また、\( \epsilon _t \sim \mathcal{N} (0, \sigma ^2 ) \)と、\( y_t= c+ \phi y_{t-i} +\epsilon _t \) から³、
$$\begin{eqnarray}
p(y_t | y_{t-1} )= \mathcal{N}(y_t | c+ \phi y_{t-i} , \sigma ^2 )
\end{eqnarray}$$
です。正規分布に従う確率変数の線形変換については以下の記事を参考にしてください。

正規分布の性質

正規分布の導入をします。確率になっていること、平均の値、分散の値を具体的に計算します。

masamunetogetoge.com

2019.06.26

対数を取っても関数の勾配の正負は変わらないので、尤度より計算が楽な対数尤度を使う事が多いです。
$$\begin{eqnarray}
L &=& \sum _{t=1}^T \log p(y_t | y_{t-1} ) \\
&=& -\frac{T}{2} \log 2\pi – \frac{T}{2} \log \sigma ^2 – \sum \frac{(y_t – c – \phi y_{t-1} )^2}{2\sigma ^2}
\end{eqnarray}$$
\(c, \phi , \sigma ^2 \) で偏微分して0と置くことで、パラメーターの推定値が分かります。計算すると、\(c, \phi \)は最小二乗法と同じ値になることが分かります。
$$\begin{eqnarray}
\frac{\partial L}{\partial \sigma ^2 } &=& – \frac{T}{2\sigma ^2} + \sum \frac{(y_t – c – \phi y_{t-1} )^2}{2\sigma ^4}
\end{eqnarray}$$
なので、\(\sigma ^2 \)の推定値は以下のようになります。
$$\begin{eqnarray}
\hat{\sigma }^2 &=& \frac{1}{T} \sum_{t=1}^{T} \left( y_t – c- \phi y_{t-1} \right) ^2
\end{eqnarray}$$
計算を見て分かるように、\(AR(p) \)モデルの時の\( \sigma ^2 \)の推定値は、以下のようになります。

MAモデル

MAモデルの場合は、時系列データ同士が直接等式で結ばれていないので、少し複雑です。
最尤法を考えましょう。ARモデルの時と同ような事をして、対数尤度は
$$\begin{eqnarray}
L &=& \sum _{t=2}^T \log p(y_t | y_{t-1}、\cdots , y_1 ) +\log p(y_1 )
\end{eqnarray}$$
と表せます。\(MA(p) \) は複雑なので、\( MA(1) \)モデルを考えます。⁴
$$\begin{eqnarray}
y_t = c+\epsilon _t + \theta \epsilon_{t-1}
\end{eqnarray}$$
ホワイトノイズは、独立な正規分布に従うとしましょう。また、\( \epsilon _0 =0 \)が初期値として得られていると仮定します。このとき、
$$\begin{eqnarray}
y_1 =c+ \epsilon _1
\end{eqnarray}$$
より、
$$\begin{eqnarray}
p(y_1) = \mathcal{N} (y_1 | c, \sigma ^2 )
\end{eqnarray}$$
が成り立ちます。また、
$$\begin{eqnarray}
y_2 =c+ \epsilon _2 + \epsilon _1 = c+\epsilon _2 +\theta y_1 -\theta c
\end{eqnarray}$$
です。これから、\(y_1 \)が分かっている時は、
$$\begin{eqnarray}
p(y_2 |y_1 ) = \mathcal{N} (y_ | \theta y_1 +c-\theta c , \sigma ^2 )
\end{eqnarray}$$
となります。これを続けると、
$$\begin{eqnarray}
y_n &=& \epsilon _n +\theta \sum_{i=0}^{n-2} \left((-\theta )^i (y_{n-i-1} -c) \right) +c \\
p(y_n|y_{n-1} , \cdots , y_0 ) &=& \mathcal{N}(y_n| \theta \sum_{i=0}^{n-2} \left((-\theta )^i (y_{n-i-1} -c) \right) +c , \sigma ^2 )
\end{eqnarray}$$
となることが分かります。例えば、\(n=3 \)では、
$$\begin{eqnarray}
p(y_3|y_2 , y_1 )= \mathcal{N}(y_n| \theta y_2 – \theta ^2 y_1 – \theta c + \theta ^2 \mu +c, \sigma ^2 )
\end{eqnarray}$$
となります。 対数尤度は以下のようになります。
$$\begin{eqnarray}
L&=& -\frac{T}{2} \log 2\pi – \frac{T}{2} \log \sigma ^2 -\sum_{t=2}^T \left( y_t -( \theta \sum_{i=0}^{t-2} (-\theta )^i (y_{t-i-1} -c ) +c) \right) ^2 /2\sigma ^2-\frac{(y_1 -c)^2 }{2\sigma ^2 }
\end{eqnarray}$$
これから、各パラメーターで偏微分して0とおくことで、パラメーターの推定値を求めることが出来ます。
綺麗にパラメーターを求める方法があれば教えてください。

pythonによるパラメーター推定

ARモデルを適当に生成して、最尤法と最小二乗法でパラメーターを求めてみます。

import numpy as np
import pandas as pd
np.random.seed(seed=1)
#最小二乗法
def OLS(Y):
    Y_1 =pd.Series(Y).shift(1).fillna(0)
    phi = np.dot(Y-np.mean(Y), Y_1 -np.mean(Y_1) )/ np.sum( (Y_1 -np.mean(Y_1))**2)
    c = np.mean(Y) - phi*np.mean(Y_1)
    return phi, c
#最尤法
def MLM(Y):
    phi_hat, c= OLS(Y)
    Y_1  =pd.Series(Y).shift(1).fillna(0)
    sig =np.sum( (Y-c- phi_hat*Y_1)**2 )/len(Y)
    return sig

\( \phi =0.6 , c=0, \sigma ^2 =1 \)の\(AR(1) \)モデルを生成します。データの数が多くなると予測が正確になるはずなので、データの数を 増やした時に正確な予測が出来るか調べてみます。

for T in [100,1000,10000]:
    Y=np.zeros(T)
    for t in range(1,T):
        Y[t]=phi*Y[t-1] +epsilon[t]
    phi_hat,c =OLS(Y)
    sig = MLM(Y)
    print("T={}".format(T))
    print("phi = {}".format(phi_hat))
    print("c = {}".format(c,T))
    print("σ^2 = {}".format(sig))
"""
T=100
phi = 0.6015330784508672
c = 0.03145012413073635
σ ^2 = 1.1263396263956051
T=1000
phi = 0.6018348509325396
c = 0.05636292389898365
σ ^2 = 1.0478865879775259
T=10000
phi = 0.611607654506957
c = -0.00029865192568464345
σ ^2 = 1.0154715330598576
"""

データ数\(T \)が大きくなるにつれて、\( c , \sigma ^2 \)は正確に予測できるようになっています。尤度を比べてみましょう。尤度の式を再渇します。
$$\begin{eqnarray}
L =-\frac{T}{2} \log 2\pi – \frac{T}{2} \log \sigma ^2 – \sum \frac{(y_t – c – \phi y_{t-1} )^2}{2\sigma ^2}
\end{eqnarray}$$
この式をそのまま使うと、\(T \)が大きくなった時にペナルティが付いてる感じなので、\(T \)で割った式を使います。
$$\begin{eqnarray}
L’=L/T
\end{eqnarray}$$
\(L’ \)が小さくなってるはずです。

#尤度の計算
def Likelihood(Y):
    T=len(Y)
    phi_hat , c = OLS(Y)
    sig = MLM(Y)
    Y_1 =pd.Series(Y).shift(1).fillna(0)
    x_t = (Y-c-phi_hat *Y_1)**2
    L= -T/2 *np.log(2*np.pi) - T/2 * np.log(sig) - np.sum(x_t**2)/(2*sig)
    return L

for T in [100,1000,10000]:
    Y=np.zeros(T)
    for t in range(1,T):
        Y[t]=phi*Y[t-1] +epsilon[t]
    L=Likelihood(Y)/len(Y)
    print("T={}".format(T))
    print("L = {}".format(L))

"""
T=100
L = -2.8361167328901393
T=1000
L = -2.5598749200521143
T=10000
L = -2.443939503925576
"""

尤度が増えている事が確認できました。パラメーターの推定としての精度は上がっているようです。

まとめ

• ARモデルのパラメーター推定の方法として、最小二乗法を紹介した
• 最尤法を紹介した
• MAモデルの計算が出来なかった
• python でパラメーターの推定の計算をした

1. \( \sum(y_t – \bar{y_t} )^2 = \sum (y_t )^2 -T\bar{y_t} ^2 \) などです。
2. ベイズの定理を忘れた人は解説記事をどうぞ
3. \(y_t \)たちは、\(\epsilon _t \)によって確率変数になります。
4. 計算が間違っていたらすいません。