操作変数法と２段階最小二乗法 - 知らないものは知らないを通す

概要

『空間統計学 - 自然科学から人文・社会科学まで』(朝倉書店）という書籍を読んでいる. その中の数学的準備において, 個人的に馴染みのない複数の統計分析手法・モデルが出てきた.　本当は「空間統計学」の学びを共有したいと思っているが, その前提知識が欠けているため, 勉強がてら知識の整理を行う

今回扱うモデルは, 以下のものである.

1. 操作変数法(Instrument variables; IV)法
2. 二段階最小二乗法(two stage least squares; 2SLS)
3. 一般化モーメント法(generalized method of moments; GMM)(ポストが長くなったので次回です😢)

https://www.asakura.co.jp/detail.php?book_code=12831

回帰モデルの課題

まず基本的な型として, 線形回帰モデルを復習する. いま, 観測データ $y_i (i = 1,2,...,h-1)$に対し, 次の回帰式を考える.

$$y_i = \beta_0 + \sum_{k=1}^{h-1} x_{h,i} \beta_h + \epsilon_i$$

ここで, $y_i$は従属変数, $x_{h,i} (h=1,2,...,k-1)$は（外生）説明変数, $\beta_0$は定数項, $\beta_h$は$x\_{h,i}$に対応する回帰係数, そして$\epsilon_i$は誤差項である. さらに, 誤差項$\epsilon_i$は, 平均0,分散$\sigma^{2}$の正規分布に従うとする.$\epsilon ~ N(0, \sigma^{2}$これを観測データ全体を並べた式として, 行列表現で表す.

$$\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \beta_0 \begin{pmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{pmatrix} + \begin{pmatrix} x_{1,1} & x_{2,1} & \cdots & x_{k-1,1} \\ x_{1,2} & x_{2,2} & \cdots & x_{k-1,2} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1,n} & x_{2,n} & \cdots & x_{k-1,n} \end{pmatrix} \begin{pmatrix} \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix}$$

この式を簡略化し, $y = \beta_0 1 + x \beta_1 + \epsilon$と表現することもある. しかし, もう少し便利な書き方として, 次の書き方に直す.

$$\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_{1,1} & x_{2,1} & \cdots & x_{k-1,1} \\ 1 & x_{1,2} & x_{2,2} & \cdots & x_{k-1,2} \\ \vdots & \vdots & \ddots & \vdots \\ 1 & x_{1,n} & x_{2,n} & \cdots & x_{k-1,n} \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_n \end{pmatrix} + \begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix}$$

すると, この回帰式は, $y = X\beta + \epsilon$と, とてもシンプルに書ける.

次に, この式を用いて, 観測データ$y$から, 説明変数$X$と回帰係数$\beta$を推定することを考える.

最小二乗法では, 回帰式において回帰係数$\beta$の推定値$\hat{\beta}$を使い, 予測誤差$\hat{\epsilon}$の二乗和を最小にすることを考える.
つまり, $\hat{\epsilon}^{T}\hat{\epsilon} = \hat{\sigma}^{2} = \parallel y - X\hat{\beta} \parallel_{2}$を求める. 最適化の一階条件により, 両辺を$\hat{\beta}$で微分した式を 0とおくと, $\hat{\beta}$を求めることができる.

$$\begin{align} \hat{\sigma}^{2} &= \parallel y - X\hat{\beta} \parallel_{2} \\ &= y^{T}y - y^{T}X\hat{\beta} - (X\hat{\beta})^{T}y + X\hat{\beta}(X\hat{\beta})^{T} \\ &= y^{T}y - y^{T}X\hat{\beta} - \hat{\beta}^{T}X^{T}y + X\hat{\beta}\hat{\beta}^{T}X^{T} \\ \frac{\partial \hat{\sigma}^{2}}{\partial \hat{\beta}} &= -X^{T}y - X^{T}y + 2X^{T}X\hat{\beta} = 0 \Rightarrow X^{T}y &= X^{T}X\hat{\beta} \\ \Rightarrow \hat{\beta} &= (X^{T}X)^{-1}X^{T}y \end{align}$$

特にこの基本モデルでは次の仮定がおかれている(一部だけ提示）

1. $X$は確定的である:　値が固定されている
2. $X$が与えられたとき, $y$の条件付き期待値$E[y|X]$は, $X\beta$であり, $ \epsilon]の期待値は0, つまり$E[\epsilon|X] = 0$
3. 誤差項$\epsilon$は, 均一分散性と共分散がゼロである: 各観測データに対する誤差項$\epsilon_i$は互いに無相間

しかし,実際のデータに触れると,これらの仮定がすべて成り立つケースは少ない. そのため, そのまま最小二乗法を適用することが不適切な場合が多くなる.
そのための対策が必要になってくる, つまり, この欠落変数バイアスをどのように回避するか,,,

操作変数法(IV法）

上記の仮定のうち, 特に二番目の説明変数$:X$と 誤差項$\epsilon$が無相関の場合, 先ほどのOLS推定量は, 一致性も不遍性ももたない^[1]

ここでは, なかでも欠落変数バイアスと呼ばれる, 未観測の説明変数の存在によるバイアスの影響を考える.
以下の問題設定は,著書『因果推論入門 ~ ミックステープ：基礎から現代的アプローチまで』(技術評論社）を参照している.

欠落変数バイアス

次の問題を考えたい. 「教育年数$S$」と「収入$Y$」の因果関係について知りたいが, 実はこの「教育年数」が未観測の変数(例：個人の能力$A$）によって内生的である, という古典的な労働経済学の問題である.

これを, シンプルな回帰式で表現すると以下の通り(これを真のモデルとする）.

$$Y_i = \alpha + \delta S_i + \gamma A_i + \epsilon_i$$

$\epsilon$は, 「教育年数$S$」と「個人の能力$A$」とは無相間の誤差項である. つまり, $E[ \varepsilon S] = 0, E[\varepsilon A] = 0$である.

しかし, 何らかの原因により, 変数$A$が未観測となった場合, （仕方なく）以下のモデルで代替したとする.

$$Y_i = \alpha + \delta S_i + \eta_i$$

ここで,$\eta_i$は, 真のモデルでの$\gamma A_i + \varepsilon_i$に等しい複数の項からなる誤差項である.

教育年数$S$と, 個人の能力$A$には相関があると仮定しているため, $\eta$と$S$も相関しているとみなす.
よって, 最小二乗法の導出から, 教育年数$S$にかかる係数$\delta$は,

$$\hat{\delta} = \frac{C(Y,S)}{V(S)} = \frac{E[YS] - E[Y]E[S]}{V(S)}$$

で求まる.
ここで, $Y$を, 真のモデルで表現した場合の式を代入すると,

$$\begin{align} \hat{\delta} &= \frac{E[\alpha + S + S^2\delta + \gamma SA + \varepsilon S] - E(S)E[\alpha + \delta S + \gamma A + \varepsilon]}{V(S)} \\ &= \frac{\delta E(S^2) - \delta E(S)^2 + \gamma E(AS) - \gamma E(S)E(A) + E(\varepsilon S) - E(S)E(\varepsilon)}{V(S)} \\ &=\delta + \gamma \frac{C(A,S)}{V(S)} \\ \end{align}$$

この式をみればわかるが, 代替モデル（変数が少ないほう）を用いた場合の$\hat{\delta}$は, 真のモデルで求められる$\delta$より, $\gamma \frac{C(A,S)}{V(S)}$だけ, バイアスがかかっていることがわかる. 実際は, $A$が未観測であるため, このバイアスの値は計算できないが, 理論上, 未観測変数と, （内生）変数とが強い相関をもつほど, より大きなバイアスがかかった結果を推定しかねない, ということである.

操作変数の導入

上記のバイアスを回避しつつ, うまく$\delta$を計算する方法はなりだろうか？　ここでひとつのテクニックとして, 操作変数(nstrument variables; IV) を導入することを考える. この操作変数$Z$は, 未観測変数（個人の能力$A$)や,誤差項$\varepsilon_i$とは独立で, 他方で, 教育年数$S$と相関をもつ, とても都合のよいものが発見されたとする. この操作変数を使って, 収入$Y$との共分散を計算してみると,

$$\begin {align} C(Y,Z) &= C(\alpha + \delta S + \gamma A + \varepsilon, Z) \\ &= E[(\alpha + \delta S + \gamma A + \varepsilon)Z] + E(Y) E(Z) \\ &= \{\alpha E(Z) - \alpha E(Z) \} + \delta \{ E(SZ) - E(S)E(Z)\} \\ &+ \gamma \{E(AZ) - E(A)E(Z)\} + \{E(\varepsilon Z) - E(\varepsilon) E(Z) \} \\ &= \delta C(S, Z) + \gamma C(A, Z) + C(\varepsilon, Z) \\ \end {align}$$

右辺の第二項,　第三項は, 操作変数$Z$が個人の能力$A$と誤差項$\varepsilon_i$とは独立,であるため 0 となる. 結果として, 第一項のみが残る.
この式から$\delta$を計算するのはとても簡単で,

$$\hat{\delta} = \frac{C(Y,Z)}{C(S,Z)}$$

で求めることができる. 式をみるとわかるが, 未観測変数$A$が含まれていないことから, バイアスを含まない推定値を, 未観測変数に依存せずに求めることができるのだ.

このテクニックの肝は, 以下の仮定を満たす操作変数をうまく探すことにかかっている

1. 未観測変数および誤差項とは独立であること
2. 内生変数と強い相関があること((弱い相関でもよいとなれば, それは単なる乱数でも機能してしまうため))

２段階最小二乗法

操作変数法(IV法）により, 未観測変数による欠落変数バイアスの影響を無視できることがわかった, 次に, ２段階最小二乗法(two stage least squares; 2SLS) により, この操作変数推定量（操作変数を導入した場合の回帰係数の推定量）を求めていく.

さきほどの例にもとづき, 教育年数$S$, 収入$Y$,そして　操作変数$Z$のデータの標本があるとし, 次の２つの式により各種データが生成されているとする.

$$\begin{align} Y_i &= \alpha + \delta S_i + \varepsilon_I \\ S_i &= \gamma + \beta Z_i + \epsilon_i \\ \end{align}$$

ここで仮定として, 1) $C(Z, \varepsilon) = 0$であること, 2) $\beta \neq 0$である, こととしている.

二つの式を, さきほどの $\hat{\delta}$に真のモデル$Y_i = \alpha + \delta S_i + \gamma A_i + \varepsilon_i$代入すると,

$$\begin{align} \hat{\delta} &= \frac{C(Y,Z)}{C(S,Z)} =\frac{C(\alpha + \delta S + \gamma A + \varepsilon, Z}{S, Z}\\ &=\frac{1}{C(S,Z)} [\beta C(S, Z) L \gamma C(A, Z) + C(\varepsilon, Z) ] \\ &= \beta \end{align}$$

となる. 最後の式は, 除外制約($C(A,Z) =0$ かつ $C(\epsilon, Z) = 0$)により.

そして, 上式に$Y_i$に関する回帰式を代入すると,

そして $\hat{\delta}$は次のように表現することもできる

$$\begin{align} \hat{\delta} &= \frac{C(Y,Z)}{C(S,Z)} \\ &= \frac{\frac{C(Z,Y)}{V(Z)}} {\frac{C(Z, S)}{V(Z)} } \\ \end{align}$$

つまりは, 操作変数法っで得られる推定量は, 二つの回帰係数の比で表される^[2].

ここで, 分母は, $S_i = \gamma + \beta Z_i + \epsilon_i$から, $\hat{\beta}$と等しくなる. つまり,

$$\begin{align} \hat{\beta} &= \frac{C(Z, S)}{V(Z)} \\ \hat{\beta}V(Z) &= C(Z, S) \\ \end{align}$$

この関係式を先ほどの式に代入すると

$$\begin{align} \hat{\delta} &= \frac{C(Z,Y)} {C(Z, S)} \\ &= \frac{\hat{\beta}C(Z,Y)} {\hat{\beta}C(Z, S)} \\ &= \frac{\hat{\beta}C(Z,Y)} {\hat{\beta}^{2}V(Z)} \\ &= \frac{C(\hat{\beta}Z,Y)} {V(\hat{\beta}Z)} (\text{それぞれ, 共分散, 分散の乗法の関係を利用})\\ \end{align}$$

最後は, $V(aX) = a^{2} V(X), C(aX, Y) = aC(X,Y)$($a$は定数, $X,Y$は変数）の関係式を利用.
ここで, $\hat{\beta}$は何であるかというと, 操作変数$Z$を用いた回帰式$S_i = \gamma + \beta Z_i + \epsilon_i$で推定された$\hat{\beta}$を用いた変数$S$への当てはめ値である. つまり, $\hat{S} = \hat{\gamma} + \hat{\beta}Z$を用いると,

$$\begin{align} \hat{\delta}_{IV} &= \frac{C(\hat{\beta}Z,Y)} {V(\hat{\beta}Z)} \\ &= \frac{C(\hat{S},Y)} {V(\hat{S}Z)} \end{align}$$

となっている.
式を導出すると,

$$\begin{align} C(\hat{S}, Y) &= E[\hat{S} Y] - E [\hat{S} ] E[Y ] \\ &= E(Y[\hat{\gamma} + \hat{\beta}Z ]) - E(Y)E(\hat{\gamma} + \hat{\beta}) \\ &= \hat{\gamma}E(Y) + \beta E(YZ) - \hat{\gamma}E(Y) - \beta E(Y)E(Z) \\ &= \hat{\beta}(E(Y,Z) - E(Y)E(Z)) \\ &= \hat{\beta} C(Y,Z) \end{align}$$

から

まとめ

ここでようやく２段階最小二乗法(2SLS) の意味や手順が判明した.
つまり,

1. 操作変数$Z$を用いて, ある内生変数$S$に対する回帰式モデルを作成し, その内生変数に対する推定値$\hat{S}$をつくる(１段階目)
2. 推定値$\hat{S}$を使った, 本来の従属変数$Y$への回帰式モデルを作成し, 回帰係数$ \beta_{2SLS}$を推定する　（２段階目）

操作変数法(IV法）と比較すると,

1. 操作変数法(IV法)では, 操作変数$Z$を使った２つの回帰式モデルで推定した回帰係数の比を使う
2. ２段階最小二乗法(2SLS)では, 操作変数$Z$を使った内生変数の推定値(1段階目)をもちいた,直接的な回帰係数の推定を行う（２段階目）

という違いがあると 個人的に理解している.

手法の利用Tipsとして, 操作変数の数 >= 内生変数の数であることが求められる. 特に, 操作変数と内生変数の数が一致する場合には, 操作変数法(IV法)と２段階最小二乗法（2SLS)の推定値は一致する^[3]

あと, ２段階最小二乗法（2SLS)をプログラムで実装する場合,１段階目と２段階目を別々に計算すると, バイアスのない回帰係数の推定値を得られるが, 標準誤差が正しくないらしく, 厳密には解が一致しない

余談だが, 計量経済学などで有名な2SLSだが, 空間統計学では代表的な空間計量経済学空間ラグモデル(spatial lag model; SLM) の推定法のひとつとなっている

この2SLSのより一般化したものとして, 一般化モーメント法(generalized method of moments; GMM) がある. これについては次回の記事で扱うとする.

参考サイト

以下のQiitaの記事は, 操作変数法のイメージがわかりやすい

https://qiita.com/icchi-ioo/items/2aa1f36dc73675932953

https://qiita.com/icchi-ioo/items/74da42b24988a5d18e65

２段階最小二乗法を含め, 全体的な操作変数法の枠組みをわかりやすく解説

https://yukiyanai.github.io/jp/classes/econometrics2/contents/slides/metrics2_topic09_slides.pdf

Rで実装する場合は, semライブラリの$tsls$関数や, AERパッケージの$ivreg$関数が利用されているようだ.
Rでのプログラム例は以下のものが参考にできる

https://user.keio.ac.jp/~nagakura/R/R_2OLS.pdf

https://info-zebra.com/2-stage-least-squares

脚注

1. このような問題がおきる例として, $x$が観測誤差を持つ場合, モデルに考慮されていない除外変数の存在が大きい場合, $x$と$y$の同時性がある場合などがある

2. たとえば, 次の二つの回帰モデルから, 分子($\rho$)と分母($\phi$)を推定できる. $S_i = \alpha_1 + \phi Z_i + \varepsilon_i, Y_i = \alpha_0 + \rho Z_i + \epsilon_i$

3. https://yukiyanai.github.io/jp/classes/econometrics2/contents/slides/metrics2_topic09_slides.pdf より