はじめに
Rで自己回帰モデル（ARモデル）の導出を、ユールウォーカー法と最小二乗法で実施すると、下のように出力されます。ただし、次数は2で固定しています。
> # 人工データの作成（x_i=a_1*x_{i-1}+a2*x_{i-2}+epsilon_i,x_0=0,x_1=1,a_1=-1/4,a_2=1/8）
> x = c(0, 0)
> n = 100
> set.seed(0)
> for(i in 3:n){
+ x = c(x, - 1 / 4 * x[i-1]+ 1 / 8 * x[i-2] + rnorm(1, mean=0, sd=1/2))
+ }
>
> # ユールウォーカー法
> ar(x, aic = FALSE, order.max=2, method='yule-walker')
Call:
ar(x = x, aic = FALSE, order.max = 2, method = "yule-walker")
Coefficients:
1 2
-0.1398 0.1659
Order selected 2 sigma^2 estimated as 0.1941
>
> # 最小二乗法
> ar(x, aic = FALSE, order.max=2, method='ols')
Call:
ar(x = x, aic = FALSE, order.max = 2, method = "ols")
Coefficients:
1 2
-0.1432 0.1701
Intercept: 0.0002221 (0.04426)
Order selected 2 sigma^2 estimated as 0.1919
この結果の導出方法を、関数arを使わないRプログラムと数式で示します。
ユールウォーカー法
自己回帰モデルを数式で表します。
x_i = a_1 x_{i-1} + a_2 x_{i-2} + \epsilon_i,\ \ \ \epsilon_i\sim N(0,\sigma^2)
両辺に$x_{i-k}(k>0)$をかけます。
x_i x_{i-k} = a_1 x_{i-1}x_{i-k} + a_2 x_{i-2}x_{i-k} + \epsilon_ix_{i-k}
$E(x)=0$を仮定し両辺期待値をとると、
\gamma_k = a_1 \gamma_{k-1} + a_2 \gamma_{k-2}
となります。ただし、$\gamma_k$を$k$次の自己共分散とします。
$k=1,2$を代入すると、
\begin{align}
\gamma_1 &= a_1 \gamma_{0} + a_2 \gamma_{1} \\
\gamma_2 &= a_1 \gamma_{1} + a_2 \gamma_{0}
\end{align}
一方、自己回帰モデルの式の両辺に$x_i$をかけ、期待値をとると、
\gamma_0 = a_1 \gamma_1 + a_2 \gamma_2 + \sigma^2
これらより、
\begin{pmatrix}
a_1\\
a_2
\end{pmatrix}
=\begin{pmatrix}
\gamma_{0} & \gamma_{1} \\
\gamma_{1} & \gamma_{0}
\end{pmatrix} ^ {-1}
\begin{pmatrix}
\gamma_1\\
\gamma_2
\end{pmatrix}\\
\sigma^2 = \gamma_0 - a_1 \gamma_1 + a_2 \gamma_2
となります。
Rプログラムは下の通りです。
> # ユールウォーカー法
> gamma0 = sum((x - mean(x))^2) / n # 分散
> gamma1 = sum((x[2:n] - mean(x)) * (x[1:(n-1)] - mean(x))) / n # 1次の自己共分散
> gamma2 = sum((x[3:n] - mean(x)) * (x[1:(n-2)] - mean(x))) / n # 2次の自己共分散
> a = solve(matrix(c(gamma0, gamma1, gamma1, gamma0), nrow=2, ncol=2)) %*% c(gamma1, gamma2)
> a
[,1]
[1,] -0.139830
[2,] 0.165902
> (gamma0 - a[1] * gamma1 - a[2] * gamma2) / (n - 3) * n
[1] 0.19411
この結果はar(method='yule-walker')関数の結果と一致します。
最小二乗法
\boldsymbol{y} = \begin{pmatrix}
x_3\\
x_4\\
\vdots\\
x_n
\end{pmatrix}
,\ X=
\begin{pmatrix}
1 & x_2 & x_1\\
1 & x_3 & x_2\\
\vdots & \vdots & \vdots\\
1 & x_{n-1} & x_{n-2}
\end{pmatrix}
,\ \boldsymbol{a}=
\begin{pmatrix}
a_0\\
a_1\\
a_2
\end{pmatrix}
とします。
\boldsymbol{\epsilon} = \boldsymbol{y} - X \boldsymbol{a}
とするとき、$||\boldsymbol{\epsilon}||$を最小とする$\boldsymbol{a}$は、
\boldsymbol{a} = (X^TX)^{-1}X^T\boldsymbol{y}
で与えられます。
$\boldsymbol{\epsilon}\sim N(0,\sigma^2I)$としたとき、$\sigma^2$は下式で推定されます。
\hat{\sigma}^2 = \frac{||\boldsymbol{y}-X\boldsymbol{a}||}{n-2}
$V(a_0)$は下式で推定されます。
V(a_0)=\hat\sigma^2(X^TX)^{-1}
Rプログラムは下の通りです。
> # 最小二乗法
> X = cbind(1, x[2:(n-1)]-mean(x),x[1:(n-2)]- mean(x))
> y = x[3:n] - mean(x)
>
> a = solve(t(X) %*% X) %*% t(X) %*% y
> a
[,1]
[1,] 0.0002221244
[2,] -0.1431719722
[3,] 0.1700645819
> sigma2 = sum((y - X %*% a)^2) / (n - 2)
> sigma2
[1] 0.1918545
> sqrt(solve(t(X) %*% X)[1, 1] * sigma2)
[1] 0.0442592
この結果はar(method='ols')関数の結果と一致します。