はじめに

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')関数の結果と一致します。