はじめに

「Pythonによる計量経済学入門」を読みました。
そこで操作変数法が紹介されていたので、書籍をベースに自分で調べた内容をメモします。

流れ

OLS推定: まず、OLSを用いてモデルを推定します
Breusch-Pagan検定: OLSの誤差項に異分散性があるかを確認します。異分散性が認められればGMM、認められなければ2SLSを使います
IV推定: 内生性が疑われる場合は、Breusch-Pagan検定の結果を参考に2SLSまたはGMMを用いて推定します
Hausman検定: OLSとIV法の推定結果を比較し、内生性があるかどうかを判断します

そもそも操作変数(IV)法とは

操作変数法とは、モデルに内生性の問題がある場合に、それを解決するための手法です。

ここでの内生性とは、説明変数が誤差項と相関していることを指します。この問題があると、通常の最小二乗法（OLS）ではバイアスのある推定量が得られてしまいます。

内生性の問題

計量経済学では、モデルにおいて説明変数（独立変数）が誤差項と無相関であることを前提に、OLS推定を行います。しかし、実際のデータでは説明変数と誤差項が相関していることがあります。これが内生性の問題です。

内生性が生じる理由には、以下のようなものがあります。

省略変数バイアス: モデルに含めるべき重要な変数が省略されていると、その影響が誤差項に含まれ、説明変数と相関が生じます
逆因果関係: 説明変数が従属変数に影響を与えるのではなく、その逆の因果関係がある場合です
同時性バイアス: 説明変数と従属変数が同時に決定されると、両者に相互の影響が生じます
測定誤差: 説明変数が誤って測定されている場合、その誤差が誤差項に含まれるため、相関が発生します

操作変数法

こうした内生性の問題が疑われる場合、操作変数法を使います。
操作変数法には、代表的な方法が2つあります

モデル	特徴
2段階最小二乗法(2SLS)	シンプル。誤差項に異分散性がない場合に使える
一般化モーメント法(GMM)	操作変数法を一般化して、異分散性がある場合にも使えるようにしたもの

実装例

コードはcolabでも公開しています

ライブラリーのインポート

import numpy as np
import scipy.stats as st
import statsmodels.api as sm

from statsmodels.sandbox.regression.gmm import IV2SLS, IVGMM

データの準備
今回はMrozデータセットを使います。このデータセットはMroz Family Income and Education Dataとして知られており、教育年数や収入に関する情報を含んでいます。
```
#%% RdatasetsからMrozの読み込み
mroz = sm.datasets.get_rdataset('Mroz', 'Ecdat')
mroz.data = mroz.data[mroz.data['hearnw'] > 0]

mroz.data.head(5)
```
なお今回使用する列は、以下の通りです。

列名概要

hearnw 労働収入

educw 自分自身の教育年数

educwf 父親の教育年数

educwm 母親の教育年数

単回帰分析(OLS)を実施

#%% 収入を教育年数で説明する単回帰モデル
y = np.log(mroz.data['hearnw']) #労働収入
x = mroz.data[['educw']] #自分自身の教育年数
X = sm.add_constant(x)
results_ols = sm.OLS(y, X).fit(use_t=False)

print(results_ols.summary())

▼統計的に有意な関係がある一方で、R²の低さや他の変数の影響を考慮すると、教育年数以外にも収入に影響を与える要因もありそうです( -> 内生性の問題を疑う)

Breusch-Pagan検定
単回帰分析の誤差項に異分散性があるかを確認します。
異分散性の有無によって、どの操作変数法を使うべきかが変わります

from statsmodels.stats.diagnostic import het_breuschpagan

# OLS結果を使ってBreusch-Pagan検定を行う
bp_test = het_breuschpagan(results_ols.resid, results_ols.model.exog)
bp_stat = bp_test[0]
bp_p_value = bp_test[1]

print(f'Breusch-Pagan 検定統計量: {bp_stat}')
print(f'p値: {bp_p_value}')

# p値が小さい場合は異分散性が疑われるため、GMMを使うことを検討
if bp_p_value < 0.05:
    print("異分散性があるため、GMMを検討すべきです。")
else:
    print("異分散性が見られないため、2SLSで問題ない可能性があります。")

▼今回のデータだと、異分散性がないので2SLSで良さそうです

Breusch-Pagan 検定統計量: 0.30858650702479684
p値: 0.5785488590370306
異分散性が見られないため、2SLSで問題ない可能性があります。

2段階最小二乗法(2SLS)を行う
OLSでは本人の教育年数を用いて、労働収入を推定しました。しかし本人の教育年数は、両親の教育年数で推定できるような気がします
```
#%% 父母の教育年数を操作変数として使う2SLS
z = mroz.data[['educwf', 'educwm']]
Z = sm.add_constant(z)
results_iv = IV2SLS(y, X, instrument=Z).fit()
print(results_iv.summary())
```

Hausman検定
OLSと2SLSの推定結果を比較し、内生性があるかどうかを判断します

# 2段階最小二乗法（2SLS）とOLSの差を比較するための関数を定義
def hausman_test(ols_results, iv_results):
    """
    ols_results: OLS推定結果
    iv_results:  2SLS推定結果
    """
    # OLSとIVの推定量の差を計算
    b_ols = ols_results.params
    b_iv = iv_results.params
    diff = b_iv - b_ols
    
    # OLSの推定量の共分散行列
    cov_ols = ols_results.cov_params()
    # IVの推定量の共分散行列
    cov_iv = iv_results.cov_params()
    
    # 差の共分散行列を計算
    cov_diff = cov_iv - cov_ols
    
    # Hausman検定の統計量を計算
    stat = np.dot(np.dot(diff.T, np.linalg.inv(cov_diff)), diff)
    
    # 自由度
    df = len(b_ols) - 1
    
    # p値を計算
    p_value = st.chi2.sf(stat, df)
    
    return stat, p_value

# OLSと2SLSの結果を使ってHausman検定を実行
hausman_stat, p_value = hausman_test(results_ols, results_iv)

print(f'Hausman検定統計量: {hausman_stat}')
print(f'p値: {p_value}')

# p値に基づく判定
if p_value < 0.05:
    print("帰無仮説を棄却: OLS推定量はバイアスがあり、2SLSを使用すべきです。")
else:
    print("帰無仮説を採択: OLS推定量にバイアスはなく、OLSを使用することが適切です。")

▼今回の場合、2SLSを使用した方が良いようです

Hausman検定統計量: 4.087958454157072
p値: 0.04318977974246603
帰無仮説を棄却: OLS推定量はバイアスがあり、2SLSを使用すべきです。

最後に

今回はOLSを行ったあと、内生性の問題を疑い、操作変数法を試してみました。
Hausman検定によれば、OLS推定量にはバイアスがあるため2SLSを使用することが適切です。

R^2値に注目すればOLSの方が高いものの、OLSにはバイアスが含まれるため、きちんと2SLSを使用しましょう!!