What is heteroskedasticity? How to spot it and what to do

A regression error is heteroskedastic when its variance differs across observations. The typical picture is a fan: the points hug the regression line at one end and spread out at the other. Heteroskedasticity does not bias ordinary least squares (OLS) coefficients, but it makes the usual standard errors, t statistics and confidence intervals unreliable.

This article shows what heteroskedasticity looks like, why it matters, how to detect it and what to do about it, using a simulation you can run with OpenEconometrics. The glossary entry gives the short definition.

What heteroskedasticity looks like

Think of household saving. A household earning 1,000 a month has little room to save or to dissave. A household earning 10,000 can save a large share of its income or spend all of it. Average saving may rise steadily with income, while the spread around that average grows too.

In a regression yᵢ = β₀ + β₁xᵢ + uᵢ, the errors are homoskedastic when their variance is the same for every observation and heteroskedastic when it is not:

homoskedastic:    Var(uᵢ | xᵢ) = σ²
heteroskedastic:  Var(uᵢ | xᵢ) = σᵢ²
in this article:  σᵢ = 0.1 · xᵢ

The simulated data in this article follow the last two lines. Saving rises by 0.15 for each extra unit of income, and the standard deviation of the error is one tenth of income.

Saving spreads out as income rises in 500 simulated households

A worked simulation

The script draws 500 synthetic households with right-skewed incomes (median about 3,000 a month), fits OLS and runs the diagnostics discussed below. Save it as heteroskedasticity.py, because the later scripts import its simulate() function.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=500, seed=2026):
    """Synthetic households whose saving varies more at higher incomes."""
    rng = np.random.default_rng(seed)
    income = np.exp(np.log(3) + 0.4 * rng.standard_normal(n))  # thousands a month
    error = 0.1 * income * rng.standard_normal(n)               # spread grows with income
    saving = -0.2 + 0.15 * income + error
    return pd.DataFrame({"income": income, "saving": saving})


if __name__ == "__main__":
    df = simulate()
    fit = oe.ols(data=df, y="saving", x=["income"])
    print(fit.summary())

    # Residual spread in the lowest and highest quarter of fitted values
    fitted = fit.predict(kind="xb")["xb"].to_numpy()
    resid = fit.predict(kind="residuals")["residuals"].to_numpy()
    quarter = pd.qcut(fitted, 4, labels=False)
    print("Residual SD by fitted-value quarter:",
          [round(float(resid[quarter == q].std()), 3) for q in range(4)])

    for name, test in [("Breusch-Pagan", fit.hettest()),
                       ("Koenker", fit.hettest(method="iid")),
                       ("White", fit.white_test())]:
        print(f"{name}: chi2({test['df']}) = {test['statistic']:.2f}, p = {test['p_value']:.3g}")

    robust = oe.ols(data=df, y="saving", x=["income"], covariance="HC3")
    for label, model in [("classical", fit), ("HC3", robust)]:
        b = model.coefficients[1]
        print(f"{label:>9}: slope {b.estimate:.4f}, se {b.std_error:.4f}, "
              f"95% CI [{b.ci_low:.4f}, {b.ci_high:.4f}]")

With openecon 0.3.20a1 from PyPI, the regression part of the output is:

OLS — saving
Observations: 500  |  Covariance: nonrobust  |  Confidence: 95%

Term        Estimate  Std. error         t     P>|stat|   CI lower    CI upper
---------  ---------  ----------  --------  -----------  ---------  ----------
Intercept  -0.138576   0.0420278  -3.29724   0.00104615  -0.221149  -0.0560021
income       0.13217    0.011586   11.4078  6.11283e-27   0.109407    0.154934

r_squared: 0.207179  |  adjusted_r_squared: 0.205587  |  rmse: 0.379918  |  df_model: 1
df_resid: 498  |  ss_resid: 71.8801  |  ss_model: 18.7837  |  ss_total: 90.6638
log_likelihood: -224.567  |  condition_number_scaled: 1  |  aic: 453.134  |  bic: 461.564
Model F test: F(1, 498) = 130.137, p = 6.113e-27

The estimated slope is 0.132, against a true value of 0.15. The classical standard error of 0.0116 makes the estimate look precise. The next section shows that it is less precise than that.

Why OLS stays unbiased but the standard errors go wrong

OLS is unbiased when the error has mean zero given the regressors. That assumption says nothing about the error's variance, so heteroskedasticity on its own does not bias the coefficients. What changes is the variance of the estimator. In general it is

Var(β̂ | X) = (X′X)⁻¹ X′ΩX (X′X)⁻¹
Ω = diag(σ₁², …, σₙ²)

Only when every σᵢ² equals a common σ² does this collapse to the classical formula σ²(X′X)⁻¹. The classical formula replaces each observation's own variance with an average. The slope is learned mostly from households far from the mean income, and here the high-income ones are also the noisiest. Averaging their large variances with the small ones understates the uncertainty.

Heteroskedasticity also means OLS is no longer the most efficient linear unbiased estimator. An estimator that gives less weight to noisy observations can do better, as the remedies below show.

A Monte Carlo experiment measures the size of the problem. The script keeps the 500 incomes fixed, draws new errors 1,000 times and refits the regression each time: OLS with classical standard errors, OLS with HC3 standard errors and weighted least squares (WLS), explained further down.

import numpy as np
import pandas as pd
import openecon as oe
from heteroskedasticity import simulate

if __name__ == "__main__":
    df = simulate()
    income = df["income"].to_numpy()
    df["w"] = 1 / income ** 2
    fits = {"OLS, classical": {}, "OLS, HC3": {"covariance": "HC3"},
            "WLS, classical": {"weights": "w", "weight_type": "aweight"}}
    rng = np.random.default_rng(7)
    rows = []
    for rep in range(1000):  # same incomes, new errors each time
        df["saving"] = -0.2 + 0.15 * income + 0.1 * income * rng.standard_normal(len(df))
        for label, options in fits.items():
            b = oe.ols(data=df, y="saving", x=["income"], **options).coefficients[1]
            rows.append({"fit": label, "slope": b.estimate, "se": b.std_error,
                         "covers": b.ci_low <= 0.15 <= b.ci_high})
    summary = pd.DataFrame(rows).groupby("fit", sort=False).agg(
        mean_slope=("slope", "mean"), sd_slope=("slope", "std"),
        mean_se=("se", "mean"), coverage=("covers", "mean"))
    print(summary.round(4).rename_axis(None))

                mean_slope  sd_slope  mean_se  coverage
OLS, classical      0.1506    0.0179   0.0110     0.794
OLS, HC3            0.1506    0.0179   0.0187     0.950
WLS, classical      0.1506    0.0107   0.0110     0.959

Classical standard errors understate how much the slope varies

Across the simulations the OLS slope averages 0.1506, close to the true 0.15, so there is no bias. But the estimates vary with a standard deviation of 0.0179, while the classical standard error averages 0.0110, about 38% too small. Nominal 95% confidence intervals therefore contain the true slope in only 79.4% of samples. The HC3 standard error averages 0.0187, and its intervals cover the true slope in 95.0% of samples.

How to detect heteroskedasticity

Plot the residuals

The first check is visual. Plot the residuals against the fitted values, and against each regressor you suspect. Under constant variance the cloud is a horizontal band of even width. A fan or funnel suggests the variance changes with the level of the outcome.

OLS residuals fan out as fitted saving rises

The first script puts a number on the funnel. The residual standard deviation rises from 0.197 in the quarter of households with the lowest fitted values to 0.579 in the quarter with the highest:

Residual SD by fitted-value quarter: [0.197, 0.258, 0.362, 0.579]

The Breusch–Pagan test and Koenker's version

The Breusch–Pagan test regresses the squared residuals, divided by their mean, on variables that might drive the variance and tests whether they explain any of it (Breusch and Pagan 1979; Cook and Weisberg 1983). By default fit.hettest() uses the fitted values as that variable. fit.hettest(rhs=True) uses all regressors instead, and fit.hettest(variables=[...]) picks specific terms.

The default statistic assumes normally distributed errors. Koenker (1981) proposed a studentized version, n·R² from the same auxiliary regression, which does not rely on normality. Select it with method="iid".

White's test

White (1980) proposed a test that does not require you to name the variables. It regresses the squared residuals on the regressors, their squares and their cross-products, and compares n·R² with a chi-squared distribution. With one regressor the auxiliary regression has two slope terms, income and income², so the statistic has 2 degrees of freedom.

The first script prints all three tests:

Breusch-Pagan: chi2(1) = 377.74, p = 3.86e-84
Koenker: chi2(1) = 132.58, p = 1.12e-30
White: chi2(2) = 188.12, p = 1.42e-41

All three reject constant variance decisively, as they should for data built with a strong fan. Real data are rarely this clear, and the tests have two blind spots. A rejection can come from a misspecified mean, such as a missing squared term, because the squared residuals then vary systematically with the regressors. A large p-value does not prove constant variance either, since a test that does not reject may simply lack power.

What to do about it

Use heteroskedasticity-robust standard errors

The simplest remedy keeps the OLS coefficients and replaces the standard errors with White's (1980) heteroskedasticity-consistent estimator or one of its refinements. The first script does this with covariance="HC3", and its last two lines compare the results:

classical: slope 0.1322, se 0.0116, 95% CI [0.1094, 0.1549]
      HC3: slope 0.1322, se 0.0245, 95% CI [0.0841, 0.1803]

The coefficient is identical; only the uncertainty changes. In this sample the HC3 standard error, 0.0245, is larger than its Monte Carlo average of 0.0187. A robust standard error is itself an estimate, and it can be noisy when a few observations have high leverage.

Long and Ervin (2000) recommend correcting for heteroskedasticity whenever there is reason to suspect it, without letting a screening test decide, and using HC3 when n ≤ 250. The guide to robust and clustered standard errors explains HC0 to HC3 and clustering, and the OLS method page lists every covariance option.

Weighted least squares when you know the variance

If you know how the variance changes, you can do better than OLS. Weighted least squares gives each observation a weight proportional to 1/σᵢ², so noisy observations count less. Here σᵢ is proportional to income, so the weight is 1/income². OpenEconometrics takes such weights as analytic weights, weight_type="aweight".

import numpy as np
import pandas as pd
import openecon as oe
from heteroskedasticity import simulate

if __name__ == "__main__":
    df = simulate()
    ols = oe.ols(data=df, y="saving", x=["income"])

    # Known variance form: sd(u) is proportional to income, so weight by 1 / income²
    df["w"] = 1 / df["income"] ** 2
    wls = oe.ols(data=df, y="saving", x=["income"], weights="w", weight_type="aweight")
    wls_hc3 = oe.ols(data=df, y="saving", x=["income"], weights="w", weight_type="aweight",
                     covariance="HC3")

    # Unknown form: model log(û²) as a linear function of log(income), then reweight
    u = ols.predict(kind="residuals")["residuals"].to_numpy()
    aux = pd.DataFrame({"log_u2": np.log(u ** 2), "log_income": np.log(df["income"])})
    variance = oe.ols(data=aux, y="log_u2", x=["log_income"])
    print(f"Variance model: log(û²) slope on log(income) = {variance.coefficients[1].estimate:.3f}")
    df["w_hat"] = 1 / np.exp(variance.predict(kind="xb")["xb"].to_numpy())
    fgls = oe.ols(data=df, y="saving", x=["income"], weights="w_hat", weight_type="aweight",
                  covariance="HC3")

    for label, model in [("WLS, 1/income²", wls), ("WLS, 1/income², HC3", wls_hc3),
                         ("FGLS, HC3", fgls)]:
        b = model.coefficients[1]
        print(f"{label:>20}: slope {b.estimate:.4f}, se {b.std_error:.4f}, "
              f"95% CI [{b.ci_low:.4f}, {b.ci_high:.4f}]")

Variance model: log(û²) slope on log(income) = 2.106
      WLS, 1/income²: slope 0.1491, se 0.0112, 95% CI [0.1271, 0.1712]
 WLS, 1/income², HC3: slope 0.1491, se 0.0122, 95% CI [0.1252, 0.1731]
           FGLS, HC3: slope 0.1483, se 0.0124, 95% CI [0.1239, 0.1727]

The WLS slope is 0.1491. Its classical standard error, 0.0112, is about as small as the OLS classical one, but this time it is honest. In the Monte Carlo the WLS slope varied with a standard deviation of 0.0107, against 0.0179 for OLS, and its 95% intervals covered the true slope in 95.9% of samples. The right weights cut the sampling standard deviation by about 40%.

Model the variance when you do not know it

Usually the form of the variance is unknown. Feasible generalized least squares (FGLS) estimates it first. If σᵢ is proportional to incomeᵢ raised to some power γ, then log σᵢ² is linear in log(incomeᵢ) with slope 2γ. The script regresses log(û²) on log(income) and gets a slope of 2.106, close to the true value of 2. It then weights each household by the inverse of its fitted variance.

A variance model can be wrong, so pair the weighted fit with robust standard errors. Romano and Wolf (2017) show that WLS combined with heteroskedasticity-consistent standard errors gives asymptotically valid inference even when the weighting model is misspecified. Here the FGLS slope is 0.1483 with an HC3 standard error of 0.0124.

Transform the outcome

Sometimes the error acts multiplicatively: a positive outcome equals a systematic part times a random factor that does not depend on the regressors. The spread of y then grows with its level, and modelling log(y) instead of y turns the random factor into an additive error with constant variance. The transformation changes the question: the coefficients then describe approximate proportional changes, and predictions of y need care when converted back from logs. Logs also require a positive outcome. They are not an option for the saving data here, where many households dissave, as the fan chart shows.

What can go wrong

• Treating robust standard errors as a fix for bias. They repair the uncertainty, not the coefficients. If the mean is misspecified or a regressor is endogenous, the estimate is still wrong; see correlation versus causation.
• Letting a pretest choose the standard errors. A non-significant Breusch–Pagan test is weak evidence of constant variance, and Long and Ervin (2000) advise against making the decision this way.
• Reading every rejection as heteroskedasticity. An omitted nonlinearity can produce the same pattern in the squared residuals. Plot the residuals against each regressor and try a more flexible mean before concluding.
• Trusting WLS weights blindly. Weights that depend only on the regressors do not bias WLS under the usual mean assumption, but wrong weights make its classical standard errors wrong. Report robust standard errors unless the weights are known.
• Ignoring dependence. Heteroskedasticity-robust standard errors assume independent observations. Grouped data need clustered standard errors and time series need HAC standard errors.

A short checklist

1. Plot the residuals against the fitted values and against key regressors.
2. Run fit.hettest() (or Koenker's method="iid" when errors may not be normal) and fit.white_test(), and read any rejection alongside the plots.
3. Report heteroskedasticity-robust standard errors by default, and use HC3 when n ≤ 250.
4. If you know or can model the variance, consider WLS or FGLS for precision, and keep robust standard errors.
5. Consider logs for a positive outcome whose spread grows with its level, and interpret the coefficients accordingly.
6. State the covariance estimator in every table you report.

References

• Breusch, T. S. and Pagan, A. R. (1979). A simple test for heteroscedasticity and random coefficient variation. Econometrica, 47(5), 1287–1294. doi:10.2307/1911963
• Cook, R. D. and Weisberg, S. (1983). Diagnostics for heteroscedasticity in regression. Biometrika, 70(1), 1–10. doi:10.1093/biomet/70.1.1
• Koenker, R. (1981). A note on studentizing a test for heteroscedasticity. Journal of Econometrics, 17(1), 107–112. doi:10.1016/0304-4076(81)90062-2
• Long, J. S. and Ervin, L. H. (2000). Using heteroscedasticity consistent standard errors in the linear regression model. The American Statistician, 54(3), 217–224. doi:10.1080/00031305.2000.10474549
• Romano, J. P. and Wolf, M. (2017). Resurrecting weighted least squares. Journal of Econometrics, 197(1), 1–19. doi:10.1016/j.jeconom.2016.10.003
• White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica, 48(4), 817–838. doi:10.2307/1912934