計量の教科書には載っているものの非線形モデルには使えず、内生性を制御できないので顧みられることの少ないSUR *1 をふと思い出したので、手順を確認してみました。
1. 2つの方程式からなるSURモデル
理屈はGreeneのEconometric Analysisに詳しく載っているので参照して欲しいのですが、SURの概要を説明します。
以下の連立方程式を同時推定することを考えます*2。
は被説明変数のベクトル、は説明変数の行列、は係数、は誤差項のベクトルです。添字の数字は方程式を表します。
一つ一つ推定しても良い気がしますが、の場合は推定量を改善する余地があります。
例えば、天候データが無いときに2つの品種の収穫高と肥料の関係を推定するとして、日照量や気温など2つの方程式に同時に同じような影響を与える要因がある場合、という情報も使う方が推定結果は真の値に近づきます。
さて、最初はの情報が無いので、以下をOLSで推定します。
これで、以下の分散共分散行列がつくれるようになりました。
これをウェイトに使って、不均一分散の補正をした回帰(FGLS)をかけます。
つまり、以下のように表記を簡素化して、
以下のようにGLS推定量を求めます。
ここで、はクロネッカー積、はn行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)とすると、の巨大行列ですが、がどうなっているか分かります。