Breusch-Pagan and Koenker Score Diagnostics

Why Test for Variance Dependence?

A residual plot can suggest that error variability changes with the predictors, but a visual pattern does not quantify how much evidence there is against constant conditional variance. The Breusch–Pagan and Koenker tests turn that diagnostic question into a joint test: do the squared residuals systematically vary with a specified set of covariates? Their shared auxiliary-regression structure makes them convenient, while their different scaling reflects different assumptions about the error distribution.

These are diagnostics for the variance specification, not tests of whether the conditional mean is correctly specified. Let the mean model be \(y_i=x_i'\beta+u_i\), estimated by OLS, and let \(\hat u_i\) be its residuals. The covariate vector \(z_i\) contains the variables whose relationship with the error variance is being tested. It may, but need not, contain all the regressors in the mean model.

The exponential form ensures a positive variance and provides a convenient parametric model for deriving the score test. The null is the substantive target: the test asks whether the chosen variance covariates provide evidence of departures from constant variance. It does not test every possible form of heteroskedasticity.

The Auxiliary Regression

First estimate the conditional mean by OLS and calculate \(\hat u_i^2\). Then regress these squared residuals on an intercept and the chosen variance covariates:

Let \(R^2_{\text{aux}}\) denote the coefficient of determination from this auxiliary regression, and let \(q\) be the number of nonconstant covariates in \(z_i\). Under the null and suitable regularity conditions, the Koenker statistic is \(nR^2_{\text{aux}}\), which is compared with a chi-square distribution with \(q\) degrees of freedom. The conventional Breusch–Pagan score test uses normal-theory scaling; its familiar large-sample form is also \(nR^2_{\text{aux}}\). A finite-sample score expression makes the difference in scaling explicit, as derived below.

The intercept is included to allow the average squared residual to be nonzero. It is not counted among the \(q\) tested coefficients. The test is joint: the null is that all \(q\) slopes are zero, not that each slope is assessed by a separate test.

Deriving the Breusch–Pagan Score

The classical derivation assumes independent normal errors conditional on the covariates, with variance \(\sigma^2\exp(z_i'\gamma)\). Under \(H_0:\gamma=0\), the score for \(\gamma\), if the errors and \(\sigma^2\) were known, is

Because the intercept in \(z_i\) would describe the overall scale \(\sigma^2\), rather than a change in variance with covariates, the scale is a nuisance parameter. Removing that nuisance direction centers the tested covariates. Write \(Z_c\) for the matrix whose rows are \(z_i-\bar z\), where \(\bar z=n^{-1}\sum_i z_i\). Under normality, the efficient information for the centered variance slopes is \(\tfrac12 Z_c'Z_c\). Thus the score quadratic form is

This expression shows why the test is related to an auxiliary regression: \(r'P_{Z_c}r\) is the explained sum of squares when the centered, scaled squared errors are projected onto the centered variance covariates. In practice, the errors and \(\sigma^2\) are unknown, so OLS residuals and a residual-based scale estimate are used. Under the null, normality gives \(\operatorname{Var}(u_i^2/\sigma^2)=2\). Consequently the total variation in scaled squared errors is approximately \(2n\), and the score statistic is asymptotically equivalent to \(nR^2_{\text{aux}}\). With \(q\) tested slopes, its null reference distribution is \(\chi^2_q\).

Koenker’s Studentized Variant

The normal-theory score derivation uses the fact that a normal error has fourth moment \(3\sigma^4\), so the variance of its square is \(2\sigma^4\). If errors are not normal, the variation of squared errors need not have that value. Koenker’s studentized version adjusts the scaling using the observed variation in squared residuals. In its standard auxiliary-regression form, the statistic is \(nR^2_{\text{aux}}\), with a \(\chi^2_q\) reference distribution under the null.

The studentization makes this version less sensitive to nonnormality than the classical normal-theory score test. It does not remove the need for assumptions: the usual large-sample justification relies on independent observations, a suitable conditional-mean model, constant conditional variance under the null, and finite moments sufficient for the squared residuals and their variation to behave regularly. In particular, a heavy-tailed distribution with no finite fourth moment can undermine the usual chi-square approximation.

Both tests use the same chosen variance covariates and both are asymptotic diagnostics. Their statistics need not be identical in a finite sample because the classical BP score uses normal-theory scaling, whereas Koenker’s version studentizes using the empirical variation in squared residuals. The choice between them should reflect the error-distribution assumptions that are credible for the application, not a search for whichever p-value is smaller.

Worked Examples

Interpretation, Limits, and a Common Pitfall

A small p-value indicates that the observed auxiliary-regression association would be unusual under the constant-variance null and the assumptions used for the reference distribution. It is evidence against constant variance in the direction represented by the chosen \(z_i\). A large p-value is a failure to detect that pattern, not confirmation of homoskedasticity. The tests also do not measure the practical size of heteroskedasticity; inspect the residual scale and consider the substantive importance of any detected pattern.

A common mistake is to regress squared residuals on every available predictor, see a significant result, and conclude that the variance model is wrong in every possible way. The test only addresses its chosen covariate set and its asymptotic approximation. Another mistake is to regard the Koenker version as assumption-free because it is studentized. It relaxes the normality requirement, but still relies on independent observations and appropriate moment and model conditions.

Finally, a rejection does not invalidate OLS coefficient estimates by itself. Under the appropriate conditional-mean assumptions, OLS can still estimate the mean-model coefficients, while conventional variance formulas may be unsuitable under heteroskedasticity. Robust covariance estimation is a separate inferential response; the diagnostic tests do not substitute for deciding how uncertainty should be estimated.