dとsの式でpが内生変数、zが操作変数、iが個体を表す番号ですが、こんなんで。
noi <- 20 # 個体数 t <- 5 # 観測期間 obs <- t*noi i <- rep(1:(noi), each=t) fe <- runif(noi, min=0, max=100) a0 <- rep(fe, each=t) a1 <- -1 b0 <- 2 b1 <- 3 b2 <- 4 z <- runif(obs, min=0, max=3) mu <- rnorm(obs, mean=0, sd=2) nu <- rnorm(obs, mean=0, sd=1) p <- (b0 - a0 + b2*z + nu - mu)/(a1 - b1) d <- a0 + a1*p + mu s <- b0 + b1*p + b2*z + nu
within変換をかけてIVで推定
within_transfer <- function(x, i){ m <- tapply(x, i, mean) x - rep(m, each=t) } w_d <- within_transfer(d, i) w_z <- within_transfer(z, i) w_p <- within_transfer(p, i) zm <- matrix(c(c(w_z)), obs, 1) xm <- matrix(c(c(w_p)), obs ,1) estimated_a1 <- solve(t(zm) %*% xm) %*% (t(zm) %*% w_d)
観測数obsが2000ぐらい無いと誤差多し
Stata風に切片項を計算する
estimated_mu <- w_d - estimated_a1 %*% w_p ma_d <- mean(d) mi_d <- rep(tapply(d, i, mean), each=t) ma_z <- mean(z) mi_z <- rep(tapply(z, i, mean), each=t) ma_p <- mean(p) mi_p <- rep(tapply(p, i, mean), each=t) estimated_a0 <- (d - mi_d + ma_d) - (p - mi_p + ma_p)*c(estimated_a1) - c(estimated_mu)
切片項のP値が欲しい場合
StataのFAQを見ると、within変換をかけた変数に全体の平均値を加算してから推定をすることで、切片項の有意性を出していました。
w_d <- within_transfer(d, i) + mean(d) w_z <- within_transfer(z, i) + mean(z) w_p <- within_transfer(p, i) + mean(p) zm <- matrix(c(rep(1, obs), c(w_z)), obs, 2) xm <- matrix(c(rep(1, obs), c(w_p)), obs ,2) estimated_a <- solve(t(zm) %*% xm) %*% (t(zm) %*% w_d)
分散共分散行列と標準誤差の計算は以下です。自由度の計算でパネルの数（個体種類）も加味しないと、within推定とdummy variable estimatorの標準誤差が一致しないことに注意しましょう。
df <- (obs - 2) - (noi - 1) ssr <- sum((w_d - xm %*% estimated_a)^2) s2 <- ssr/df vcov <- s2*solve( t(xm) %*% zm %*% solve(t(zm) %*% zm) %*% t(zm) %*% xm ) se <- sqrt(diag(vcov))
one-wayクラスター頑強標準誤差の計算
気づくと一般化したone-wayのロバスト標準誤差（ここでは個体ごとの不均一分散を調整）の場合は以下です。「切片項のP値が欲しい場合」の続きとして実行できます。
# 時点を示す番号（two-way clustering用） # ti <- rep(1:t, noi) # 残差からウェイトΩを計算 residuals <- (w_d - xm %*% estimated_a) omega <- matrix(0, obs, obs) # Ωは対角成分としてt×tの細胞をnoi個とる for(j in 1:obs){ omega[,j] <- residuals[j]*residuals*(i[j]==i) # two-way clusteringのときは以下にする # omega[,j] <- residuals[j]*residuals*(i[j]==i|ti[j]==ti) } # plmパッケージのvcovHC(model, type = 'HC0') と同じ値 vcov_hc0 <- solve( t(xm) %*% zm ) %*% (t(zm) %*% omega %*% zm) %*% solve( t(zm) %*% xm ) # 自由度を調整 df <- nrow(zm) - ncol(zm) dfcw <- df / (df - (noi - 1)) dfc <- (noi / (noi - 1))*((obs - 1)/df) vcov <- dfc*vcov_hc0*dfcw # 標準偏差を計算 se <- sqrt(diag(vcov))
IVとFGLSのあわせ技になっていますが、xmとzmを同じにすればIVでないのと同じになります。
なお、太田 (2013)と、パッケージを使わないようにコードは変えましたが、ストックホルム大学のMahmood Arai教授のレクチャーノートを参考にしました。
plmパッケージでのone-wayクラスター頑強標準誤差の計算の仕方
もっとも需要が大きそうな話を書き忘れていたので4年経過していますが追記します。Stata風に切片項を推定する方法は無さそうですが、自由度調整をすれば同様にクラスター頑強標準誤差になります。
library(plm) # データフレームをつくる df01 <- data.frame(i, t = rep(1:5, noi), d, s, p, z) # パネルデータ分析用データフレームにする pdf01 <- pdata.frame(df01, index=c("i", "t")) # within推定をする r_plm <- plm(d ~ p | z, data = pdf01, model = "within") # 自由度を調整 # df <- with(r_plm, nrow(model)- (length(coefficients)) + 1) # dfcw <- df / (df - (noi - 1)) # noi:観測数 # dfc <- (noi / (noi - 1))*((obs - 1)/df) # obs: サンプルサイズ # vcov <- dfc*vcovHC(r_plm, type = 'HC0')*dfcw # # 補正VCOVで標準誤差を計算 # summary(r_plm, vcov = vcov) # summary(r_plm, vcov = vcovHC(r_plm, type = 'sss')) # これでStata互換になるそうです：https://blog.theleapjournal.org/2016/06/sophisticated-clustered-standard-errors.html
