【Rで時系列解析】最小２乗法で自己回帰過程の当てはめ - ケィオスの時系列解析メモランダム

観測された時系列に最小２乗法で自己回帰過程を当てはめる計算のメモです．どんな形の過程を仮定しても，基本的な考え方を理解しておけば，係数を求めることができます．楽をしたければ，arとか，varsパッケージなどを使えば良いと思います．しかし，お気楽なコマンドに頼っていると，理屈や計算が難しいそうという先入感で思考が停止し，問題解決のための応用パターンや拡張パターンの発見を見逃してしまうかもしれません．ここでは，モデルの形を自由に決めたいし，理解を深めたいので，行列の計算でモデルのパラメタを推定してみます．以下に，いくつかのパターンを示しておきます．

単変量の自己回帰過程
- １次自己回帰過程
- 2次自己回帰過程
多変量の場合 (完全な自己回帰過程ではありません)
- 自己ともう一つ別の変数に依存する過程
- 瞬時的因果がある過程
まとめ
参考にしたページ

単変量の自己回帰過程

１次自己回帰過程

　時系列 $y_t$ が次の１次自己回帰過程に従うとします．

　 $y_t = a_0 + a_1 y_{t-1} + \varepsilon_t$

観測データとして長さ $N$ のデータ $\{y_t\}$ がある場合に最小２乗法で係数 $\{a_0, a_1\}$ の推定値 $\{\hat{a}_0, \hat{a}_1\}$ をRで計算してみます．データ $\{y_t\}$ を仮定の式に代入すると，

　 $\begin{eqnarray} y_N &=& a_0 + a_1 y_{N-1} + \varepsilon_N \\ y_{N-1} &=& a_0 + a_1 y_{N-2} + \varepsilon_{N-1} \\ &\vdots& \\ y_{2} &=& a_0 + a_1 y_{1} + \varepsilon_{2} \\ \end{eqnarray}$

のように， $N-1$ 個の式になります．これを行列で表せば，

　 $\begin{bmatrix} y_N \\ y_{N-1} \\ \vdots\\ y_{2} \\ \end{bmatrix} = \begin{bmatrix} 1 & y_{N-1} \\ 1 & y_{N-1} \\ \vdots & \vdots \\ 1 & y_{1} \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ \end{bmatrix} + \begin{bmatrix} \varepsilon_{N} \\ \varepsilon_{N-1} \\ \vdots \\ \varepsilon_{2} \\ \end{bmatrix}$

となります．さらに，

　 $\vec{y} = \begin{bmatrix} y_N \\ y_{N-1} \\ \vdots\\ y_{2} \\ \end{bmatrix}, \quad X = \begin{bmatrix} 1 & y_{N-1} \\ 1 & y_{N-2} \\ \vdots & \vdots \\ 1 & y_{1} \\ \end{bmatrix}, \quad \vec{a} = \begin{bmatrix} a_0 \\ a_1 \\ \end{bmatrix} , \quad \vec{\varepsilon} = \begin{bmatrix} \varepsilon_{N} \\ \varepsilon_{N-1} \\ \vdots \\ \varepsilon_{2} \\ \end{bmatrix}$

で表せば，

　 $\vec{y} = X \vec{a} + \vec{\varepsilon}$

です．ベクトルはすべて縦書きの列ベクトルとします．残差としての $\{\varepsilon_{t}\}$ の2乗和が最小になるように (最小２乗法で) $\vec{a }$ を計算すると，

　 $\vec{a } = \left({}^t X X \right)^{-1} {}^t X \, \vec{y}$

となります． ${}^t X$ は， $X$ の転置行列を表します．ということで，今回のポイントは，

$\vec{y} = X \vec{a} + \vec{\varepsilon}$ で ${}^{t} \vec{\varepsilon}\, \vec{\varepsilon}$ が最小になるとき，係数ベクトル $\vec{a}$ は，

　 $\vec{a } = \left({}^t X X \right)^{-1} {}^t X \, \vec{y}$

で求められる．

ということです．

【補足】係数を求める公式の導出
　 $\{\varepsilon_{t}\}$ の2乗和は，

　 $\begin{eqnarray} {}^{t} \vec{\varepsilon}\, \vec{\varepsilon} &=& {}^{t} \! \left(\vec{y} - X \vec{a}\right) \left(\vec{y} - X \vec{a}\right) \\ &=& {}^{t} \vec{y} \, \vec{y} - {}^{t} \vec{y} \left(X \vec{a} \right) - {}^{t} \! \left( X \vec{a} \right) \vec{y} + {}^{t} \! \left( X \vec{a } \right) X \vec{a} \\ &=& {}^{t} \vec{y} \, \vec{y} - {}^{t} \vec{y}\, X \,\vec{a} - {}^{t} \vec{a}\, {}^{t}\! X \,\vec{y} + {}^{t} \vec{a}\, {}^{t} \! X X \, \vec{a} \\ &=& {}^{t} \vec{y} \, \vec{y} - {}^{t} \vec{y}\, X \,\vec{a} - {}^{t} \! \left({}^{t} \vec{y}\, X \,\vec{a} \right) + {}^{t} \vec{a}\, {}^{t} \! X X \, \vec{a} \end{eqnarray}$ 　

となります．さらに，この式の各項はすべてスカラーなので， ${}^{t} \vec{y}\, X \,\vec{a} = {}^{t} \! \left({}^{t} \vec{y}\, X \,\vec{a} \right)$ が成り立ちます．ですので，

　 $\begin{eqnarray} {}^{t} \vec{\varepsilon}\, \vec{\varepsilon} &=& {}^{t} \vec{y} \, \vec{y} - 2 {}^{t} \vec{y}\, X \,\vec{a} + {}^{t} \vec{a}\, {}^{t} \! X X \, \vec{a} \end{eqnarray}$ 　

　最小2乗法ということで，係数で偏微分したものがゼロになればよい．つまり，

　 $\displaystyle \frac{\partial \left( {}^{t} \vec{\varepsilon}\, \vec{\varepsilon}\right)}{\partial \vec{a}} = \vec{0}$

を仮定して計算します． ${}^{t} \vec{\varepsilon}\, \vec{\varepsilon} = {}^{t} \vec{y} \, \vec{y} - 2 {}^{t} \vec{y}\, X \,\vec{a} + {}^{t} \vec{a}\, {}^{t} \! X X \, \vec{a}$ の各項を微分してみると，

　 $\displaystyle \frac{\partial \left( {}^{t} \vec{y} \, \vec{y} \right)}{\partial \vec{a}} = \vec{0}$

　 $\displaystyle \frac{\partial \left( 2 {}^{t} \vec{y}\, X \,\vec{a} \right)}{\partial \vec{a}} = {}^{t} \!\! \left(2\, {}^{t} \vec{y}\, X \right)=2\, {}^{t} \! X \, \vec{y}$

　 $\displaystyle \frac{\partial \left( {}^{t} \vec{a}\, {}^{t} \! X X \, \vec{a} \right)}{\partial \vec{a}} =\left({}^{t} \! X X +\, {}^{t} \!\! \left({}^{t} \! X X \right) \right) \vec{a} =2\, {}^{t} \! X X \, \vec{a}$

となります．まとめると，

　 $\displaystyle \frac{\partial \left( {}^{t} \vec{\varepsilon}\, \vec{\varepsilon}\right)}{\partial \vec{a}} = - 2\, {}^{t} \! X \, \vec{y} + 2\, {}^{t} \! X X \, \vec{a} = \vec{0}$

になるので，後は， ${}^{t} \! X X \, \vec{a} = {}^{t} \! X \, \vec{y}$ を $\vec{a}$ について解くだけです．
【補足はここまで】

　Rスクリプトでの実装例は以下です．

# 長さNの時系列 y[1], y[2], ... , y[N] があるとき，
Y <- matrix(y[2:N],ncol=1)　# ncolで列数のみ指定
X <- matrix(c(rep(1,(N-1)),y[1:(N-1)]),ncol=2)　# ncolで列数のみ指定
#####################
# 係数の推定
A <- solve(t(X) %*% X) %*% (t(X) %*% Y)
# 残差の２乗平均
sig2 <- sum((Y - X %*% A)^2)/nrow(Y)

　Aに係数 $\{\hat{a}_0, \hat{a}_1\}$ が入っています．sig2には残差の２乗平均が入っています．

　Rのarコマンドを使って

ar(y, order.max=1, method="ols")

とすれば、同じ結果が得られますが，平均値が引かれた後の $\hat{a}_1$ のみが出てきます．

　今回は詳しく説明しませんが，自己回帰過程の次数は，Final Prediction Errorとか赤池情報量基準で決められます (元の時系列が自己回帰過程のサンプルであれば正しく決められます)．

2次自己回帰過程

　次は２次自己回帰過程

　 $y_t = a_0 + a_1 y_{t-1} + a_2 y_{t-2} +\varepsilon_t$

を仮定します．この場合，

　 $\begin{bmatrix} y_N \\ y_{N-1} \\ \vdots\\ y_{3} \\ \end{bmatrix} = \begin{bmatrix} 1 & y_{N-1} & y_{N-2} \\ 1 & y_{N-2} & y_{N-3} \\ \vdots & \vdots & \vdots \\ 1 & y_{2} & y_{1} \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ \end{bmatrix} + \begin{bmatrix} \varepsilon_{N} \\ \varepsilon_{N-1} \\ \vdots \\ \varepsilon_{3} \\ \end{bmatrix}$

となるので，これも， $\vec{y} = X \vec{a} + \vec{\varepsilon}$ の形になってます．Rスクリプトでは以下のようになります．

# 長さNの時系列 y[1], y[2], ... , y[N] があるとき，
Y <- matrix(y[3:N],　ncol=1) # ncolで列数のみ指定
X <- matrix(c(rep(1,(N-2)),y[2:(N-1)],y[1:(N-2)]),　ncol=3) # ncolで列数のみ指定
#####################
# 係数の推定
A <- solve(t(X) %*% X) %*% (t(X) %*% Y)
# 残差の２乗平均
sig2 <- sum((Y - X %*% A)^2)/nrow(Y)

多変量の場合 (完全な自己回帰過程ではありません)

自己ともう一つ別の変数に依存する過程

　多変量の場合も，計算の形・流れは単変量と変りません．ここでは，２変量 $y_t$ ， $z_t$ が観測可能で， $y_t$ が以下の式に従うことを仮定します．

　 $y_t = a_0 + a_1 y_{t-1} + b_1 z_{t-1} +\varepsilon_t$

　今回は， $z_t$ の過程を与えていないので，これは自己回帰とは呼べないと思います．

　このとき，長さNの時系列があれば，

　 $\begin{bmatrix} y_N \\ y_{N-1} \\ \vdots\\ y_{2} \\ \end{bmatrix} = \begin{bmatrix} 1 & y_{N-1} & z_{N-1} \\ 1 & y_{N-2} & z_{N-2} \\ \vdots & \vdots & \vdots \\ 1 & y_{1} & z_{1} \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ b_1 \\ \end{bmatrix} + \begin{bmatrix} \varepsilon_{N} \\ \varepsilon_{N-1} \\ \vdots \\ \varepsilon_{2} \\ \end{bmatrix}$

が準備できます．これも， $\vec{y} = X \vec{a} + \vec{\varepsilon}$ の形になってます．係数の推定値 $\{\hat{a}_0, \hat{a}_1, \hat{b}_1\}$ を求めるための式変形は，単変量のときと全く同じです．Rスクリプトは以下のようになります．

# 長さNの時系列 y[1], y[2], ... , y[N] があるとき，
Y <- matrix(y[2:N],　ncol=1) # ncolで列数のみ指定
X <- matrix(c(rep(1,(N-1)),y[1:(N-1)],z[1:(N-1)]),　ncol=3) # ncolで列数のみ指定
#####################
# 係数の推定
A <- solve(t(X) %*% X) %*% (t(X) %*% Y)
# 残差の２乗平均
sig2 <- sum((Y - X %*% A)^2)/nrow(Y)

瞬時的因果がある過程

　次は，瞬時的因果がある場合です．２変量 $y_t$ ， $z_t$ が観測可能で， $y_t$ が以下の式に従うことを仮定します (自己回帰してませんが)．

　 $y_t = a_0 + a_1 y_{t-1} + b_0 z_{t} + b_1 z_{t-1} +\varepsilon_t$

このとき，長さNの時系列があれば，

　 $\begin{bmatrix} y_N \\ y_{N-1} \\ \vdots\\ y_{2} \\ \end{bmatrix} = \begin{bmatrix} 1 & y_{N-1} & z_{N} & z_{N-1} \\ 1 & y_{N-2} & z_{N-1} & z_{N-2} \\ \vdots & \vdots & \vdots \\ 1 & y_{1} & z_{2} & z_{1} \\ \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ b_0 \\ b_1 \\ \end{bmatrix} + \begin{bmatrix} \varepsilon_{N} \\ \varepsilon_{N-1} \\ \vdots \\ \varepsilon_{2} \\ \end{bmatrix}$

が準備できます．これも， $\vec{y} = X \vec{a} + \vec{\varepsilon}$ の形になってます．やはり，係数を計算するための式変形は，他の場合と全く同じです．Rスクリプトは以下のようになります．

# 長さNの時系列 y[1], y[2], ... , y[N] があるとき，
Y <- matrix(y[2:N],　ncol=1) # ncolで列数のみ指定
X <- matrix(c(rep(1,(N-1)),y[1:(N-1)],z[2:N],z[1:(N-1)]),　ncol=4) # ncolで列数のみ指定
#####################
# 係数の推定
A <- solve(t(X) %*% X) %*% (t(X) %*% Y)
# 残差の２乗平均
sig2 <- sum((Y - X %*% A)^2)/nrow(Y)