In [ ]:
% matplotlib inline
import statsmodels.api as sm
import statsmodels.discrete.discrete_model as smds
import numpy as np
import pandas as pd
内生性の問題を無視して推計すればその推定量は一致生を欠く。
例)
$ y_{1i}^* = a_0 + a_1y_{2i}+b_1x_{1i}+e_{1i}$
$ y_{2i} = b_0 + b_1x_{1i}+e_{2i}$
は観測できる変数。
さらに2本の連立方程式モデルで、双方の式に相手方の被説明変数が明示的に説明変数に含まれるというケースもある。
例)
$ y_{1i}^ = a_0 + a1y{2i}^+b1x{1i}+e{1i}$
$ y{2i}^ = b_0 + a2y{1i}^*+b2x{1i}+e_{2i}$
Maddala[1983]のモデル1〜モデル6
①被説明変数がいずれも直接観察される連続変数であれば、二段階最小自乗法や三段階最小自乗法により推計を行うことができる。
それ以外の場合、
②連続変数と切断された変数。
③連続変数と二値的変数。
④いずれも切断された変数。
⑤切断された変数と二値的変数。
⑥いずれも二値的変数。
は二段階推定などで内生性の問題を解く必要がある。
例)
$ y_{1i}^* = b_1x_{1i} + b_2y_{2i}^*+u_{1i}$
$ y_{1i}=y_{1i}^* $
$ y_{2i} = cz_i + v_i $
$y_{2i}$が内生変数か外生変数かを確かめたい。$u_i$と$v_i$が相関していなければ$y_{2i}$は外生変数である。
それを確かめるには、$u_i = \theta v_i + e_i$を考えなければならない。θ=０かどうかを調べれば良い。もしviが観測できれば、
$ y{1i}^ = b1x{1i} + b2y{2i}^+\theta vi + e{i}$
をOLSで推定し、t=0かどうかをt検定で検定すればよい。ただし、それはこの問題設定ではできないので、
$ y_{2i} = cz_i + v_i $
をOLSで推定し、cのOLS推定量$ \hat{c}$と$\hat{vi}=y{2i}-\hat{c}zi$を得て、
$ y{1i} = bix{1i}+b2y{2i}+\theta \hat{v_i} + e_i + \theta(v_i - \hat{v_i})$
について考える。θ=0という帰無仮説のもとでは第五項を無視してよいので、上式をOLSで推定し、θ=0という仮説をt検定で行える。
$ y_{1i}^* = b_1x_{1i} + b_2y_{2i}^*+u_{1i}$
$ y_{1i}=\begin{cases}
1 ~~ if y_{1i}^* > 0\\
0 ~~ otherwise
\end{cases} $
$ y_{2i} = cz_i + v_i $
$ y_{1i}^* = b_1x_{1i} + b_2y_{2i}^*+u_{1i}$
$ y_{1i}=\begin{cases}
y_{1i}^* ~~ if y_{1i}^* > 0\\
0 ~~ otherwise
\end{cases} $
$ y_{2i} = cz_i + v_i $
を考える。
$ y_{2i} = cz_i + v_i $
をOLSで推定し、cのOLS推定量$ \hat{c}$と$\hat{v_i}=y_{2i}-\hat{c}z_i$を得て、これを説明変数として追加して、
$ y_{1i} = b_1x_{1i}+b_2y_{2i}+b_3 \hat{v_i} + u_i $
をトービット・モデルで推計し、b_3=0という仮説をt検定で行う。
Mroz[1987]のデータセットを用いて、上記の理論を女性の労働市場参加・不参加と労働時間について適用する。
$ inlf_i^* = a_1 + c_1 infant1_i+c_2infant2_i+c_3age_i+c_4nwifeinc_i+c_5exper_i + c_6expersq_i + c_7educ_i + w_i $
$ inlf_i = \begin{cases}
1 ~~ if inlf_i^* > 0\\
0 ~~ otherwise
\end{cases} $
ここで$w_i$はN(0, 1)に従うとする。女性の就業決定に関し教育歴(educ)が内生変数か外生変数かを検定する。これは女性が将来の就業(inlf)を考慮して教育歴を決定しているという仮説の検証で用いられる検定方法である。教育歴educは母親と父親の学歴に依存する。
$educ_i = a_0 + b_1motheduc_i + b_2 fatheduc_i + u_i $
これをOLSで推定し、その残差を用いて
$ inlf_i^* = a_1 + c_1 infant1_i+c_2infant2_i+c_3age_i+c_4nwifeinc_i+c_5exper_i + c_6expersq_i + c_7educ_i + w_i $
に追加し、その残差(educres)が有意かどうか下記の式において検証する。
$ inlf_i^* = a_1 + c_1 infant1_i+c_2infant2_i+c_3age_i+c_4nwifeinc_i+c_5exper_i + c_6expersq_i + c_7educ_i + c_8educres_i + w_i $
In [7]:
mroz = pd.read_csv("mroz.csv")
In [8]:
X = mroz['inlf']
In [9]:
X = np.column_stack((mroz['inlf'], mroz['inlf'] ))
In [10]:
infant1 = lambda x : 1 if x == 1 else 0
infant2 = lambda x : 1 if x >= 2 else 0
child1 = lambda x : 1 if x <= 2 else 0
child2 = lambda x : 1 if x >=3 else 0
infant1 = pd.DataFrame([infant1(x) for x in mroz['kidslt6']],columns=['infant1'], index=[i for i in range(0, 753) ])
infant2 = pd.DataFrame([infant2(x) for x in mroz['kidslt6']],columns=['infant2'], index=[i for i in range(0, 753)])
child1 = pd.DataFrame([child1(x) for x in mroz['kidsge6']],columns=['child1'], index=[i for i in range(0, 753)])
child2 = pd.DataFrame([child2(x) for x in mroz['kidsge6']],columns=['child2'], index=[i for i in range(0, 753)])
In [11]:
add = pd.concat([infant1, infant2, child1, child2], axis=1)
In [13]:
mroz = pd.concat([mroz, add], axis=1)
In [16]:
X = np.column_stack((mroz['motheduc'], mroz['fatheduc']))
In [17]:
X = sm.add_constant(X)
y = mroz['educ']
In [18]:
model = sm.OLS(y, X)
results = model.fit()
print(results.summary())
In [20]:
educres = results.resid_pearson
In [22]:
X = np.column_stack((mroz['infant1'], mroz['infant2'], mroz['age'], mroz['nwifeinc'], mroz['exper'], mroz['expersq'], mroz['educ'], educres))
In [23]:
X = sm.add_constant(X)
In [24]:
inlf= lambda x : 1 if x > 0 else 0
y = [inlf(x) for x in mroz['inlf']]
In [25]:
model = smds.Probit(y, X)
results = model.fit()
print(results.summary())
結果を見ると、educresの係数は-0.0364と負であるが、漸近的t値は-0.339、p値は0.735である。
よって、教育歴が女性就業の外生変数であるという帰無仮説は、つまりeducresの係数が0であるという仮説は、棄却されない。
教育歴 = 外生変数 ＋ educres
女性就業を教育歴を含む変数で回帰した時の残差 = θeducres ＋ 残差
女性就業 = 学生変数 + 教育歴 + θeducres ＋ 残差
トービット・モデルの場合も同様に検定を行うことができる。
$(hours/1000)_i=
\begin{cases}
(hours/1000)_i & if (hours/1000)_i > 1\\
0 & otherwise
\end{cases}
$
・プロビット・モデルに連続説明変数が内戦変数として含まれるとする。
$ y_{1i}^* = b_1x_{1i} + b_2y_{2i}^*+u_{1i}$
$ y_{1i}=\begin{cases}
1 ~~ if y_{1i}^* > 0\\
0 ~~ otherwise
\end{cases} $
$ y_{2i} = cz_i + v_i $
$ (u_i, v_i) \sim NIID ((0.0), (1, \sigma^2), \rho)$
y_2iとu_iは相関を持つ。＝プロビット推計は一致性をもたない。
線形射影より
$u_i(Cov(u_i, v_i)/\sigma^2)\cdot v_i + e_i = \theta _i + e_i $
ここで$e_i$は$x_i、v_i、y_{2i}$と独立。
In [ ]: