BP 异方差检验

Breusch-Pagan 异方差检验。

核心代码

    @staticmethod
    def breusch_pagan_test(
        df: pd.DataFrame,
        y_var: str,
        x_vars: List[str],
        decimals: int = 4
    ) -> Dict[str, Any]:
        """
        Breusch-Pagan检验：检验异方差性

        参数:
            df: 数据框
            y_var: 因变量名
            x_vars: 自变量列表
            decimals: 小数位数

        返回:
            包含检验结果的字典
        """
        # 准备数据
        cols = [y_var] + x_vars
        df_clean = df[cols].dropna()

        if len(df_clean) < len(x_vars) + 2:
            raise ValueError(_('样本量不足以进行Breusch-Pagan检验'))

        # 1. 估计OLS模型
        formula = build_patsy_formula(y_var, x_vars)
        model = smf.ols(formula, data=df_clean)
        result = model.fit()

        # 2. 获取残差
        residuals = result.resid

        # 3. 计算残差平方
        resid_squared = residuals ** 2

        # 4. 用残差平方对自变量回归
        df_aux = df_clean.copy()
        df_aux['resid_sq'] = resid_squared

        aux_formula = build_patsy_formula('resid_sq', x_vars)
        aux_model = smf.ols(aux_formula, data=df_aux)
        aux_result = aux_model.fit()

        # 5. 计算LM统计量 = n * R²
        n = len(df_clean)
        r_squared = aux_result.rsquared
        lm_stat = n * r_squared

        # 6. 自由度 = 自变量个数
        df_test = len(x_vars)

        # 7. 计算p值（卡方分布）
        p_value = 1 - chi2.cdf(lm_stat, df_test)

        # 8. 判断结果
        if p_value < 0.01:
            conclusion = _('强烈拒绝原假设，存在显著异方差')
            stars = "***"
        elif p_value < 0.05:
            conclusion = _('拒绝原假设，存在异方差')
            stars = "**"
        elif p_value < 0.1:
            conclusion = _('弱拒绝原假设，可能存在异方差')
            stars = "*"
        else:
            conclusion = _('不能拒绝原假设，不存在显著异方差')
            stars = ""

        return {
            'test_name': 'Breusch-Pagan Test for Heteroskedasticity',
            'null_hypothesis': _('误差项具有同方差性（不存在异方差）'),
            'alternative_hypothesis': _('误差项具有异方差性'),
            'lm_statistic': round(lm_stat, decimals),
            'chi2_statistic': round(lm_stat, decimals),
            'df': df_test,
            'p_value': round(p_value, decimals),
            'significance': stars,
            'conclusion': conclusion,
            'n_obs': n,
            'aux_r_squared': round(r_squared, decimals)
        }