6. Residual diagnostics and heteroskedasticity

Level 7: Rigorous Statistical Diagnostics

We are rarely sure of having specified the model correctly, and we have already seen that R squared can be artificially high or low regardless of the specification. Hence the need for diagnostic tests. The most critical assumptions are those in UR.3 and UR.4: constant variance, no autocorrelation and normality.

First, a fundamental subtlety: residuals are not errors. Errors are unobservable; residuals are what remains after projecting Y onto the space spanned by the regressors. Their variance-covariance matrix turns out to be the error variance multiplied by the identity minus the projection matrix (the so-called hat matrix). The consequence is that residuals can have different variances from one another even if the true errors all share the same variance: the variance of the fitted values is smaller near the mean and larger at the extremes.

Diagnostic plots. The most direct way to examine model adequacy is to visualize the residuals: against the fitted values, against the explanatory variable, and to look at their histogram, their box plot and their Q-Q plot. If there are no violations of UR.2 or UR.3, the plots should reveal no pattern; the histogram and the Q-Q plot should look approximately normal for UR.4 to hold. A Q-Q plot displays the quantiles of the data against those of a theoretical distribution: if the data come from that distribution, the points form approximately a straight line.

import statsmodels.api as sm
import scipy.stats as stats
from matplotlib import pyplot as plt

fig = sm.qqplot(model.resid.values, stats.t, distargs=(6,), fit=True, line='45')
plt.show()

Since there is no guarantee that we picked the right functional form, residual plots are a great clue about what we have missed. Transforming variables can ease the problem: if there are non-linearities, deal with them first and only then check for constant variance; transforming Y mainly attacks problems with the error term, while transforming X mainly attacks non-linearity. Sometimes transforming Y fixes the non-linearity but ruins normality, and sometimes transforming X achieves nothing at all.

Heteroskedasticity. If the error variance conditional on X is the same for every observation, the residuals are homoskedastic; if not, they are heteroskedastic. The consequences are:

• The OLS parameters remain unbiased.
• The estimates stop being efficient: OLS gives equal weight to every observation when in fact those with larger variance carry less information.
• The estimate of the residual variance is biased and therefore the standard errors are biased, which drags bias into the test statistics and intervals.
• Because of that bias we may fail to reject the null of a zero coefficient when it is actually false: a type II error.

Possible corrections: take logarithms, switch method (for instance weighted least squares, which assigns different weights to each observation), or use White's robust standard errors, which leave the coefficients untouched.

Tests. Goldfeld-Quandt splits the sample into two subsets ordered by a known explanatory variable and compares their residual variances with an F test; it requires the variance to be a monotone function of that variable and does not distinguish between heteroskedasticity and a specification problem. Breusch-Pagan regresses the squared residual on the explanatory variables and tests whether that coefficient is zero; the statistic is N times the R squared of the auxiliary regression, with a chi squared distribution. White is more general because it allows non-linear effects (squares and cross products), although it loses power with many regressors.

import statsmodels.stats.diagnostic as sm_diagnostic
print(sm_diagnostic.het_breuschpagan(model.resid, model.model.exog))
print(sm_diagnostic.het_white(model.resid, model.model.exog))