Rで計量経済学の変量効果モデルを推定

だいぶ前にパッケージを使わず固定効果モデルで推定する例は書いていたのですが、変量効果モデルで推定する例も書いてみました。計量経済学で言う変量効果モデルは、観測値ごとの誤差 $u_{it}$ に加えて、個体ごとの誤差 $c_{i}$ がある線形回帰モデル $y_{it} = X\beta + c_{i} + u_{it}$ です。 $y$ は被説明変数、 $X$ は説明変数、 $\beta$ は係数、 $i$ は個体を示す添字、 $t$ は時点を示す添字になります。実用ではパッケージを使うべきですが、説明でもパッケージを使うとブラックボックス感が出てしまうので。

1.データセット

定番テキストのWooldridgeではjtrainデータセット *1が例になっていたりするのですが、実データだと真のDGPとそのパラメーターが分からないので、乱数から生成します。

mkdata <- function(n_of_t, n_of_i){
    # n_of_t: 時点の数, n_of_i: 個体の数

    # 個体番号をつくる
    id <- rep(1:n_of_i, each = n_of_t, min=-0.5, max=0.5)

    # 時点変数をつくる
    t <- rep(1:n_of_t, n_of_i)

    # 説明変数x,zをつくる
    x <- runif(n_of_i*n_of_t)
    z <- runif(n_of_i*n_of_t)

    # 誤差項をつくる
    e_ind <- rep(rnorm(n_of_i, sd = 2), each = n_of_t)
    e_ids <- rnorm(n_of_i*n_of_t)
    e <- e_ind + e_ids

    # 従属変数yをつくる
    y <- 1 + 2*x - 3*z + e

    data.frame(id = as.factor(id), t = as.factor(t), y, x, z)
}

set.seed(1303)
df01 <- mkdata(5, 100)

2. 行列演算などで推定

定番テキストのWooldridgeのGSLとRandom Effect Modelの説明に沿って計算します。OLSの推定結果から、個体の分散 $\sigma^2_c$ と観測ごとの分散 $\sigma^2_u$ を推定して、ウェイト行列 $\Sigma$ をつくってGLSにかけます。

$\Sigma = \begin{bmatrix} \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c & \sigma^2_c \\ \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_c & \sigma^2_u + \sigma^2_c \\ \end{bmatrix}$

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

自由度がややこしいので注意してください。

frml <- y ~ x + z
X <- model.matrix(frml, data = df01)
mf <- model.frame(frml, data = df01)
y <- model.response(mf)
beta_ols <- solve(t(X)%*%X) %*% t(X) %*% y
residuals <- X %*% beta_ols - y
df <- nrow(X) - ncol(X)
sigma2_ols <- sum(residuals^2)/df
T <- length(levels(df01$t))
n <- length(levels(df01$id))

# total error variance
sigma2_v <- sum(residuals^2)/df

# individual error variance
sigma2_c <- with(df01, {
    sum(sapply(levels(id), function(i){
        r <- residuals[id == i]
        cp <- outer(r, r)
        sum(cp[upper.tri(cp, diag = FALSE)])
        # 以下と同じ結果
        # sum <- 0
        # for(t in 1:(T-1)){
        #     for(k in (t+1):T){
        #         sum <- sum + r[t]*r[k]
        #     }
        # }
        # sum
    }))/((n*T*(T-1)/2) - ncol(X))
})

# idiosyncratic error variance
sigma2_u <- sigma2_v - sigma2_c

omega <- matrix(sigma2_c, T, T) + sigma2_u*diag(T)
iomega <- diag(n) %x% solve(omega)
beta_gls <- solve(t(X) %*% iomega %*% X) %*% t(X) %*% iomega %*% y

# Usual Variance Matrix Estimator
usual_vcov <- solve(t(X) %*% iomega %*% X)

# Robust Variance Matrix Estimator
rbst_vcov <- with(df01, {
    dfs <- split(df01, id)
    so <- solve(omega)
    A <- matrix(0, ncol(X), ncol(X))
    B <- matrix(0, ncol(X), ncol(X))
    n <- length(dfs)
    for(i in 1:n){
        X <- model.matrix(frml, data = dfs[[i]])
        y <- dfs[[i]]$y
        u <- y - X %*% beta_gls
        A <- A + t(X) %*% so %*% X
        B <- B + t(X) %*% so %*% u %*% t(u) %*% so %*% X
    }
    solve(A) %*% B %*% solve(A)
})

3. パッケージで推定

比較のために、plmパッケージで推定します。

library(plm)
pd01 <- pdata.frame(df01, index=c("id"))
r_plm <- plm(frml, model = "random", data = pd01)

4. 推定結果を比較

係数の推定量と標準誤差を比較してみましょう。

beta <- matrix(
    c(1, 2, -3, beta_ols, beta_gls, coef(summary(r_plm))[, 1]),
    ncol(X), dimnames = list(rownames(beta_gls), c("true", "OLS", "GLS", "plm/random")))

sigma <- matrix(
    c(sqrt(sigma2_ols*diag(solve(t(X) %*% X))), sqrt(diag(usual_vcov)), sqrt(diag(rbst_vcov)), coef(summary(r_plm))[, 2]),
    ncol(X), dimnames = list(rownames(beta_gls), c("OLS", "GLS/Usual", "GLS/Robust", "plm/random")))

係数が

	true	OLS	FGLS	plm/random
(Intercept)	1.0000	0.7393	1.1459	1.1471
x	2.0000	2.6266	2.0862	2.0845
z	-3.0000	-2.7875	-3.0411	-3.0419

標準誤差が

	OLS	FGLS/Usual	FGLS/Robust	plm/random
(Intercept)	0.28772	0.25141	0.25710	0.25087
x	0.36308	0.18589	0.20435	0.18237
z	0.35117	0.17542	0.15863	0.17208

となります。FGLSで推定した係数は、OLSよりも真の値に近くなっており、係数も標準誤差もplmパッケージとほとんど同じです。

ほとんど同じと言うところに引っかかると思いますが、plmパッケージでは以下の連立方程式を解いて個体内変化と個体間相違を計算して、idiosyncratic error varianceとindividual error varianceにしており、分散分解の方法が教科書のものと異なっているためです。