一般化モーメント法と操作変数法

一般化モーメント法と操作変数法#

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

単回帰モデルYi=α+βX+γZi+uiのとき、

βIV=Cov(Zi,Yi)Cov(Zi,Xi),αIV=E(Yi)−βIVE(Xi)

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

IV推定量（行列表記）#

同時方程式モデル

Y1=β1Y2+γ1Z1+ε1Y2=β2Y1+γ2Z2+ε2

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

行列形式でβ1=(β1,γ1)Tのようにまとめて

Y1=[Y2|Z1]β1+ε1=X1β1+ε1Y2=[Y1|Z2]β2+ε2=X2β2+ε2

と書くことにする。外生変数からなる行列をZ=[Z1|Z2]と表し、転置行列を左から掛けると

ZTY1=ZTX1β1+ZTε1

ZTX1の逆行列の存在を仮定して掛けると

(ZTX1)−1ZTY1=β1+(ZTX1)−1ZTε1

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

(ZTX1n)−1(ZTY1n)=β1+(ZTX1n)−1(ZTε1n)

ZTX1/nは定数に、ZTε1/nはZが誤差項εと無相関であれば0に確率収束する。

b1(IV)=(ZTX1)−1ZTY1

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

Y1=X1β1+ε1を代入すると

b1(IV)=(ZTX1)−1ZT(X1β1+ε1)=(ZTX1)−1ZTX1β1+(ZTX1)−1ZTε1=β1+(ZTX1)−1ZTε1

となる。確率極限をとるとplimn→∞ZTε1/n=0なので

plimn→∞b1(IV)=β1

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

nを大きくしていったとき、標本平均∑iXi/nは母平均μに一致する。これを確率極限といい、plim(∑iXi/n)=μのように表す。

 plim (∑iXi2n)= plim (xTxn)=μ2−σ2

（Var[X]=E[X2]−E[X]2なのでμ2が残っている)

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

 plim (∑iXiYin)= plim (xTyn)=μXμY−σXY

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

 plim (∑iXiεin)= plim (xTεn)=μXμε−σXε=0

となる

# 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 0x7f8ec6d18280>
../../_images/1a7b9810a1360ccbb27ce8f8a361c7f43e23a392d8afde0df947bd973fad28c5.png 
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"])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 18
     16 import seaborn as sns
     17 import pandas as pd
---> 18 df = pd.DataFrame(dict(y=y, x=x, z=z, u=u, e=e))
     19 sns.pairplot(df, height=1.5, y_vars=["y", "x"])

NameError: name 'e' is not defined

# 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
IV-2SLS Estimation Summary
Dep. Variable: y R-squared: 0.8714
Estimator: IV-2SLS Adj. R-squared: 0.8711
No. Observations: 500 F-statistic: 2402.9
Date: Sun, Jul 30 2023 P-value (F-stat) 0.0000
Time: 02:17:35 Distribution: chi2(1)
Cov. Estimator: robust
Parameter Estimates
Parameter Std. Err. T-stat P-value Lower CI Upper CI
x 0.7702 0.0157 49.020 0.0000 0.7394 0.8010


Endogenous: x
Instruments: z
Robust Covariance (Heteroskedastic)
Debiased: False
id: 0x7fe0f2d91f10

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

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

Contents