操作変数法の流れをまとめてみた（単回帰分析からの2段階最小二乗法）

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

https://amzn.to/3MSTwAV

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

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

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

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

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

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

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

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

# 実装例
コードはcolabでも公開しています

https://colab.research.google.com/drive/1nmSvTQSVLHgnOitY1JQWRwoeZMdCxiUK#scrollTo=Oe7D_1aFbsPa

1. ライブラリーのインポート
    ```python
    import numpy as np
    import scipy.stats as st
    import statsmodels.api as sm
    
    from statsmodels.sandbox.regression.gmm import IV2SLS, IVGMM
    ```
2. データの準備
    今回はMrozデータセットを使います。このデータセットはMroz Family Income and Education Dataとして知られており、教育年数や収入に関する情報を含んでいます。
    ```python
    #%% 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|母親の教育年数|
    
3. 単回帰分析(OLS)を実施
    ```python
    #%% 収入を教育年数で説明する単回帰モデル
    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²の低さや他の変数の影響を考慮すると、教育年数以外にも収入に影響を与える要因もありそうです( -> 内生性の問題を疑う)
    ![スクリーンショット 2024-09-21 11.23.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/846404/552f4091-c7f8-0520-a695-eeebb0d2fde3.png)
    

4. Breusch-Pagan検定
    単回帰分析の誤差項に異分散性があるかを確認します。
    異分散性の有無によって、どの操作変数法を使うべきかが変わります
    ```python
    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で問題ない可能性があります。

5. 2段階最小二乗法(2SLS)を行う
    OLSでは本人の教育年数を用いて、労働収入を推定しました。しかし本人の教育年数は、両親の教育年数で推定できるような気がします

    ```python
    #%% 父母の教育年数を操作変数として使う2SLS
    z = mroz.data[['educwf', 'educwm']]
    Z = sm.add_constant(z)
    results_iv = IV2SLS(y, X, instrument=Z).fit()
    print(results_iv.summary())
    ```

    ![スクリーンショット 2024-09-21 11.35.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/846404/ec5d59f7-1e56-1651-9485-e09372ddde9e.png)

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

    ```python
    # 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
    ```
    
    ```python
    # 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を使用しましょう!!