餡子付゛録゛

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

FE-IV練習用データ生成

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)