regression-diagnostics 0.2.0

use super::{lm_heteroskedasticity_test, non_intercept_columns};
use crate::OlsFit;

/// Result of the Breusch-Pagan heteroskedasticity test (LM form).
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct BreuschPagan {
    /// Lagrange-multiplier statistic `n·R²_aux`.
    pub statistic: f64,
    /// Degrees of freedom (number of regressors in the auxiliary regression).
    pub df: usize,
    /// Upper-tail p-value under χ²(df).
    pub p_value: f64,
}

/// Breusch-Pagan test for heteroskedasticity.
///
/// Regresses the **squared residuals** on the model's original regressors and
/// forms the Lagrange-multiplier statistic `LM = n·R²_aux`, asymptotically
/// χ²(k) with `k` the number of non-intercept regressors. A small p-value is
/// evidence that residual variance depends on the predictors (non-constant
/// variance), violating the homoskedasticity assumption.
///
/// Breusch-Pagan tests specifically for variance that is a **linear** function of
/// the regressors; for a more general alternative (including nonlinear terms) see
/// [`white_test`](super::white_test).
pub fn breusch_pagan(fit: &OlsFit) -> BreuschPagan {
    let n = fit.n_observations();
    let cols = non_intercept_columns(fit);
    let x = fit.design_matrix();

    let regressors: Vec<Vec<f64>> = cols.iter().map(|&j| x.column(j).to_vec()).collect();

    let resid_sq: Vec<f64> = fit.residuals().iter().map(|e| e * e).collect();

    let (statistic, df, p_value) = lm_heteroskedasticity_test(n, &regressors, &resid_sq);
    BreuschPagan {
        statistic,
        df,
        p_value,
    }
}