問題設定

賃金($w$万円/年)と学歴(教育年数:$X$)との関係を考える。

学歴が1年伸びると賃金が25万円/年増加するとしよう。ただし、総賃金は本人の能力$Q\in[-5,5]$と性別$S\in{-1,1}$によって就ける業種や職種が違うため、賃金は大きく異なる。すなわち、$Q=-5$に比べて$Q=5$の人は平均で200万円ほど高い仕事に就けるし、男性($S=1$)の方が女性($S=-1$)よりも平均で200万円高い賃金を得ている。この関係式は次のようになるものとしよう。

$$w=300+25X+20Q+100S+\varepsilon$$

ただし、$\varepsilon\sim N(0,20)$。

ところが、ここで注意すべきは、学歴$X$の選択が内生的に決まるということ。すなわち、$Q$や$S$によって期待できる賃金$w$が違うし、学歴$X$年にするには相応の教育コストがかかる。したがって、$Q$ の高い人ほど長い学歴を選択するだろうし、男性ほど学歴が長くなる可能性が高い。もちろん、学歴は家庭環境($K\sim N(0,2)$)にもよるから、その関係は、

$$X=10+0.5Q+5S+K+\eta$$

と表せる。ただし、$\eta\sim N(0,1)$。

サンプルデータの生成

 N<-1000
 Q<-runif(N,-5,5)
 S<-sample(c(-1,1),N,T)
 K<-rnorm(N,0,2)
 eta<-rnorm(N,0,1)
 X<-15+0.5*Q+2.5*S+K+eta
 eps<-rnorm(N,0,5)
 w<-300+25*X+20*Q+100*S+eps

OLS

内生性があろうとなかろうと、$plim({^tX}\varepsilon)=0$でありさえすれば、OLSで一致性は得られる。

> o01<-lm(w~X+Q+S)
> summary(o01)

Call:
lm(formula = w ~ X + Q + S)

Residuals:
Min       1Q   Median       3Q      Max
-14.4623  -3.3285   0.2185   3.2389  17.3528

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 300.65518    1.04473   287.8   <2e-16 ***
X            24.96559    0.06871   363.3   <2e-16 ***
Q            20.14302    0.06390   315.2   <2e-16 ***
S           100.23460    0.23223   431.6   <2e-16 ***

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.933 on 996 degrees of freedom
Multiple R-squared: 0.9993,     Adjusted R-squared: 0.9993
F-statistic: 5.088e+05 on 3 and 996 DF,  p-value: < 2.2e-16

o01<-lm(w~X+Q+S)
 summary(o01)

問題なのは、$Q$や$S$が抜け落ちた場合。$Q$や$S$が抜けても、それらが$X$と相関していなければ、誤差項が$20Q+200S+\varepsilon$となるだけで、$plim{^tX}(20Q+200S+\varepsilon)=0$だから、一致性は得られる。ところが、$X$が$Q$や$S$と相関していると、これが欠落することで、$plim^tX(20Q+200S+\varepsilon)\neq 0$となり、一致性は得られない。

> o01<-lm(w~X)
>  summary(o01)
Call:
lm(formula = w ~ X)

Residuals:
Min       1Q   Median       3Q      Max
-221.347  -49.288    2.695   51.354  220.723

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -72.6660     9.6440  -7.535 1.09e-13 ***
X            49.5551     0.6324  78.364  < 2e-16 ***

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 72.15 on 998 degrees of freedom
Multiple R-squared: 0.8602,     Adjusted R-squared: 0.8601
F-statistic:  6141 on 1 and 998 DF,  p-value: < 2.2e-16

o01<-lm(w~X)
summary(o01)

操作変数法

$X$を決める変数の中で、家庭環境$K$だけは賃金$w$に影響しないから、これを操作変数に用いる。

 > o02<-lm(X~K)
> Xh<-predict(o02)
> X1<-cbind(matrix(1,N,1),X)
> Xh1<-cbind(matrix(1,N,1),Xh) > cov(Xh1,0.5*Q+5*S+eta)
[,1]
0.00000000
Xh 0.03686849
> solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
[,1]
301.79073
X 24.28418

o02<-lm(X~K)
Xh<-predict(o02)
X1<-cbind(matrix(1,N,1),X)
Xh1<-cbind(matrix(1,N,1),Xh)
cov(Xh1,0.5*Q+5*S+eta)
solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w

実際、$w$を分解してみると、

>   solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
       [,1]
301.79073
X  24.28418
> solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(300+25*X)
  [,1]
300
X   25
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(20*Q)
       [,1]
16.106584
X -1.432243
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(100*S)
         [,1]
-15.0677828
X   0.7604274
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%eps
         [,1]
0.75192925
X -0.04400716

 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(300+25*X)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(20*Q)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(100*S)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%eps

となるので、一致性が得られそうである。ちなみに、後の考察のために、次の値も確認しておく。

> solve(t(Xh1)%*%X1)
                          Xh
0.055824241 -0.0036999162
X -0.003699916  0.0002496958
>  t(Xh1)%*%(20*Q)
         [,1]
-5115.967
Xh -81542.804
>  t(Xh1)%*%(100*S)
        [,1]
-3800.00
Xh -53261.84
>  t(Xh1)%*%eps
        [,1]
99.8445
Xh 1303.2224

 solve(t(Xh1)%*%X1)
 t(Xh1)%*%(20*Q)
 t(Xh1)%*%(100*S)
 t(Xh1)%*%eps

弱相関の問題

ただし、操作変数$K$と$X$との相関が弱かった場合、一致性が得られなくなる。$X$の説明変数から$K$が脱落している(係数が0)の場合を考えよう。

> X<-10+0.5*Q+5*S+0.1*K+eta
>  w<-300+25*X+20*Q+100*S+eps
>  o02<-lm(X~K)
>  Xh<-predict(o02)
>
>  X1<-cbind(matrix(1,N,1),X)
>  Xh1<-cbind(matrix(1,N,1),Xh)
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
       [,1]
354.12528
X  18.50304

 X<-10+0.5*Q+5*S+0.1*K+eta
 w<-300+25*X+20*Q+100*S+eps
 o02<-lm(X~K)
 Xh<-predict(o02)
 
 X1<-cbind(matrix(1,N,1),X)
 Xh1<-cbind(matrix(1,N,1),Xh)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w

バイアスが生じている。$w$を分解してみると、

> solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(300+25*X)
  [,1]
300
X   25
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(20*Q)
       [,1]
120.81933
X -12.99935
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(100*S)
        [,1]
-70.663384
X   6.901801
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%eps
        [,1]
3.9693363
X -0.3994184

 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(300+25*X)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(20*Q)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(100*S)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%eps

$Q$や$S$の項のバイアスが大きい。$\varepsilon$の項のバイアスもそこそこ生じている。

さらに、各式の前段と後段を分けてみると、

> solve(t(Xh1)%*%X1)
                      Xh
1.9315115 -0.19927210
X -0.1992721  0.02056935
>  t(Xh1)%*%(20*Q)
         [,1]
-5115.967
Xh -50194.523
>  t(Xh1)%*%(100*S)
        [,1]
-3800.00
Xh -36478.16
>  t(Xh1)%*%eps
       [,1]
99.8445
Xh 947.8570

 solve(t(Xh1)%*%X1)
 t(Xh1)%*%(20*Q)
 t(Xh1)%*%(100*S)
 t(Xh1)%*%eps

であり、後段に比べて、前段がずいぶん大きくなっていることがわかる。

一般式

一般的に、回帰方程式

$$Y=X\beta+Z\gamma+\varepsilon$$

を考える。ただし、$X$は内生変数、$Z$は外生変数で$E(Z)=0$。
ここで、$Z$が観測できなかった場合を考える。
この場合、$Z\gamma+\varepsilon$を誤差項として扱わなければならないことになる。

$X$と$Z$の間に相関がなければ、OLS推定量は、
$$\hat{\beta}=({^tX}X)^{-1}{^tX}Y$$

もともと$Y=X\beta+Z\gamma+\varepsilon$なのだから、
$$\hat{\beta}=\beta+({^tX}X)^{-1}{^tX}Z+({^tX}X)^{-1}{^tX}\varepsilon$$

したがって、$plime\hat{\beta}$は、後ろの2項のplimが0なら$\beta$に等しくなり、一致性を持つ。

$plim{^tX}\varepsilon=0$のはずだから問題ないが、$plim{^tX}Z=0$とは限らない。

もし、$plim{^tX}Z=0$なら一致性を持つから、$X$が内生変数だろうとなかろうと気にすることはない。

ところが、$plim{^tX}Z\neq 0$なら、そのぶんがバイアスとなる。

$X$と$Z$に相関があるのなら、$X$を$Z$で回帰してみる。
$$X=Z\delta_1+V\delta_2+\eta$$

ただし、$V$は$Y$と無関係な$X$の説明変数。
$Z$は観察できないが、$V$なら観察できるとすると、操作変数法が使える。
$X$を$V$でOLS回帰して予測値$\hat{X}$を得る。
$$\hat{X}=V\hat{\delta}$$

これを使って
$$\hat{\beta}_{IV}=({^t\hat{X}}X)^{-1}{^t\hat{X}}Y$$

が得られる。このとき、
$$plim\hat{\beta}_{IV})=\beta+plim({^t\hat{X}}X)^{-1}{^t\hat{X}}Z\gamma+plim({^t\hat{X}}X)^{-1}{^t\hat{X}}\varepsilon
$$

であるから、$Z$と$\hat{X}$(つまり$V$)が無相関で、$\varepsilon$と$\hat{X}$も無相関なら、一致推定量が得られる。ただし、$V$と$X$の相関が弱い場合は、微妙になる。

このことを、$X$が定数項+1変数の場合で考えてみる。
$$
({^t\hat{X}}X)^{-1}
=\left(
\left[
\begin{array}{ccc}
1 & \cdots & 1 \
\hat{X}_1 & \cdots & \hat{X}_N \
\end{array}
\right]
\left[
\begin{array}{cc}
1 & X_1 \
\vdots & \vdots \
1 & X_N \
\end{array}
\right]
\right)^{-1}
$$

$$
=\left[
\begin{array}{cc}
N & \sum_iX_i \
\sum_i\hat{X}_i & \sum_i\hat{X}_iX_i \
\end{array}
\right]^{-1}
$$

$$
=\frac{1}{N\sum_i\hat{X}_iX_i-\sum_iX_i\sum_i\hat{X}_i}
\left[
\begin{array}{cc}
\sum_i\hat{X}_iX_i & -\sum_iX_i \ -\sum_i\hat{X}_i & N \
\end{array}
\right]
$$

$X$と$\hat{X}$の相関が0に近いと、分母の

$$
N\sum_i\hat{X}_iX_i-\sum_iX_i\sum_i\hat{X}_i
$$

が0に近いということになるので、$({^t\hat{X}}X)^{-1}$の値が大きくなってしまう。したがって、${^t\hat{X}}Z\gamma$や${^t\hat{X}}\varepsilon$が小さくても、$plim({^t\hat{X}}X)^{-1}{^t\hat{X}}(Z\gamma+\varepsilon)$全体では0にならないことがあり、それがそのままバイアスとなる。

二段階最小二乗法(TSLS)

データ$Q$も使えるなら、$X$を$Q$と$K$で回帰して、

$$
\hat{X}=\hat{\delta}_0+\hat{\delta}_1Q+\hat{\delta}_3K
$$

を使うこともできるが、弱相関の$K$ひとつを操作変数とするよりもましになるというわけでもなさそう。

> o02<-lm(X~Q+K)
>  Xh<-predict(o02)
>
>  X1<-cbind(matrix(1,N,1),X,Q)
>  Xh1<-cbind(matrix(1,N,1),Xh,Q)
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
       [,1]
250.02832
X  29.70317
Q  17.23184
>
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(300+25*X)
          [,1]
3.000000e+02
X 2.500000e+01
Q 3.672444e-12
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(200*S)
         [,1]
-104.990877
X    9.894105
Q   -6.014914
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%eps
        [,1]
2.5237632
X -0.2438844
Q  0.2392950

 o02<-lm(X~Q+K)
 Xh<-predict(o02)
 
 X1<-cbind(matrix(1,N,1),X,Q)
 Xh1<-cbind(matrix(1,N,1),Xh,Q)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
 
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(300+25*X)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%(200*S)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%eps

だからといって、$K$が不要なわけではない。
$K$がなくて、$W=300+25X+...$の式の中にある$Q$だけで$X$を回帰してしまうと、$\hat{X}$と$Q$が完全に相関してしまうので、${^t\hat{X}}X$の逆行列が解けなくなってしまう。

> o02<-lm(X~Q)
>  Xh<-predict(o02)
>
>  X1<-cbind(matrix(1,N,1),X,Q)
>  Xh1<-cbind(matrix(1,N,1),Xh,Q)
>  solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
 以下にエラー solve.default(t(Xh1) %*% X1) :
システムは数値的に特異です:条件数の逆数 = 4.29402e-018
>
>  cor(Xh,Q)
[1] 1

 o02<-lm(X~Q)
 Xh<-predict(o02)
 
 X1<-cbind(matrix(1,N,1),X,Q)
 Xh1<-cbind(matrix(1,N,1),Xh,Q)
 solve(t(Xh1)%*%X1)%*%t(Xh1)%*%w
 
 cor(Xh,Q)

これは$Q$と$S$両方使っても同じこと。$X$の説明変数には$w$の説明変数以外のものも含まれていないとうまくいかんということ?

tslsコマンドで解いてみる

> library(sem)
>  summary(tsls(w~X+Q,~Q+K))

2SLS Estimates

Model Formula: w ~ X + Q

Instruments: ~Q + K

Residuals:
Min.   1st Qu.    Median      Mean   3rd Qu.      Max.
-9.54e+01 -7.39e+01 -6.12e+01  1.98e-11  7.91e+01  9.95e+01

Estimate Std. Error t value  Pr(>|t|)
(Intercept)   250.03     83.205   3.005 2.722e-03
X              29.70      8.481   3.502 4.815e-04
Q              17.23      4.037   4.268 2.158e-05

Residual standard error: 77.1049 on 997 degrees of freedom

>  summary(tsls(w~X+Q,~Q))
 以下にエラー chol.default(XtZ %*% invZtZ %*% t(XtZ)) :
次数 3 の主対角行列が正定値ではありません

 library(sem)
 summary(tsls(w~X+Q,~Q+K))
 summary(tsls(w~X+Q,~Q))