はじめに

簡単な単回帰の例で外生性（exogeneity）と内生性（endogeneity）の違いを実感する．

\[Y=\alpha+\beta X + \varepsilon\] \(X\)が外生的であるとは \[E[\varepsilon|X]=0 \Longrightarrow E[\varepsilon]=0, E[\varepsilon X]=0 \Longrightarrow Cov(\varepsilon,X)=0\Longrightarrow cor(\varepsilon,X)=0.\] 内生的であるとは\(E[\varepsilon|X]\neq 0\)．

外生的な例

相関のない二次元正規分布から\(X\)と\(\varepsilon\)のサンプリング. 真のモデルとして以下を想定 \[Y=0+1* X + \varepsilon\]

data<-as.data.frame(mvrnorm(100, c(10,0), matrix(c(10, 0, 0, 1), 2, 2)))
colnames(data)<-c("x","e")
data$y<-0+1*data$x+data$e
plot(data)

\(X_i\)と\(\varepsilon_i\)の内積をとってみる

data$x %*% data$e

##           [,1]
## [1,] -150.2438

回帰係数は

lm(y ~ x ,data)

## 
## Call:
## lm(formula = y ~ x, data = data)
## 
## Coefficients:
## (Intercept)            x  
##    -0.02147      0.98681

サンプリングでゆらぐので，1000回繰り返して標本分布を見てみる．

ip_ex<-c()
beta_ex<-c()
for (i in 1:1000) {
  ## random sampling from non correlated multi normal dist.
  data<-as.data.frame(mvrnorm(100, c(10,0), matrix(c(1, 0, 0, 1), 2, 2)))
  colnames(data)<-c("x","e")
  data$y<-0+1*data$x+data$e

  ## inner product
  data$x %*% data$e -> ip_ex[i]
  
  ## beta
  unname(lm(y ~ x ,data)$coef[2]) -> beta_ex[i]
  
}

内積は

hist(ip_ex,col="skyblue")

mean(ip_ex)

## [1] 0.3522053

0に近い？

回帰係数は

hist(beta_ex,col="skyblue")

mean(beta_ex)

## [1] 1.000352

平均的に理論値の1．

内生的な例

相関0.7の二次元正規分布から\(X\)と\(\varepsilon\)のサンプリング. 真のモデルとして以下を想定 \[Y=0+1* X + \varepsilon\]

data2<-as.data.frame(mvrnorm(100, c(10,0),
                             matrix(c(10, 0.7*sqrt(10)*1,  0.7*sqrt(10)*1, 1), 2, 2)))
colnames(data2)<-c("x","e")
data2$y<-0+1*data2$x+data2$e
plot(data2)

\(X_i\)と\(\varepsilon_i\)の内積をとってみる

data2$x %*% data2$e

##          [,1]
## [1,] 345.8846

回帰係数は

lm(y ~ x ,data2)

## 
## Call:
## lm(formula = y ~ x, data = data2)
## 
## Coefficients:
## (Intercept)            x  
##      -1.883        1.198

サンプリングでゆらぐので，1000回繰り返して標本分布を見てみる．

ip_end<-c()
beta_end<-c()
for (i in 1:1000) {
  ## random sampling from non correlated multi normal dist.
  data<-as.data.frame(mvrnorm(100, c(10,0),
                              matrix(c(10, 0.7*sqrt(10)*1,  0.7*sqrt(10)*1, 1), 2, 2)))
  colnames(data)<-c("x","e")
  data$y<-0+1*data$x+data$e
  
  ## inner product
  data$x %*% data$e -> ip_end[i]
  
  ## beta
  unname(lm(y ~ x ,data)$coef[2]) -> beta_end[i]
  
}

内積は

hist(ip_end,col="skyblue")

mean(ip_end)

## [1] 229.8664

0からはずれる．

回帰係数は

hist(beta_end,col="skyblue")

mean(beta_end)

## [1] 1.222172

平均的に理論値の1から外れる．

操作変数（IV）

\(X, Y, \varepsilon\)のほかに，\(X\)と相関する外生変数\(Z\)があったら，これを操作変数（instrumental variable）として推定してやればいいだろう． \[E[\varepsilon | Z]=0, E[X|Z]\neq0 \Longrightarrow Cov(\varepsilon, Z)=0, Cov(X,Z)\neq 0\] このとき，操作変数法による\(X\)の係数\(\beta\)の推定量は， \[\hat{\beta}_{\mathrm{IV}}=\frac{\sum_i(Y_i-\bar{Y})(Z_i-\bar{Z})}{\sum_i(X_i-\bar{X})(Z_i-\bar{Z})}\]

以下，\(cor(X,Z) = 0.7, cor(X,\varepsilon) =0.7, cor(Z,\varepsilon) =0\)を仮定した同時正規分布からサンプリング． 真のモデルとして以下を想定 \[Y=0+1* X + \varepsilon,\quad X=\gamma_0+\gamma_1Z\] ただし，\(\gamma_0=E[X]-\gamma_1E[Z], \gamma_1=0.7\sqrt{V[X]/V[Z]}\)．

data3<-as.data.frame(mvrnorm(100, c(10,5,0),
                             matrix(c(10, 0.7*sqrt(10)*sqrt(5),  0.7*sqrt(10)*1,
                                      0.7*sqrt(10)*sqrt(5),5,0,
                                      0.7*sqrt(10)*1,0,1), 3, 3)))
colnames(data3)<-c("x","z","e")
data3$y<-0+1*data3$x+data3$e
plot(data3)

まず，ふつうにOLSをしてみる．

cov(data3$y,data3$x)/cov(data3$x,data3$x)

## [1] 1.199835

lm(y ~ x ,data3)

## 
## Call:
## lm(formula = y ~ x, data = data3)
## 
## Coefficients:
## (Intercept)            x  
##      -2.024        1.200

1000回繰り返して推定値の標本分布を見てみる．

beta_nonIV<-c()
for (i in 1:1000) {
  ## random sampling from non correlated multi normal dist.
  data<-as.data.frame(mvrnorm(100, c(10,5,0),
                              matrix(c(10, 0.7*sqrt(10)*sqrt(5),  0.7*sqrt(10)*1,
                                       0.7*sqrt(10)*sqrt(5),5,0,
                                       0.7*sqrt(10)*1,0,1), 3, 3)))
  colnames(data)<-c("x","z","e")
  data$y<-0+1*data$x+data$e
  
  ## beta
  unname(lm(y ~ x ,data)$coef[2]) -> beta_nonIV[i]
  
}
hist(beta_nonIV,col="skyblue")

mean(beta_nonIV)

## [1] 1.22103

\(X\)の効果が過大評価されている．

つぎに，操作変数を入れて\(X\)の係数を推定する．AERパッケージのivregを使う（点推定値だけなら，共分散の計算で求まる）．

cov(data3$y,data3$z)/cov(data3$x,data3$z)

## [1] 1.02678

ivreg(y ~ x | z, data = data3)

## 
## Call:
## ivreg(formula = y ~ x | z, data = data3)
## 
## Coefficients:
## (Intercept)            x  
##      -0.164        1.027

いいじゃないかいいじゃないか． 1000回繰り返して推定値の標本分布を見てみる．

beta_IV<-c()
for (i in 1:1000) {
  ## random sampling from non correlated multi normal dist.
  data<-as.data.frame(mvrnorm(100, c(10,5,0),
                              matrix(c(10, 0.7*sqrt(10)*sqrt(5),  0.7*sqrt(10)*1,
                                       0.7*sqrt(10)*sqrt(5),5,0,
                                       0.7*sqrt(10)*1,0,1), 3, 3)))
  colnames(data)<-c("x","z","e")
  data$y<-0+1*data$x+data$e
  
  ## beta
  unname(ivreg(y ~ x | z, data = data)$coef[2]) -> beta_IV[i]
  
}
hist(beta_IV,col="skyblue")

mean(beta_IV)

## [1] 0.9946882

こういうのでいいんだよこういうので．