一般化モーメント法と操作変数法 - データサイエンス関連+αのメモ

母集団モーメントを利用する母回帰係数の表記

単回帰モデル $Y_i = \alpha + \beta X + \gamma Z_i + u_i$ のとき、

\beta_{IV} = \frac{Cov(Z_i, Y_i)}{ Cov(Z_i, X_i) } , \hspace{2em} \alpha_{IV} = E(Y_i) - \beta_{IV} E(X_i)

(1)

モーメント条件から導出するOLS¶

IV推定量（行列表記）¶

同時方程式モデル

Y_1 = \beta_1 Y_2 + \gamma_1 Z_1 + \varepsilon_1\\ Y_2 = \beta_2 Y_1 + \gamma_2 Z_2 + \varepsilon_2

(2)

があるとする。外生変数が $Z$ である。

行列形式で $\boldsymbol{\beta}_1 = (\beta_1, \gamma_1)^T$ のようにまとめて

\renewcommand{\b} when command \b does not yet exist; use \newcommand

\renewcommand{\b}[1]{\boldsymbol{#1}}
\b{Y}_1 = [\b{Y}_2|\b{Z}_1] \b{\beta}_1 + \b{\varepsilon}_1 = \b{X}_1 \b{\beta}_1 + \b{\varepsilon}_1\\
\b{Y}_2 = [\b{Y}_1|\b{Z}_2] \b{\beta}_2 + \b{\varepsilon}_2 = \b{X}_2 \b{\beta}_2 + \b{\varepsilon}_2

と書くことにする。外生変数からなる行列を\b{Z} = [\b{Z}_1|\b{Z}_2]と表し、転置行列を左から掛けると

Undefined control sequence: \b at position 1: \̲b̲{Z}^T \b{Y}_1 =…

\b{Z}^T \b{Y}_1 = \b{Z}^T \b{X}_1 \b{\beta}_1 + \b{Z}^T \b{\varepsilon}_1

\b{Z}^T \b{X}_1の逆行列の存在を仮定して掛けると

Undefined control sequence: \b at position 2: (\̲b̲{Z}^T \b{X}_1)^…

(\b{Z}^T \b{X}_1)^{-1} \b{Z}^T \b{Y}_1 =  \b{\beta}_1
+ (\b{Z}^T \b{X}_1)^{-1} \b{Z}^T \b{\varepsilon}_1

標本モーメントを使って表すと

Undefined control sequence: \b at position 15: \left( \frac{ \̲b̲{ Z^T X}_1 }{n}…

\left( \frac{ \b{ Z^T X}_1 }{n} \right) ^{-1}
\left( \frac{ \b{ Z^T Y}_1 }{n} \right)
=  \b{\beta}_1
+ \left( \frac{ \b{ Z^T X}_1 }{n} \right)^{-1}
  \left( \frac{ \b{Z}^T \b{\varepsilon}_1 }{n} \right)

\b{Z^T X}_1/nは定数に、\b{Z}^T \b{\varepsilon}_1/nは\b{Z}が誤差項\b{\varepsilon}と無相関であれば0に確率収束する。

Undefined control sequence: \b at position 1: \̲b̲{b}_1^{(IV)} = …

\b{b}_1^{(IV)} = (\b{ Z^T X}_1)^{-1} \b{ Z^T Y}_1

を**操作変数推定量（instrumental variable estimator）**という。

\b{Y}_1 = \b{X}_1 \b{\beta}_1 + \b{\varepsilon}_1を代入すると

Undefined control sequence: \b at position 15: \begin{align}
\̲b̲{b}_1^{(IV)}
&=…

\begin{align}
\b{b}_1^{(IV)}
&= (\b{ Z^T X}_1)^{-1} \b{Z}^T (\b{X}_1 \b{\beta}_1 + \b{\varepsilon}_1)\\
&= (\b{ Z^T X}_1)^{-1} \b{Z}^T \b{X}_1 \b{\beta}_1
  + (\b{ Z^T X}_1)^{-1} \b{Z}^T \b{\varepsilon}_1\\
&= \b{\beta}_1 + (\b{ Z^T X}_1)^{-1} \b{Z}^T \b{\varepsilon}_1\\
\end{align}
\DeclareMathOperator*{\plim}{plim}

となる。確率極限をとると\plim_{n\to\infty} \b{Z}^T \b{\varepsilon}_1/n = 0なので

Undefined control sequence: \DeclareMathOperator at position 1: \̲D̲e̲c̲l̲a̲r̲e̲M̲a̲t̲h̲O̲p̲e̲r̲a̲t̲o̲r̲*{\plim}{plim}
…

\DeclareMathOperator*{\plim}{plim}
\plim_{n\to\infty} \b{b}_1^{(IV)}
= \b{\beta}_1

標本モーメントの確率極限

$n$ を大きくしていったとき、標本平均 $\sum_i X_i / n$ は母平均 $\mu$ に一致する。これを確率極限といい、 $\text{plim}(\sum_i X_i / n) = \mu$ のように表す。

\text{ plim }\left( \frac{\sum_i X_i^2}{n} \right) = \text{ plim }\left( \frac{\boldsymbol{x}^T \boldsymbol{x}}{n} \right) = \mu^2 - \sigma^2

(10)

（ $Var[X] = E[X^2] - E[X]^2$ なので $\mu^2$ が残っている)

2つの変数の2次の標本モーメントは

\text{ plim }\left( \frac{\sum_i X_i Y_i}{n} \right) = \text{ plim }\left( \frac{\boldsymbol{x}^T \boldsymbol{y}}{n} \right) = \mu_X \mu_Y - \sigma_{XY}

(11)

となる。平均がゼロで $X$ と無相関の $\varepsilon$ との2次の標本モーメントは

\text{ plim }\left( \frac{\sum_i X_i \varepsilon_i}{n} \right) = \text{ plim }\left( \frac{\boldsymbol{x}^T \boldsymbol{\varepsilon}}{n} \right) = \mu_X \mu_{\varepsilon} - \sigma_{X \varepsilon} = 0

(12)

となる

# Generate Data
import numpy as np
np.random.seed(0)
n = 10
b = np.array([0.3, 0.5])

z = np.array([])
x = np.array([])
y = np.array([])
for i in range(10):
    z_ = np.random.uniform(low=0, high=100, size=n)
    if i == 0:
        x_ = np.zeros(shape=(n,))
    else:
        x_ = y_

    y_ = b[0] * x_ + b[1] * z_ + np.random.normal(size=n)
    z = np.concatenate([z, z_])
    x = np.concatenate([x, x_])
    y = np.concatenate([y, y_])

Z = z.reshape(n*10, -1)
X = x.reshape(n*10, -1)

import seaborn as sns
import pandas as pd
df = pd.DataFrame(dict(y=y, x=x, z=z))
sns.pairplot(df, height=1.5, y_vars=["y", "x"])

<seaborn.axisgrid.PairGrid at 0x7fe0f848a070>

np.linalg.inv(Z.T @ X) @ (Z.T @ y)

array([1.4245943])

# Generate Data
import numpy as np
n = 500
np.random.seed(0)
b_zx = 0.5
b_xy = 0.8
z = np.random.uniform(low=0, high=100, size=n)
u = np.random.normal(scale=10, size=n) # unobserved variable
x = u + b_zx * z + np.random.normal(scale=5, size=n) # Z -> X
y = u + b_xy * x + np.random.normal(size=n) # X -> Y; U -> {X, Y}

X = x.reshape(n, -1)
Z = z.reshape(n, -1)


import seaborn as sns
import pandas as pd
df = pd.DataFrame(dict(y=y, x=x, z=z, u=u, e=e))
sns.pairplot(df, height=1.5, y_vars=["y", "x"])

<seaborn.axisgrid.PairGrid at 0x7fe0f41ef550>

# OLS
np.linalg.inv(X.T @ X) @ (X.T @ y)

array([0.87772517])

# IV
z2y = np.linalg.inv(Z.T @ Z) @ (Z.T @ y)
z2x = np.linalg.inv(Z.T @ Z) @ (Z.T @ x)
z2y / z2x

array([0.77019506])

z2x

array([0.49124653])

# IV
np.linalg.inv(Z.T @ X) @ (Z.T @ y)

array([0.77019506])

from linearmodels.iv import IV2SLS
model = IV2SLS.from_formula("y ~ [x ~ z]", df).fit()
model

Loading...

2SLS推定量とIV推定量の同値性¶

2SLS推定量とIV推定量は同値である（『新しい計量経済学』p.216）