餡子付゛録゛

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

FE-IV練習用データ生成

dとsの式でpが内生変数、zが操作変数ですが、こんなんで。

noi <- 20 # 個体数
t <- 5 # 観測期間
obs <- t*noi
i <- c(t(replicate(t, 1:(noi))))
fe <- runif(noi, min=0, max=100)
a0 <- c(t(replicate(t, fe)))
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 - c(t(replicate(t, m)))
}

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 <- c(t(replicate(t,tapply(d, i, mean))))

ma_z <- mean(z)
mi_z <- c(t(replicate(t,tapply(z, i, mean))))

ma_p <- mean(p)
mi_p <- c(t(replicate(t,tapply(p, i, mean))))

estimated_a0 <- (d - mi_d + ma_d) - (p - mi_p + ma_p)*c(estimated_a1) - c(estimated_mu)