不均一分散一致標準誤差(HC Standard Error)によるt検定

回帰分析のロバストなt検定
1. ロバストな標準誤差
2. ロバストなt値
Pythonで不均一分散にロバストなt検定
1. Statsmodels
2. sklearn
参考

回帰分析のロバストなt検定

線形回帰分析では、誤差項が均一分散であることが仮定されている。
誤差項が不均一分散であるとき、標準誤差に誤りが生じる。

誤差項が不均一分散の線形回帰分析では、誤った標準誤差から、回帰係数のt値を求めることになる。

ロバストな標準誤差

不均一分散一致標準誤差(heteroskedasticity consistent standard error)
$$
\begin{aligned} v_{\mathrm{HCE}}\left[\widehat{\beta}_{\mathrm{OLS}}\right] &=\frac{1}{n}\left(\frac{1}{n} \sum_{i} X_{i} X_{i}^{\prime}\right)^{-1}\left(\frac{1}{n} \sum_{i} X_{i} X_{i}^{\prime} \widehat{u}_{i}^{2}\right)\left(\frac{1}{n} \sum_{i} X_{i} X_{i}^{\prime}\right)^{-1} \\ &=\left(\mathrm{X}^{\prime} \mathrm{X}\right)^{-1}\left(\mathrm{X}^{\prime} \operatorname{diag}\left(\widehat{u}_{1}^{2}, \ldots, \widehat{u}_{n}^{2}\right) \mathrm{X}\right)\left(\mathrm{X}^{\prime} \mathrm{X}\right)^{-1} \end{aligned}
$$

引用:https://en.wikipedia.org/wiki/Heteroscedasticity-consistent_standard_errors

memo

$$
HC_{se} = \sqrt{\frac{\sum{\{(X_i – \bar{X}) ^2\hat{u_i}^2\}}}{\{\sum{(X_i – \bar{X})^2}\}^2}}
$$

参考:44の例題で学ぶ計量経済学

ロバストなt値

$$
t_{hce} = \frac{\hat{\beta}}{HC_{se}}
$$

Pythonで不均一分散にロバストなt検定

上の式から、Pythonでロバストなt検定を実行する。

通常なt検定は、こちらから

ボストン住宅価格のデータセットを使う。
このデータにおける、重回帰分析モデルの誤差項は、不均一分散であった。

BP検定で、不均一分散かどうか検定した。
heteroskedasticityの検定-breusch-paganbp検定

Statsmodelsを使って、不均一分散にロバストなt検定を行います。

sklearn+自作でも試してるが、Statsmodelsと一致しないので、まだ勉強中です。

Statsmodels

#import
from sklearn.datasets import load_boston
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
#data_loading and data_setting
boston = load_boston()
data_X = boston.data
data_y = boston.target

results = sm.OLS(data_y, sm.add_constant(data_X)).fit(cov_type='HC0')

pd.DataFrame({'t_hc_se':results.tvalues[1:]}, index=boston.feature_names).T

sklearn

#import
from sklearn.datasets import load_boston
import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression

#data_loading and data_setting
boston = load_boston()
data_X = boston.data
data_y = boston.target

#回帰分析
lr = LinearRegression()
lr.fit(data_X, data_y)
coef_value = lr.coef_
y_pred = lr.predict(data_X)

X = data_X

e = (data_y - y_pred)
e = np.diag((e.T**2).reshape(-1))

a = np.linalg.inv(np.dot(X.T, X))
b = np.dot(np.dot(X.T, e), X)

hc_se = np.sqrt(np.diagonal(np.dot(np.dot(a, b), a)))
t_hc_se = coef_value /  hc_se

pd.DataFrame({'t_hc_se': t_hc_se}, index=boston.feature_names).T

参考