Python中的OLS-Breusch-Pagan测试

Question

我使用statsmodels包来估计OLS回归。现在我要Breusch Pagan test。我将pysal包用于此测试，但此函数返回一个错误：

import statsmodels.api as sm
import pysal

model = sm.OLS(Y,X,missing = 'drop')
rs = model.fit()
pysal.spreg.diagnostics.breusch_pagan(rs)

返回错误：

AttributeError: 'OLSResults' object has no attribute 'u'

我该怎么办？

Answer 1

由于我不打算使用statsmodels库，所以创建了一个Python函数来执行Breusch Pagan测试。它使用SciKit learn中的多元线性回归。

import numpy as np
from sklearn.linear_model import LinearRegression
from scipy.stats import chisqprob


def breusch_pagan_test(x, y):
    '''
    Breusch-Pagan test for heteroskedasticity in a linear regression model:
    H_0 = No heteroskedasticity.
    H_1 = Heteroskedasticity is present.

    Inputs:
    x = a numpy.ndarray containing the predictor variables. Shape = (nSamples, nPredictors).
    y = a 1D numpy.ndarray containing the response variable. Shape = (nSamples, ).

    Outputs a list containing three elements:
    1. the Breusch-Pagan test statistic.
    2. the p-value for the test.
    3. the test result.
    '''

    if y.ndim != 1:
        raise SystemExit('Error: y has more than 1 dimension.')
    if x.shape[0] != y.shape[0]:
        raise SystemExit('Error: the number of samples differs between x and y.')
    else:
        n_samples = y.shape[0]

    # fit an OLS linear model to y using x:
    lm = LinearRegression()
    lm.fit(x, y)

    # calculate the squared errors:
    err = (y - lm.predict(x))**2

    # fit an auxiliary regression to the squared errors:
    # why?: to estimate the variance in err explained by x
    lm.fit(x, err)
    pred_err = lm.predict(x)
    del lm

    # calculate the coefficient of determination:
    ss_tot = sum((err - np.mean(err))**2)
    ss_res = sum((err - pred_err)**2)
    r2 = 1 - (ss_res / ss_tot)
    del err, pred_err, ss_res, ss_tot

    # calculate the Lagrange multiplier:
    LM = n_samples * r2
    del r2

    # calculate p-value. degrees of freedom = number of predictors.
    # this is equivalent to (p - 1) parameter restrictions in Wikipedia entry.
    pval = chisqprob(LM, x.shape[1])

    if pval < 0.01:
        test_result = 'Heteroskedasticity present at 99% CI.'
    elif pval < 0.05:
        test_result = 'Heteroskedasticity present at 95% CI.'
    else:
        test_result = 'No significant heteroskedasticity.'
    return [LM, pval, test_result]

Answer 2

问题是statsmodels的回归结果实例与pysal中的实例不兼容。

您可以使用statsmodels中的breushpagan，它接受OLS残差和异方差解释变量的候选者，因此它不依赖于模型的特定模型或实现。

文件： http://statsmodels.sourceforge.net/devel/generated/statsmodels.stats.diagnostic.het_breushpagan.html

这里有例子http://statsmodels.sourceforge.net/devel/examples/notebooks/generated/regression_diagnostics.html

我不知道是否有任何实质性的差异，在实施布雷乌什异教测试。

在statsmodels中，这个名字似乎拼错了。