餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

行列を使った計算で、SURの感じを掴んでみよう

計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ないSUR *1 をふと思い出したので、手順を確認してみました。

1. 2つの方程式からなるSURモデル

理屈はGreeneのEconometric Analysisに詳しく載っているので参照して欲しいのですが、SURの概要を説明します。
以下の連立方程式を同時推定することを考えます*2

{\displaystyle 
\begin{eqnarray}
  \left\{
    \begin{array}{l}
		y_{1} = X_{1}\beta_1 + \epsilon_{1} \\
		y_{2} = X_{2}\beta_2 + \epsilon_{2}
    \end{array}
  \right.
\end{eqnarray}
}

 yは被説明変数のベクトル、 Xは説明変数の行列、 \betaは係数、 \epsilonは誤差項のベクトルです。添字の数字は方程式を表します。
一つ一つ推定しても良い気がしますが、 \mbox{COV}(\epsilon_1, \epsilon_2) \neq 0の場合は推定量を改善する余地があります。

例えば、天候データが無いときに2つの品種の収穫高と肥料の関係を推定するとして、日照量や気温など2つの方程式に同時に同じような影響を与える要因がある場合、 \mbox{COV}(\epsilon_1, \epsilon_2)という情報も使う方が推定結果は真の値に近づきます。

さて、最初は \mbox{COV}(\epsilon_1, \epsilon_2)の情報が無いので、以下をOLSで推定します。

{\displaystyle 
	\begin{pmatrix}
		y_{1} \\
		y_{2}
	\end{pmatrix}
	=
	\begin{pmatrix}
		X_{1} & 0 \\
		0 & X_{2}
	\end{pmatrix}
	\begin{pmatrix}
		\beta_{1} \\
		\beta_{2}
	\end{pmatrix}
	+
	\begin{pmatrix}
		\epsilon_1 \\
		\epsilon_2
	\end{pmatrix}
}

これで、以下の分散共分散行列 \Sigmaがつくれるようになりました。

{\displaystyle 
	\Sigma= 
	\begin{pmatrix}
		\mbox{COV}(\epsilon_1, \epsilon_1) & \mbox{COV}(\epsilon_1, \epsilon_2) \\
		\mbox{COV}(\epsilon_1, \epsilon_2) & \mbox{COV}(\epsilon_2, \epsilon_2)
	\end{pmatrix}
}

これをウェイトに使って、不均一分散の補正をした回帰(FGLS)をかけます。
つまり、以下のように表記を簡素化して、

{\displaystyle 
	y = X\beta + \epsilon
}

以下のようにGLS推定量を求めます。


\newcommand{\argmin}{\mathop{\rm arg~min}\limits}
{\displaystyle 
	\hat{\beta} = \argmin_{\beta} \space  (y - X \beta)^t ( \Sigma^{-1} \otimes I_n ) (y - X \beta)
}

ここで、 \otimesクロネッカー積、 I_nはn行n列の対角行列、 \Sigma \otimes I_nは共分散を対角成分にとる細胞で構成される行列です。 nは観測数です。後述するコマンドを見れば、何が起きるかは一目瞭然だと思います。

この式を満たすために、 \beta微分して一階条件を出して整理すると、

{\displaystyle 
	\hat{\beta} = (X^t ( \Sigma^{-1} \otimes I_n ) X)^{-1} X^t ( \Sigma^{-1} \otimes I_n ) y
}

となります。 (\Sigma \otimes I_n)^{-1} = \Sigma^{-1} \otimes I_nに注意してください。

2. RによるSUR推定

systemfitパッケージを使えば良いのですが、前節で説明したモデルとの対応関係が分からなくなるので、行列を使って計算します。

#
# データセットの作成
#
set.seed(20181103)
n <- 50 # 標本サイズが小さい方が、OLSとSURの差が出る

l <- rnorm(n, sd=2) # 共分散を持つようにする
e1 <- rnorm(n, sd=1) + l
e2 <- rnorm(n, sd=1) + l

x1 <- runif(n, min=0, max=2)
x2 <- runif(n, min=-2, max=2)
y1 <- 1 + 2*x1 + e1
y2 <- 3 + 2*x2 + e2

X1 <- matrix(c(rep(1, n), x1, rep(0, 2*n)), n, 4)
X2 <- matrix(c(rep(0, 2*n), rep(1, n), x2), n, 4)
X <- rbind(X1, X2)
y <- c(y1, y2)

#
# 以下、実際のSUR推定量の計算
#
# 1段階目のOLS推定量を計算
beta_ols <- solve(t(X)%*%X)%*%t(X)%*%y

# 誤差項を計算
r_e <- y - X %*% beta_ols
r_e1 <- r_e[1:n]
r_e2 <- r_e[(n+1):(2*n)]

# ∑を計算
sigma <- matrix(c(cov(r_e1,r_e1), cov(r_e1,r_e2), cov(r_e1,r_e2), cov(r_e2,r_e2)), 2, 2)

# SUR推定量を計算
beta_sur <- solve(t(X) %*% (solve(sigma) %x% diag(n)) %*% X) %*% t(X) %*% (solve(sigma) %x% diag(n)) %*% y

なお、beta_ols - beta_surで推定量にどれぐらい違いがあるか分かるので、nの値を変えて試してみてください。t値やP値はOLSと同様に計算できます。また、sigma %x% diag(n)とすると、 n  \times nの巨大行列ですが、 \Sigma \otimes I_nがどうなっているか分かります。

*1: すぐに同時性ガガガ…と言われるので、GMMか3SLS、もしくは全情報最尤法を使うのが一般的で、SUR自体は3SLSを学ぶ前の予備知識的な扱いが多いと思います。

*2:もっとも単純なモデルを考えましたが、二つの式で、係数が同じ値を持つような制約を置くこともできます。