はじめに
この記事では、西谷慶彦・新谷元嗣・川口大司・奥井亮（2019）『計量経済学』有斐閣 の実証例を再現するPythonコードを紹介します。
使用するデータセットは有斐閣のサポートページにて公開されています。
URL：http://www.yuhikaku.co.jp/books/detail/9784641053854
※ この記事は随時追記していきます。
※ ある章で一度 import したライブラリは、その章の間は import を省略します。
第４章 線形単回帰モデルの推定と検定
実証例 4.1 単回帰分析 (pp.128-129)
均一分散
まずは線形単回帰分析のコードです。
書籍では、自由度を修正した不均一分散に頑健な標準誤差を使用していますが、
まずは均一分散を仮定したモデルの推定を行います
import pandas as pd
import numpy as np
# データの読み込み
df = pd.read_csv('ch04_wage.csv')
# 線形回帰に使用するライブラリの読み込み
from statsmodels.formula.api import ols
# モデルの作成（均一分散）
formula = 'wage ~ productivity' # formula = y ~ x と記述する
model = ols(formula, data=df).fit() # 均一分散
# 結果を表示
print(model.summary())
''' 結果
OLS Regression Results
==============================================================================
Dep. Variable: wage R-squared: 0.963
Model: OLS Adj. R-squared: 0.962
Method: Least Squares F-statistic: 501.5
Date: Tue, 14 Sep 2021 Prob (F-statistic): 4.04e-15
Time: 22:45:26 Log-Likelihood: -96.984
No. Observations: 21 AIC: 198.0
Df Residuals: 19 BIC: 200.1
Df Model: 1
Covariance Type: nonrobust
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Intercept 276.1296 87.611 3.152 0.005 92.759 459.501
productivity 0.5468 0.024 22.395 0.000 0.496 0.598
==============================================================================
Omnibus: 1.769 Durbin-Watson: 1.481
Prob(Omnibus): 0.413 Jarque-Bera (JB): 1.105
Skew: -0.239 Prob(JB): 0.575
Kurtosis: 1.983 Cond. No. 5.59e+04
==============================================================================
'''
結果を見ると、係数は教科書の実証例と等しいですが、標準誤差とそれを基に計算されるt値、p値、信頼区間が異なります。
これは上のコードでは、誤差項の均一分散を仮定して、推定しているためです。
不均一分散
次に教科書同様、誤差項の不均一分散を仮定して、自由度を修正した標準誤差を使用して推定を行います。
これには、上のコードの
model = ols(formula, data=df).fit()
という行の、fit()のオプションとしてcov_typeを指定する必要があります。
# 不均一分散の推定
model_hetero = ols.(formula, data=df).fit(cov_type='HC1', use_t = True)
print(model_hetero)
''' 結果
OLS Regression Results
==============================================================================
Dep. Variable: wage R-squared: 0.963
Model: OLS Adj. R-squared: 0.962
Method: Least Squares F-statistic: 714.1
Date: Tue, 14 Sep 2021 Prob (F-statistic): 1.55e-16
Time: 22:50:32 Log-Likelihood: -96.984
No. Observations: 21 AIC: 198.0
Df Residuals: 19 BIC: 200.1
Df Model: 1
Covariance Type: HC1
================================================================================
coef std err t P>|t| [0.025 0.975]
--------------------------------------------------------------------------------
Intercept 276.1296 71.256 3.875 0.001 126.990 425.269
productivity 0.5468 0.020 26.722 0.000 0.504 0.590
==============================================================================
Omnibus: 1.769 Durbin-Watson: 1.481
Prob(Omnibus): 0.413 Jarque-Bera (JB): 1.105
Skew: -0.239 Prob(JB): 0.575
Kurtosis: 1.983 Cond. No. 5.59e+04
==============================================================================
'''
上のコードを実行することで、教科書と同じ標準誤差が得られます。
ちなみに、頑健な（ロバスト）標準誤差にはいくつかの種類が存在しており、上のコードではHC1を指定しています。
以下、その他のオプションになります。
・HC0: White (1980)のロバスト標準誤差
・HC1: MacKinnon and White (1985)のロバスト標準誤差（１）
・HC2: MacKinnon and White (1985)のロバスト標準誤差（２）
・HC3: MacKinnon and White (1985)のロバスト標準誤差（３）
※詳しくは Angrist and Pischke(2008) "Mostly Harmless Econometrics" の Chapter8 を参照
_$\hat{\beta}1 = 1$ の検定
以下、手計算による方法と、関数を用いる方法で $\hat{\beta}_1 = 1$ の検定を行います。
# 手計算
(model_hetero.params['productivity'] -1)/ model_hetero.bse['productivity']
''' 結果
-22.14628233152748
'''
# 自動計算
hypothesis = 'productivity = 1'
model_hetero.t_test(hypothesis)
''' 結果
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
c0 0.5468 0.020 -22.146 0.000 0.504 0.590
==============================================================================
'''
区間推定
最後に、信頼係数95%の区間推定を行います。
# 自動計算
model_hetero.conf_int(alpha=0.05) # alpha=(1-信頼係数)
''' 結果（教科書では丸めた数値を使用しているため、数値が若干異なる）
0 1
Intercept 126.989943 425.269270
productivity 0.503990 0.589649
'''
# 手計算
upper = model_hetero.params['productivity'] + 1.96*model_hetero.bse['productivity']
lower = model_hetero.params['productivity'] - 1.96*model_hetero.bse['productivity']
'''
upper
0.5869271611134851
lower
0.5067120178869308
'''
以上が、第4章の実証例になります。
第５章 重回帰モデルの推定と検定
実証例 5.1 重回帰分析 (pp.151-152)
前章では説明変数が１つでしたが、複数に増やして重回帰分析を行います。
重回帰であっても、コードに大きな変化はなく、式の記述を少し変えるだけで推定を行えます。
具体的には、
formula = 'y ~ x1 + x2 + ...'
のように記述することで分析可能です。
import pandas as pd
import numpy as np
import statsmodels.formula.api import ols
# データの読み込み
df = pd.read_csv('youdou.csv')
# 使用する変数の作成
df['lny80'] = np.log(df['y80']) # 80年GDPの対数値
df['growth_rate'] = (np.log(df['y99']) - np.log(df['y80'])) / 19 *100 # 教科書で定義されている経済成長率
# 重回帰分析（信頼）
formula = 'growth_rate ~ trust80 + education80 + lny80'
MLR1 = ols(formula, data = df).fit()
print(MLR1.params)
''' 結果
Intercept 6.048849
trust80 0.020585
education80 2.612075
lny80 -2.383094
'''
# 重回帰分析（規範）
formula = 'growth_rate ~ norm80 + education80 + lny80'
MLR2 = ols(formula, data = df).fit()
print(MLR2.params)
''' 結果
Intercept 5.290872
norm80 0.338269
education80 4.387236
lny80 -1.991136
'''
実証例 5.2 FWL定理の確認 (p.154)
次に、FWL定理による係数の算出が上の結果と一致するかを確認します。
# 前準備
formula_trust = 'trust80 ~ education80 + lny80' # 信頼をその他の変数に回帰する式
formula_y = 'growth_rate ~ education80 + lny80' # 経済成長率を信頼以外の変数に回帰する式
model_trust = ols(formula_trust, data=df).fit()
model_y = ols(formula_y, data=df).fit()
# 残差の抽出
resid_trust = model_trust.resid # 残差
resid_y = model_y.resid # 残差
# 回帰
formula = 'resid_y ~ -1 + resid_trust' # -1は定数項なしにするため
model = ols(formula, data =df).fit()
print(model.params)
''' 結果
resid_trust 0.020585
'''
確かに、実証例 5.1 の結果と一致しています。
実証例 5.3 FWL定理の別表現の確認 (pp.155-156)
次に、FWL定理の別表現である、教科書の(5.7)式に基づいた係数の算出を行います。
formula = 'growth_rate ~ -1 + resid_trust'
model = ols(formula, data =df).fit()
print(model.params)
''' 結果
resid_trust 0.020585
'''
こちらも、実証例 5.1 の結果と一致します。
※ただし、標準誤差が正しく計算されません。
実証例 5.4 標準誤差（p.163）
今までの実証例では、係数のみを表示していましたが、単回帰分析と同様の方法で標準誤差も簡単に確認できます。
なお、標準誤差は単回帰分析の時と同様に自由度修正済みロバスト標準誤差を用いています。
# 信頼についての回帰
formula = 'growth_rate ~ trust80 + education80 + lny80'
model1 = ols(formula, data = df).fit(cov_type='HC1') # 自由度修正済みロバスト標準誤差を使用
print(model1.summary())
''' 結果
OLS Regression Results
==============================================================================
Dep. Variable: growth_rate R-squared: 0.562
Model: OLS Adj. R-squared: 0.531
Method: Least Squares F-statistic: 20.21
Date: Wed, 22 Sep 2021 Prob (F-statistic): 2.53e-08
Time: 23:04:45 Log-Likelihood: -9.1361
No. Observations: 47 AIC: 26.27
Df Residuals: 43 BIC: 33.67
Df Model: 3
Covariance Type: HC1
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 6.0488 0.426 14.185 0.000 5.213 6.885
trust80 0.0206 0.076 0.272 0.786 -0.128 0.169
education80 2.6121 2.709 0.964 0.335 -2.697 7.921
lny80 -2.3831 0.491 -4.849 0.000 -3.346 -1.420
==============================================================================
Omnibus: 11.045 Durbin-Watson: 1.895
Prob(Omnibus): 0.004 Jarque-Bera (JB): 11.189
Skew: 0.945 Prob(JB): 0.00372
Kurtosis: 4.464 Cond. No. 85.9
==============================================================================
'''
# 規範についての回帰
formula = 'growth_rate ~ norm80 + education80 + lny80'
model2 = ols(formula, data = df).fit(cov_type='HC1')
print(model2.summary())
''' 結果
OLS Regression Results
==============================================================================
Dep. Variable: growth_rate R-squared: 0.639
Model: OLS Adj. R-squared: 0.614
Method: Least Squares F-statistic: 41.04
Date: Wed, 22 Sep 2021 Prob (F-statistic): 1.11e-12
Time: 23:06:37 Log-Likelihood: -4.5784
No. Observations: 47 AIC: 17.16
Df Residuals: 43 BIC: 24.56
Df Model: 3
Covariance Type: HC1
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 5.2909 0.668 7.918 0.000 3.981 6.601
norm80 0.3383 0.137 2.469 0.014 0.070 0.607
education80 4.3872 1.961 2.237 0.025 0.544 8.231
lny80 -1.9911 0.575 -3.465 0.001 -3.117 -0.865
==============================================================================
Omnibus: 2.553 Durbin-Watson: 1.934
Prob(Omnibus): 0.279 Jarque-Bera (JB): 2.053
Skew: 0.512 Prob(JB): 0.358
Kurtosis: 2.987 Cond. No. 84.9
==============================================================================
'''
おまけ：stargazer
Rにおける、計量経済分析の結果を見やすく表示する方法として、stargazer というパッケージが広く使用されています。
この stargazer ですが、実はPythonでも使えます。
ここでは本当に簡単に、コードと実際の結果について紹介します。
細かい仕様については、マニュアルを参照ください（機会があれば、新たな記事で詳述していきたいと思います）。
# 実証例 5.4 の結果を stargazerで表示
# ライブラリの読み込み
from stargazer.stargazer import Stargazer
model1_stargazer = Stargazer([model1])
model1_stargazer
結果
| Dependent variable:growth_rate | ||
| (1) | (2) | |
| Intercept | 6.049*** | 5.291*** |
| (0.426) | (0.668) | |
| education80 | 2.612 | 4.387** |
| (2.709) | (1.961) | |
| lny80 | -2.383*** | -1.991*** |
| (0.491) | (0.575) | |
| norm80 | 0.338** | |
| (0.137) | ||
| trust80 | 0.021 | |
| (0.076) | ||
| Observations | 47 | 47 |
| R2 | 0.562 | 0.639 |
| Adjusted R2 | 0.531 | 0.614 |
| Residual Std. Error | 0.307 (df=43) | 0.279 (df=43) |
| F Statistic | 20.209*** (df=3; 43) | 41.038*** (df=3; 43) |
| Note: | *p<0.1; **p<0.05; ***p<0.01 | |
非常にコンパクトに見やすくまとまります。特に、論文などで複数のモデルを作成し、比較などを行う際には重宝されるかと思います。
R同様に、LaTeXやhtml形式での表示もできます。
# LaTeX 形式で表示
print(model1_stargazer.render_latex())
実証例 5.5 多項式モデル（pp.170-171）
次に、説明変数に二乗項を入れる方法について紹介します。
I() という表現を使うと容易に二乗項を含めることが可能です。
カッコの中に計算式を記述します（ 二乗項なら I(y80*2) ）。
formula = 'growth_rate ~ y80 + I(y80**2)'
model = ols(formula, data = df).fit(cov_type='HC1')
print(model.summary())
''' 結果
OLS Regression Results
==============================================================================
Dep. Variable: growth_rate R-squared: 0.550
Model: OLS Adj. R-squared: 0.530
Method: Least Squares F-statistic: 27.39
Date: Thu, 23 Sep 2021 Prob (F-statistic): 1.88e-08
Time: 16:44:37 Log-Likelihood: -9.7476
No. Observations: 47 AIC: 25.50
Df Residuals: 44 BIC: 31.05
Df Model: 2
Covariance Type: HC1
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
Intercept 6.5187 1.385 4.705 0.000 3.803 9.234
y80 -1.2262 0.708 -1.732 0.083 -2.614 0.161
I(y80 ** 2) 0.0894 0.089 1.008 0.313 -0.084 0.263
==============================================================================
Omnibus: 6.142 Durbin-Watson: 1.908
Prob(Omnibus): 0.046 Jarque-Bera (JB): 4.966
Skew: 0.726 Prob(JB): 0.0835
Kurtosis: 3.654 Cond. No. 614.
==============================================================================
'''
または、データフレームの中に説明変数を二乗した新たな変数を作成することでも、推定可能です。
# 変数の作成
df['y80_sq'] = df['y80']**2
df.head()
# 推定
formula = 'growth_rate ~ y80 + y80_sq'
model = ols(formula, data = df).fit(cov_type='HC1')
print(model.summary())
''' 結果
OLS Regression Results
==============================================================================
Dep. Variable: growth_rate R-squared: 0.550
Model: OLS Adj. R-squared: 0.530
Method: Least Squares F-statistic: 27.39
Date: Thu, 23 Sep 2021 Prob (F-statistic): 1.88e-08
Time: 16:47:57 Log-Likelihood: -9.7476
No. Observations: 47 AIC: 25.50
Df Residuals: 44 BIC: 31.05
Df Model: 2
Covariance Type: HC1
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 6.5187 1.385 4.705 0.000 3.803 9.234
y80 -1.2262 0.708 -1.732 0.083 -2.614 0.161
y80_sq 0.0894 0.089 1.008 0.313 -0.084 0.263
==============================================================================
Omnibus: 6.142 Durbin-Watson: 1.908
Prob(Omnibus): 0.046 Jarque-Bera (JB): 4.966
Skew: 0.726 Prob(JB): 0.0835
Kurtosis: 3.654 Cond. No. 614.
==============================================================================
'''
実証例 5.6 交互作用モデル（１）（pp.173-174）
次に、二つの説明変数を掛け合わせた交差項をモデルに含める方法を紹介します。
まずは、教科書で定義されているダミー変数を作成します。
# ダミー変数の作成
df['D'] = (df['did']>0.4).astype(int)
(df['did']>0.4)によって、各行に対して TrueかFalse が返されます。
それを整数型に変換することで0 or 1 を取るダミー変数を作成しています。
次に、多項式モデルの時と同様に I() を使用して、交差項を含むモデルの推定を行います。
formula = 'growth_rate ~ D + lny80 + I(D*lny80)'
model = ols(formula, data=df).fit(cov_type='HC1')
print(model.params)
''' 結果
Intercept 5.749055
D -0.175513
lny80 -1.911200
I(D * lny80) 0.064410
'''
また、同様の推定値はダミー変数を使用せずとも、
ダミー変数が１になる標本と、０になる標本で分けて推定を行っても得ることができます。
ols()の中のオプションとして、subset()を指定することで行います。
# ダミー変数が０となる標本
formula = 'growth_rate ~ lny80'
model = ols(formula, data=df, subset=(df['D']==0)).fit(cov_type='HC1')
print(model.params)
''' 結果
Intercept 5.749055
lny80 -1.911200
'''
# ダミー変数が１となる標本
formula = 'growth_rate ~ lny80'
model = ols(formula, data=df, subset=(df['D']==1)).fit(cov_type='HC1')
print(model.params)
''' 結果
Intercept 5.573541
lny80 -1.846790
'''
実証例 5.7 交互作用モデル（２）（p.175）
次は二つのダミー変数を用いて、それらの交差項をモデルに加えます。
コードとしては、実証例 5.6 と同様です。
# ダミー変数の作成
df['D1'] = (df['did']>0.4).astype(int)
df['D2'] = 1*(df['lny80']>1.4) # この方法でもダミー変数を取得できます
formula = 'growth_rate ~ D1 + D2 + I(D1*D2)'
model = ols(formula, data = df).fit(cov_type='HC1')
print(model.params)
''' 結果
Intercept 3.454739
D1 -0.233289
D2 -0.580030
I(D1 * D2) 0.047252
'''
実証例 5.8 F検定（p.178）
5章の最後として、結合仮説を検定するためにF検定を行うコードを紹介します。
例として、多項式モデルの y80 とその二乗項がともに０であるという結合仮説を検定します。
これは実証例 4.1 とほぼ同様のコードで行えます。
# モデルの推定
formula = 'growth_rate ~ y80 + I(y80*y80)'
model = ols(formula, data=df).fit(cov_type='HC1')
# F検定
hypotheses = 'y80=0, I(y80 * y80)=0'
# I(y80 * y80)はスペースの全角半角の関係でsummary()の表からコピーしないとエラーが発生することがあります。
f_test = model.f_test(hypotheses)
print(f_test.summary)
''' 結果
<F test: F=array([[27.38610602]]), p=1.8792896482418273e-08, df_denom=44, df_num=2>>
結果は左から、
F統計量、p値、分母の自由度、分子の自由度（仮説の数）
となっています。
Reference
・西谷慶彦・新谷元嗣・川口大司・奥井亮（2019）『計量経済学』有斐閣
・Angrist, J. D. and J.-S. Pischke (2008) "Mostly Harmless Econometrics: An empiricist's Companion", Princeton University Press.