行列を使った計算で、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)]

# ∑を計算
df <- n - 2 # 自由度は複数の方程式の最も小さい値になる
sigma <- matrix(c(sum(r_e1^2)/df, sum(r_e1*r_e2)/df, sum(r_e1*r_e2)/df, sum(r_e2^2)/df), 2, 2)

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

# 分散共分散行列を計算
r <- (y - X %*% beta_sur)
vcov <- solve(t(X) %*% (solve(sigma) %x% diag(n)) %*% X) # Usual Variance Matrix (Wooldridge (2002) p.161)

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