自己造轮子——异方差的检验与WLS_基于Statsmodels

我知道用stata和R的更多，此处仅基于兴趣摸索，本文数据来自于朱顺全老师的教程，感谢本杰明·安德森教授！！！

检验异方差

bpg_test = het_breuschpagan(fit.resid, X) # Breusch-Pagan检验异方差
print("Breusch-Pagan Test:", bpg_test)

white_test = het_white(fit.resid, X) # White检验异方差
print("White Test:", white_test)

gq_test=het_goldfeldquandt(fit.resid, X) # Goldfeld-Quandt检验异方差
print("Goldfeld-Quandt Test:",gq_test)

使用statsmodels里面的三个内置函数来检验异方差（教材里手算Goldfeld-Quandt已经不需要了）

三种方法返回的数据一般都只需要关注第二位，即p值，小于0.05则认为存在异方差

WLS

wgt = 1 / fit.resid**2

wls_model = sm.WLS(y, X, weights=wgt)
wls_fit = wls_model.fit()

print(wls_fit.summary())

此处我是用残差平方的倒数作为权重（我记得李旭晖老师当时是这么教的，当然这地方有可以深挖的，所以关于异方差的处理，后面还要将继续深入）

完整代码

import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.stats.diagnostic import het_breuschpagan, het_white,het_goldfeldquandt
from patsy import dmatrices

#### 数据导入
data=pd.DataFrame()
data=pd.DataFrame(pd.read_excel('C:\\Users\\cyx94a\\Desktop\\Python\\data\\al8-2.xls'))
vars=['x','y']
df=data[vars]

#### 生成设计矩阵
y,X=dmatrices('y~x',data=df,return_type='dataframe')

#### 描述性统计和相关分析
print(data.describe())
print(data.corr())

#### OLS
model=sm.OLS(y,X)
fit=model.fit()
print(fit.summary())

#### 计算VIF
X_with_const = sm.add_constant(X) # 为矩阵X增加常数项

vif_values = [variance_inflation_factor(X_with_const.values, i) for i in range(X_with_const.shape[1])] # 计算每个特征的 VIF

print("VIF values:", vif_values) # 从列表的第二个数字及之后就是各个自变量的VIF值

####检验异方差
bpg_test = het_breuschpagan(fit.resid, X) # Breusch-Pagan检验异方差，返回列表的第二个数字是p值
print("Breusch-Pagan Test:", bpg_test)

white_test = het_white(fit.resid, X) # White检验异方差，返回列表的第二个数字是p值
print("White Test:", white_test)

gq_test=het_goldfeldquandt(fit.resid, X) # Goldfeld-Quandt检验异方差，返回列表的第二个数字是p值
print("Goldfeld-Quandt Test:",gq_test)

####用WLS修正异方差
wgt = 1 / fit.resid**2

wls_model = sm.WLS(y, X, weights=wgt)
wls_fit = wls_model.fit()

print(wls_fit.summary())