Run OLS regressions with full diagnostics: heteroscedasticity tests, robust/clustered SEs, VIF, structural breaks, and publication-ready tables via statsmodels.
Run OLS regressions with full diagnostics: heteroscedasticity tests, robust/clustered SEs, VIF, structural breaks, and publication-ready tables via statsmodels.
Ordinary Least Squares is the workhorse of empirical economics. This skill covers the full
pipeline: model specification, assumption testing, robust inference, and presentation. It
follows best practices from Angrist & Pischke (2009) and Greene (2018).
Core Concepts
Why OLS?
Under the Gauss-Markov assumptions (linearity, random sampling, no perfect multicollinearity,
zero conditional mean, homoscedasticity), OLS is BLUE — Best Linear Unbiased Estimator. In
practice, homoscedasticity almost never holds for economic cross-sectional data, so robust
standard errors are the default. Consistency requires only that E[u|X] = 0.
"""
Run Breusch-Pagan and White's tests for heteroscedasticity.
Breusch-Pagan: tests whether residual variance is a linear function of regressors.
White's test: tests against general heteroscedasticity (includes cross-products).
Returns a dict with test statistics, p-values, and plain-English interpretations.
"""
"""
Re-estimate with the requested covariance sandwich estimator.
Parameters
----------
result : OLS result from statsmodels (plain, non-robust)
cov_type : 'HC0', 'HC1', 'HC2', 'HC3' (White), or 'cluster'
cluster_var: column name for clustering (only when cov_type='cluster')
df : original DataFrame (needed for cluster variable)
Returns
-------
DataFrame with columns: coef, se_ols, se_robust, t_robust, p_robust,
ci_lower, ci_upper, stars
"""
"""
Compute Variance Inflation Factors for all continuous regressors.
VIF = 1 / (1 - R²_j), where R²_j is from regressing X_j on all other regressors.
Rule of thumb: VIF > 10 indicates serious multicollinearity.
"""
"""
Chow test for parameter stability across two sub-samples.
H0: coefficients are the same in both groups.
Statistic: F = [(RSS_pool - RSS_1 - RSS_2) / k] / [(RSS_1 + RSS_2) / (n - 2k)]
Parameters
----------
formula : Patsy formula
df : full DataFrame
break_var : column name to split on
break_value : threshold — group1 is df[break_var] < break_value
"""
len
# number of parameters
len
2
1
2
return
"F_statistic"
"p_value"
"df_numerator"
"df_denominator"
2
"n_group1"
len
"n_group2"
len
"reject_H0"
0.05
"interpretation"
f"Structural break detected at {break_var}={break_value} (p={p:.4f})"
if
0.05
else
f"No structural break at {break_var}={break_value} (p={p:.4f})"
"""
Build a text regression table comparable to R's stargazer.
Parameters
----------
results_list : list of statsmodels result objects (use .fit(cov_type='HC3'))
model_names : list of column headers, e.g. ['(1)', '(2)', '(3)']
title : table title
dep_var_label: row label for dependent variable
Returns
-------
Formatted string suitable for printing or writing to a .txt file.
"""
if
is
None
f"({i+1})"
for
in
range
len
# collect all variable names across models
for
in
for
in
if
not
in
14
"="
20
len
f"{dep_var_label:<20}"
""
f"{m:>{col_w}}"
for
in
"-"
20
len
for
in
f"{var:<20}"
f"{'':20}"
for
in
if
in
f"{float_fmt.format(coef)+star:>{col_w}}"
f"{'('+float_fmt.format(se)+')':>{col_w}}"
else
f"{'':>{col_w}}"
f"{'':>{col_w}}"
"-"
20
len
# footer stats
f"{'Observations':<20}"
""
f"{int(res.nobs):>{col_w}}"
for
in
f"{'R²':<20}"
""
f"{float_fmt.format(res.rsquared):>{col_w}}"
for
in
f"{'Adj. R²':<20}"
""
f"{float_fmt.format(res.rsquared_adj):>{col_w}}"
for
in
"Note: * p<0.1 ** p<0.05 *** p<0.01 Robust (HC3) standard errors in parentheses"
# example_a_wage_regression.py"""
Mincer wage regression with heteroscedasticity diagnostics.
Dataset: CPS-style synthetic wage data.
"""import numpy as np
import pandas as pd
from ols_regression import run_ols_full, make_regression_table, plot_diagnostics
rng = np.random.default_rng(42)
n = 2000# simulate data
educ = rng.integers(8, 21, n).astype(float)
exper = np.clip(rng.normal(20, 10, n), 0, 45)
female = rng.binomial(1, 0.48, n).astype(float)
union = rng.binomial(1, 0.15, n).astype(float)
# log-wage DGP: heteroscedastic — variance rises with education
sigma = 0.2 + 0.03 * educ
u = rng.normal(0, sigma, n)
lnwage = (1.2 + 0.10 * educ + 0.04 * exper
- 0.0006 * exper**2 - 0.22 * female + 0.12 * union + u)
df = pd.DataFrame({
"lnwage": lnwage, "educ": educ, "exper": exper,
"exper2": exper**2, "female": female, "union": union,
})
# ── Model 1: parsimonious ──
m1 = run_ols_full("lnwage ~ educ + exper + exper2 + female", df)
# ── Model 2: add union ──
m2 = run_ols_full("lnwage ~ educ + exper + exper2 + female + union", df)
# print diagnostic summaryprint(m1["summary"])
# regression table
table = make_regression_table(
[m1["result"], m2["result"]],
model_names=["(1) Base", "(2) + Union"],
dep_var_label="ln(wage)",
)
print(table)
# Marginal effect of education (partial derivative, evaluated at mean exper):# ∂ ln(wage) / ∂ educ ≈ β_educ → wage rises by ~10% per extra year of schooling# returns to education (percentage)
beta_educ = m2["result"].params["educ"]
print(f"\nReturn to education: {beta_educ*100:.1f}% per year of schooling")
# peak experience: β_exper / (2 * |β_exper2|)
b1 = m2["result"].params["exper"]
b2 = m2["result"].params["exper2"]
peak = -b1 / (2 * b2)
print(f"Peak experience (years): {peak:.1f}")
# save diagnostic plots
fig = plot_diagnostics(m2["result_ols"], title="Wage Regression Diagnostics")
fig.savefig("wage_diagnostics.png", dpi=150, bbox_inches="tight")
print("Saved wage_diagnostics.png")
Example B — OLS vs Robust SE with Heteroscedastic Data
Demonstrates how plain OLS standard errors understate uncertainty when
variance depends on a regressor.
# example_b_ols_vs_robust.py"""
Compare OLS, HC1, HC3, and clustered SEs on heteroscedastic data.
Shows that t-statistics can be inflated up to 2× with plain OLS.
"""import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from ols_regression import (
run_ols_full,
get_robust_se,
check_heteroscedasticity,
)
rng = np.random.default_rng(0)
n = 500# DGP: Var(u|x) = (1 + 2x)² — strong heteroscedasticity
x = rng.uniform(0, 5, n)
u = rng.normal(0, 1 + 2 * x, n) # heteroscedastic errors
y = 2 + 1.5 * x + u
state = rng.integers(0, 20, n) # 20 state clusters
df = pd.DataFrame({"y": y, "x": x, "state": state})
# plain OLS
res_plain = run_ols_full("y ~ x", df, cov_type="nonrobust")
het = check_heteroscedasticity(res_plain["result_ols"])
print("Breusch-Pagan:", het["breusch_pagan"]["interpretation"])
print("White's test :", het["white"]["interpretation"])
print()
# compare SE across estimators
se_ols = get_robust_se(res_plain["result_ols"], cov_type="nonrobust")
se_hc1 = get_robust_se(res_plain["result_ols"], cov_type="HC1")
se_hc3 = get_robust_se(res_plain["result_ols"], cov_type="HC3")
se_cluster = get_robust_se(
res_plain["result_ols"], cov_type="cluster",
cluster_var="state", df=df
)
comparison = pd.DataFrame({
"OLS SE": se_ols["se_ols"],
"HC1 SE": se_hc1["se_robust"],
"HC3 SE": se_hc3["se_robust"],
"Clustered SE": se_cluster["se_robust"],
}).loc[["x"]]
print("Standard error comparison for coefficient on x:")
print(comparison.to_string())
# visualise
fig, axes = plt.subplots(1, 2, figsize=(12, 5))
axes[0].scatter(x, u, alpha=0.3, s=15, color="#2c7bb6")
axes[0].axhline(0, color="red", linestyle="--")
axes[0].set_xlabel("x"); axes[0].set_ylabel("Residual u")
axes[0].set_title("Heteroscedastic residuals")
labels = ["OLS", "HC1", "HC3", "Clustered"]
values = [
float(se_ols.loc["x", "se_ols"]),
float(se_hc1.loc["x", "se_robust"]),
float(se_hc3.loc["x", "se_robust"]),
float(se_cluster.loc["x", "se_robust"]),
]
colors = ["#d7191c", "#2c7bb6", "#1a9641", "#ff7f00"]
bars = axes[1].bar(labels, values, color=colors, width=0.5)
axes[1].set_ylabel("Standard error of β̂_x")
axes[1].set_title("SE comparison: OLS vs robust estimators")
for bar, val inzip(bars, values):
axes[1].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.002,
f"{val:.4f}", ha="center", va="bottom", fontsize=9)
plt.tight_layout()
fig.savefig("se_comparison.png", dpi=150, bbox_inches="tight")
print("\nSaved se_comparison.png")
Coefficient Interpretation Guide
Specification
Interpretation of β
y = β₀ + βX + u
One unit ↑ X → β units ↑ y
ln(y) = β₀ + βX + u
One unit ↑ X → 100β% ↑ y
y = β₀ + β ln(X) + u
1% ↑ X → β/100 units ↑ y
ln(y) = β₀ + β ln(X) + u
1% ↑ X → β% ↑ y (elasticity)
y = β₀ + βX + γX² + u
∂y/∂X = β + 2γX (non-linear)
y = β₀ + β D + u (D binary)
Being in group D → β units ↑ y
Omitted Variable Bias Formula
If the true model is y = β₀ + β₁X₁ + β₂X₂ + u but you omit X₂:
plim(β̂₁^short) = β₁ + β₂ · δ₁₂
where δ₁₂ = Cov(X₂, X₁) / Var(X₁) is the regression coefficient of X₂ on X₁.
Direction of bias: positive if β₂ and Corr(X₁,X₂) have the same sign.
Checklist Before Reporting Results
Report N, R², adjusted R², F-statistic
Use HC3 or clustered SE as default (not plain OLS SE)
Check VIF — flag any variable > 10
Run RESET test for functional form misspecification
Acknowledge endogeneity if present (consider IV)
Show regression table with at least one robustness column
Inspect residual plots for patterns
Note economically meaningful effect sizes, not just p-values
References
Angrist, J. D., & Pischke, J.-S. (2009). Mostly Harmless Econometrics. Princeton UP.
Greene, W. H. (2018). Econometric Analysis (8th ed.). Pearson.
White, H. (1980). A heteroskedasticity-consistent covariance matrix estimator. Econometrica, 48(4), 817–838.
MacKinnon, J. G., & White, H. (1985). Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties. Journal of Econometrics, 29(3), 305–325.
Stock, J. H., & Watson, M. W. (2020). Introduction to Econometrics (4th ed.). Pearson.
