Statsmodels OLS Regression and Tests in Python

The Python statsmodels library provides a comprehensive framework for estimating Ordinary Least Squares (OLS) regression models and conducting rigorous statistical tests. By coupling linear algebra solvers with classical inferential statistics, statsmodels allows practitioners to fit linear models, estimate parameters, and automatically generate detailed summary outputs containing p-values, standard errors, and confidence intervals. Beyond basic model fitting, the library includes a suite of diagnostic tests to evaluate core regression assumptions, such as residual normality, homoscedasticity, and independence.

Setting Up OLS Regression in Statsmodels

To perform OLS regression, statsmodels requires the dependent variable (\(y\)) and independent variables (\(X\)) formatted as NumPy arrays or pandas DataFrames. Unlike libraries such as scikit-learn, statsmodels does not include an intercept term by default. You must explicitly add a constant column using sm.add_constant().

import numpy as np
import statsmodels.api as sm

# Sample data
np.random.seed(42)
X = np.random.rand(100, 2)
y = 2 + 3 * X[:, 0] + 1.5 * X[:, 1] + np.random.normal(0, 0.5, 100)

# Add an intercept (constant term)
X_with_const = sm.add_constant(X)

# Define and fit the OLS model
model = sm.OLS(y, X_with_const)
results = model.fit()

The fit() method uses standard matrix operations—specifically solving the normal equations via QR decomposition or singular value decomposition (SVD)—to calculate the coefficient estimates that minimize the sum of squared residuals.

Interpreting the Statistical Summary

Calling results.summary() produces a structured report divided into three main components:

Model Fit Metrics: Displays the coefficient of determination (\(R^2\)), Adjusted \(R^2\), the overall \(F\)-statistic, and its corresponding p-value (Prob (F-statistic)), which evaluates whether any independent variable significantly predicts the outcome.
Parameter Estimates: Details the estimated coefficients (coef), standard errors (std err), \(t\)-statistics (t), two-tailed p-values (P>|t|), and 95% confidence intervals for each predictor.
Residual Diagnostics: Provides baseline metrics on error behavior, including the Omnibus test, the Jarque-Bera test, and the Durbin-Watson statistic.

Hypothesis Testing on Coefficients

To test linear hypotheses about parameters, statsmodels provides dedicated methods on the fitted results object:

t-tests (t_test): Tests individual linear combinations of coefficients (e.g., whether a coefficient equals a specific non-zero value or whether two coefficients are equal).
F-tests (f_test): Tests multiple linear restrictions simultaneously, such as testing whether a subset of coefficients jointly equals zero.

# Test if the first predictor equals 3
t_test_res = results.t_test("x1 = 3")

# Joint hypothesis test (F-test)
f_test_res = results.f_test("x1 = 3, x2 = 1.5")

Residual Diagnostic Tests

Validating the Gauss-Markov assumptions is essential for ensuring that OLS estimators remain the Best Linear Unbiased Estimators (BLUE). The statsmodels.stats module includes tests for these assumptions:

1. Autocorrelation

The Durbin-Watson statistic tests for first-order autocorrelation in residuals and is printed in the standard summary. A value near 2 indicates no autocorrelation, while values approaching 0 or 4 indicate positive or negative autocorrelation, respectively. For higher-order autocorrelation, the Breusch-Godfrey test is available via sm.stats.acorr_breusch_godfrey(results).

2. Heteroscedasticity

When the variance of residuals is non-constant, standard errors are biased. statsmodels provides tests to detect heteroscedasticity:

Breusch-Pagan Test: sm.stats.het_breuschpagan(results.resid, results.model.exog) tests whether residual variance depends on the independent variables.
White Test: sm.stats.het_white(results.resid, results.model.exog) includes cross-products and squared terms for a more general check.

If heteroscedasticity is detected, robust standard errors (such as HC1, HC2, or HC3) can be calculated directly by refitting with model.fit(cov_type='HC3').

3. Normality of Residuals

Normality is required for valid hypothesis testing in small samples. The summary output includes:

Omnibus Test: Measures skewness and kurtosis jointly.
Jarque-Bera Test: Evaluates whether sample skewness and kurtosis match a normal distribution.

Both tests can be run explicitly using sm.stats.jarque_bera(results.resid).

4. Multicollinearity

High correlation between predictors inflates standard errors. The Variance Inflation Factor (VIF) measures this inflation:

from statsmodels.stats.outliers_influence import variance_inflation_factor

vifs = [
    variance_inflation_factor(X_with_const, i)
    for i in range(X_with_const.shape[1])
]

A VIF greater than 5 or 10 typically suggests problematic multicollinearity.