この記事では計量経済学の分析手法の一つである, 「操作変数法」について説明していきます。操作変数法は, 回帰モデルを正しく推定するために用いられ, 「内生性」の問題を解決するために役立つとされています。まずは内生変数と外生変数とは何か, というところから確認をしていきましょう。

内生変数と外生変数

この二つを一言で表すと, 内生変数とは誤差項と相関を持つ説明変数であり, 外生変数は誤差項と相関を持たない説明変数です。これだけでは何のことかよくわかりませんよね。数式で確認していきましょう。

以下のような単回帰モデルを考えてみます。

Y_i = \beta_0 + X_i\beta_i + \epsilon_i, i = 1, 2, ..., n\quad(1)

ここで, X_iが説明変数で, \epsilon_iは誤差項です。また, E[\epsilon_i] = 0とします。(単回帰分析の誤差の仮定より)

外生変数とは誤差項と相関を持たない変数でした。ここで, 相関係数が0ならば共分散も0となります。 X_iと\epsilon_iの共分散は以下のような式になることからも0であることが確認できます。(相関係数の出し方がわからない方はこちらの記事など参考にしてみてください)

cov(X_i, \epsilon_i)\\ =E[X_i, \epsilon_i] - E[X_i]E[\epsilon_i]\\ =E[X_i]E[\epsilon_i] - E[X_i]E[\epsilon_i] (Xとεは独立なため)\\ =0 (E[\epsilon_i] = 0より)

これが成り立っていれば, X_iは外生変数と言えます。
一方, 内生変数は共分散が0ではないということになるので,
cov(X_i, \epsilon_i) ≠ 0
となります。

内生変数によって生じる問題

内生性があることによって回帰係数が正しく推定できない, という問題が発生します。

まず, 外生変数である場合を考えてみます。
(1)の両辺の期待値をとると,

E[Y_i] = \beta_0 + \beta_1E[X_i] + E[\epsilon_i]\\ \beta_0について解くと,\\ \beta_0 = E[Y_i] -\beta_1E[X_i]\quad(E[\epsilon_i] = 0より)\\ よって, (1)は\\ Y_i = E[Y_i] -\beta_1E[X_i] + X_i\beta_i + \epsilon_i\\ Y_i - E[Y_i] = (X_i - E[X_i])\beta_1 + \epsilon_i
となります。さらに両辺に X_iをかけて, 期待値をとると,

E[X_i(Y_i - E[Y_i]) = E[X_i(X_i - E[X_i])]\beta_1 + E[X_i\epsilon_i]\\ Cov(X_i, Y_i) = \beta_1Var(X_i) + Cov(X_i, \epsilon_i)\quad(2)\\ 両辺をVar(X_i)で割ると, \\ \beta_1 = \frac{Cov(X_i, Y_i)}{Var(X_i)}\quad(Cov(X_i, \epsilon_i) = 0より)\quad(3)

よって, それぞれを標本分散, 標本共分散に置き換えることで\beta_1を実際のデータX,Yから正しく推定することができます。

では今度は内生変数の場合を考えます。先ほどと異なる点はCov(X_i, \epsilon_i) \neq 0です。同様に(2)の両辺をVar(X_i)で割ると,

\frac{Cov(X_i, Y_i)}{Var(X_i)} = \beta_1 + \frac{Cov(X_i, \epsilon_i)}{Var(X_i)}
となってしまい, 誤差が生じます。この式の中では\frac{Cov(X_i, \epsilon_i)}{Var(X_i)}が該当し, 内生変数バイアスと呼ばれています。(誤差がプラスかマイナスかは説明変数と被説明変数の相関が正か負かによります)

正しく効果を推定するにはこのバイアスをできるかぎり取り除く必要があります。そこで用いられるのが操作変数法です。

操作変数法の考え方

操作変数を一言であらわすと「内生変数との相関があり回帰モデルの誤差項と相関を持たないもの」となります。
具体例を考えてみましょう。
学歴を説明変数, 賃金を被説明変数として学歴が賃金に及ぼす影響を分析するとします。しかし学歴は内生変数である可能性が高いです。何故なら賃金のその他の決定要因は全て誤差項に含まれるわけですが, 個人の能力などは学歴と相関があると考えられます。従ってそのまま分析すると先ほどのバイアスが発生してしまいます。ここで操作変数として, 例えば「親の学歴」などを加えることが想定されます。これは本人が関与できない, という部分がポイントです。子どもの意思とは関係なく親が大学を出ていれば自分の子どもにも大学に行って欲しいと考える可能性が高そうだからです。

操作変数の導出

なんとなくイメージが掴めたところで数式的に操作変数を検討してみます。
操作変数をZ_iとします。これを式で表すと,
X_i = \gamma_0 + Z_i\gamma_1 + u_i\quad(4)\\ E[Z_i\epsilon_1] = Cov(Z_i, \epsilon_1) = 0\\ ただしu_iは誤差項でE[u_i]=0, E[Z_iu_i]=0\\
この式でZ_iがX_iに影響を与えていることと, 回帰モデルの誤差項と相関を持たないことを示せています。
さらに(1)を変形すると,
Y_i = \beta_0 + (\gamma_0+Z_i\gamma_1 + u_i)\beta_1 + \epsilon_1\\ \quad = \beta_0 + (\gamma_0+Z_i\gamma_1)\beta_1 + u_i\beta_1 + \epsilon_i\\ \quad = \beta_0 + (\gamma_0+Z_i\gamma_1)\beta_1 + e_i
となり, u_i\beta_1 + \epsilon_iは元々の誤差と説明変数を置き換えた誤差を加えたものでこれを, e_iとします。ここで仮定よりE[e_i] = 0です。

外生変数の時のようにただしく推定するには誤差項と説明変数が相関を持たないことが重要でした。つまりここでは(\gamma_0+Z_i\gamma_1)とe_iが無相関である必要があります。つまり共分散が0となればいいので, そうなるかを確認してみます。
Cov((\gamma_0+Z_i\gamma_1), e_i)\\ = E[(\gamma_0+Z_i\gamma_1)e_i]\\ = E[\gamma_0e_i]+E[Z_i\gamma_1e_i]\\ = \gamma_0E[e_i]+E[(Z_i\gamma_1)(u_i\beta_1 + \epsilon_i)]\quad(E[e_i]=0)\\ = 0 + \gamma_1\beta_1E[Z_iu_i] + \gamma_1E[Z_i\epsilon_1]\quad(E[Z_iu_i]=0, E[Z_i\epsilon_1]=0)\\ = 0

従って, 内生変数バイアスが操作変数によって回避されます。実際\beta_1は,内生変数バイアスが消えることで,
\beta_1 = \frac{Cov(\gamma_0+Z_i\gamma_1, Y_i)}{Var(\gamma_0+Z_i\gamma_1)}\quad(X_i = \gamma_0+Z_i\gamma_1より)

操作変数法の実行

操作変数法を実施するにあたり, 2段階最小二乗法という手法を用います。名前のとおり最小二乗法を2回行うわけですが, ステップは以下の通りです。

1. (4)を最小二乗法で推定して, X_iの予測値\hat{X_i} = \hat{\gamma_0}+\hat{\gamma_1}Z_iを得る
2. (1)のX_iを\hat{X_i}で置き換えて, 最小二乗法で\beta_1を推定する

これを実際のデータを用いて実行してみます。
データはWooldridgeのサイトのmorzというデータを使用します。このデータには個人の賃金や, 教育を受けた年数, 両親の教育を受けた年数などが入っています。

ここからはRを使っていきます。Rでの2段階最小2乗法はsemパッケージのtsls関数を用います。

# 必要なライブラリの読み込み
library(sem)
## stataファイルを読み込むのに使う
library(readstata13)
library(tidyverse)
mroz<-read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/mroz.dta") %>% 
  ## wage(賃金)が欠損値のデータを除去
  na.omit(wage)

ここで教育年数と職業の経験年数が賃金にどの程度影響を及ぼしているのかに関心があるとします。しかし, 先述の通り, 教育年数は内生変数である可能性が高いです。そこで今回は両親の教育年数を操作変数に設定し, 分析を行ってみます。

結果の違いを確認するために通常の回帰分析も併せて実行してみます。
まず, 全ての説明変数を外生変数として分析した結果が以下です。

 # 全て外生変数とした回帰分析
lm <- lm(wage ~ educ + exper + fatheduc + motheduc,data = mroz)
summary(lm)

教育年数の回帰係数は約0.54と算出されました。
では次に2段階最小二乗法で分析を行います。ここでは両親の教育年数(motheduc, fatheduc), 自身の職業経験年数(exper)を操作変数として設定します。

 # 2段階最小二乗法
lm_tsls <- tsls(wage ~ educ + exper, instruments = ~ exper + fatheduc + motheduc,
                data = mroz)
summary(lm_tsls)

教育年数の回帰係数は約0.33と算出されました。
以上の結果より操作変数を設定しないと効果が過剰に推定されてしまう可能性がある, ということがわかりました。

有効な操作変数の見つけ方

それでは最後に有効な操作変数の見つけ方を考えたいと思います。

操作変数は, ただ内生変数と相関があるだけでなく, 強い相関を持つほうが望ましい。ということです。
操作変数の強弱の確認方法の一つとして, (4)での\gamma_1の推定値が有意に0と異なるかどうか, があります。細かい説明は省きますが, もし\gamma_1 = 0の場合, \beta_1の推定量として意味のある値をとらなくなってしまいます。

さきほどの通常の回帰分析の結果では, 両親の教育年数(motheduc, fatheduc)は有意に0であるとはいえない結果となっているので, 弱操作変数の疑いがあります。

この場合は異なる操作変数を探す必要がある, ということになります。
適切な操作変数を見つけることはなかなか難しいことのように感じました。

参考文献

星野・田中(2016) Rによる実証分析 - 回帰分析から因果分析へ
末石(2015) 計量経済学 ミクロデータ分析へのいざない