OLS_DIAGNOSTICS

Overview

The OLS_DIAGNOSTICS function performs a suite of diagnostic tests on Ordinary Least Squares (OLS) regression residuals to assess whether the standard regression assumptions are satisfied. Validating these assumptions is critical because the reliability of OLS estimates depends on conditions including homoskedasticity, absence of autocorrelation, and normally distributed errors.

This implementation uses the statsmodels library, which provides a comprehensive set of regression diagnostic tests. The function performs three key diagnostic tests:

Heteroskedasticity: Breusch-Pagan Test

The Breusch-Pagan test, developed by Trevor Breusch and Adrian Pagan (1979), tests whether the variance of regression errors is constant across observations. The test regresses squared residuals on the original predictors and uses a Lagrange Multiplier (LM) statistic that follows a chi-squared distribution under the null hypothesis of homoskedasticity. A significant p-value (< 0.05) indicates heteroskedasticity is present, which may require weighted least squares or heteroskedasticity-consistent standard errors.

Autocorrelation: Ljung-Box Test

The Ljung-Box test, developed by Greta M. Ljung and George E. P. Box (1978), tests whether residuals exhibit serial correlation. The test statistic is:

Q = n(n+2) \sum_{k=1}^{h} \frac{\hat{\rho}_k^2}{n-k}

where n is the sample size, \hat{\rho}_k is the sample autocorrelation at lag k, and h is the number of lags tested. Under the null hypothesis of no autocorrelation, Q follows a chi-squared distribution with h degrees of freedom.

Normality: Jarque-Bera Test

The Jarque-Bera test assesses whether residuals follow a normal distribution by examining skewness (S) and kurtosis (K):

JB = \frac{n}{6}\left(S^2 + \frac{(K-3)^2}{4}\right)

The test statistic asymptotically follows a chi-squared distribution with two degrees of freedom. A significant result indicates non-normality, which may affect the validity of hypothesis tests and confidence intervals derived from the regression.

The function returns a table containing each test name, test statistic, p-value, and an interpretation based on a 0.05 significance threshold.

This example function is provided as-is without any representation of accuracy.

Excel Usage

=OLS_DIAGNOSTICS(residuals, exog, diag_test)

residuals (list[list], required): Column vector of regression residuals from an OLS model.
exog (list[list], required): Matrix of independent variables (predictors) from the OLS model.
diag_test (str, optional, default: “all”): Which diagnostic test(s) to perform.

Returns (list[list]): 2D list with diagnostic tests, or error message string.

Examples

Example 1: Demo case 1

Inputs:

residuals	exog		diag_test
0.1	1	1	all
0.2	1	2
-0.1	1	3
0.05	1	4
-0.15	1	5
0.08	1	6
-0.05	1	7
0.12	1	8
-0.08	1	9
0.03	1	10

Excel formula:

=OLS_DIAGNOSTICS({0.1;0.2;-0.1;0.05;-0.15;0.08;-0.05;0.12;-0.08;0.03}, {1,1;1,2;1,3;1,4;1,5;1,6;1,7;1,8;1,9;1,10}, "all")

Expected output:

test_name	statistic	p_value	interpretation
Breusch-Pagan	2.306	0.1289	No problem detected
Ljung-Box	4.819	0.08988	No problem detected
Jarque-Bera	0.5026	0.7778	Normal

Example 2: Demo case 2

Inputs:

residuals	exog		diag_test
0.5	1	1	heteroskedasticity
-0.3	1	2
0.2	1	3
-0.1	1	4
0.4	1	5
-0.2	1	6

Excel formula:

=OLS_DIAGNOSTICS({0.5;-0.3;0.2;-0.1;0.4;-0.2}, {1,1;1,2;1,3;1,4;1,5;1,6}, "heteroskedasticity")

Expected output:

test_name	statistic	p_value	interpretation
Breusch-Pagan	1.564	0.2111	No problem detected

Example 3: Demo case 3

Inputs:

residuals	exog		diag_test
0.01	1	2	normality
0.02	1	3
-0.01	1	4
0.03	1	5
-0.02	1	6
0.015	1	7
-0.015	1	8
0.025	1	9

Excel formula:

=OLS_DIAGNOSTICS({0.01;0.02;-0.01;0.03;-0.02;0.015;-0.015;0.025}, {1,2;1,3;1,4;1,5;1,6;1,7;1,8;1,9}, "normality")

Expected output:

test_name	statistic	p_value	interpretation
Jarque-Bera	0.8672	0.6482	Normal

Example 4: Demo case 4

Inputs:

residuals	exog		diag_test
0.1	1	1.5	autocorrelation
0.15	1	2.5
0.12	1	3.5
0.18	1	4.5
0.14	1	5.5
0.16	1	6.5
0.13	1	7.5
0.17	1	8.5

Excel formula:

=OLS_DIAGNOSTICS({0.1;0.15;0.12;0.18;0.14;0.16;0.13;0.17}, {1,1.5;1,2.5;1,3.5;1,4.5;1,5.5;1,6.5;1,7.5;1,8.5}, "autocorrelation")

Expected output:

test_name	statistic	p_value	interpretation
Ljung-Box	1.957	0.1618	No problem detected

Python Code

import math
import numpy as np
from statsmodels.stats.diagnostic import het_breuschpagan, acorr_ljungbox
from statsmodels.stats.stattools import jarque_bera

def ols_diagnostics(residuals, exog, diag_test='all'):
    """
    Performs diagnostic tests on OLS regression residuals.

    See: https://www.statsmodels.org/stable/diagnostic.html

    This example function is provided as-is without any representation of accuracy.

    Args:
        residuals (list[list]): Column vector of regression residuals from an OLS model.
        exog (list[list]): Matrix of independent variables (predictors) from the OLS model.
        diag_test (str, optional): Which diagnostic test(s) to perform. Valid options: All Tests, Heteroskedasticity, Autocorrelation, Normality. Default is 'all'.

    Returns:
        list[list]: 2D list with diagnostic tests, or error message string.
    """
    def to2d(x):
        return [[x]] if not isinstance(x, list) else x

    def validate_2d_numeric(data, name):
        # Validate 2D list structure
        if not isinstance(data, list):
            return f"Invalid input: {name} must be a 2D list."
        if not data:
            return f"Invalid input: {name} cannot be empty."
        for i, row in enumerate(data):
            if not isinstance(row, list):
                return f"Invalid input: {name} row {i} must be a list."
            if not row:
                return f"Invalid input: {name} row {i} cannot be empty."
            for j, val in enumerate(row):
                if not isinstance(val, (int, float)):
                    return f"Invalid input: {name}[{i}][{j}] must be numeric."
                if math.isnan(val) or math.isinf(val):
                    return f"Invalid input: {name}[{i}][{j}] must be finite."
        return None

    def interpret_pvalue(p_value, test_name):
        # Interpretation based on p-value threshold of 0.05
        if test_name in ['heteroskedasticity', 'autocorrelation']:
            if p_value < 0.05:
                return 'Problem detected'
            else:
                return 'No problem detected'
        elif test_name == 'normality':
            if p_value < 0.05:
                return 'Non-normal'
            else:
                return 'Normal'
        return 'N/A'

    # Normalize inputs
    residuals_2d = to2d(residuals)
    exog_2d = to2d(exog)

    # Validate inputs
    error = validate_2d_numeric(residuals_2d, 'residuals')
    if error:
        return error

    error = validate_2d_numeric(exog_2d, 'exog')
    if error:
        return error

    # Validate test_type parameter
    valid_tests = ['all', 'heteroskedasticity', 'autocorrelation', 'normality']
    if diag_test not in valid_tests:
        return f"Invalid input: diag_test must be one of {valid_tests}."

    # Convert to numpy arrays
    try:
        resid_array = np.array(residuals_2d, dtype=float)
        exog_array = np.array(exog_2d, dtype=float)
    except Exception as exc:
        return f"Invalid input: unable to convert inputs to arrays: {exc}"

    # Flatten residuals to 1D array
    if resid_array.shape[1] != 1:
        return "Invalid input: residuals must be a column vector (2D list with 1 column)."
    resid_1d = resid_array.flatten()

    # Validate dimensions
    n_obs = len(resid_1d)
    if exog_array.shape[0] != n_obs:
        return f"Invalid input: number of residuals ({n_obs}) must match number of observations in exog ({exog_array.shape[0]})."

    if n_obs < 3:
        return "Invalid input: insufficient data for diagnostic tests (need at least 3 observations)."

    # Initialize results
    results = [['test_name', 'statistic', 'p_value', 'interpretation']]

    # Determine which tests to run
    tests_to_run = []
    if diag_test == 'all':
        tests_to_run = ['heteroskedasticity', 'autocorrelation', 'normality']
    else:
        tests_to_run = [diag_test]

    # Run requested tests
    try:
        for test in tests_to_run:
            if test == 'heteroskedasticity':
                # Breusch-Pagan test
                lm_stat, lm_pvalue, f_stat, f_pvalue = het_breuschpagan(resid_1d, exog_array)
                interp = interpret_pvalue(lm_pvalue, 'heteroskedasticity')
                results.append(['Breusch-Pagan', float(lm_stat), float(lm_pvalue), interp])

            elif test == 'autocorrelation':
                # Ljung-Box test (using lag=min(10, n_obs//5))
                max_lag = min(10, max(1, n_obs // 5))
                lb_result = acorr_ljungbox(resid_1d, lags=max_lag, return_df=True)
                # Use the last lag result
                last_row = lb_result.iloc[-1]
                lb_stat = float(last_row['lb_stat'])
                lb_pvalue = float(last_row['lb_pvalue'])
                interp = interpret_pvalue(lb_pvalue, 'autocorrelation')
                results.append(['Ljung-Box', lb_stat, lb_pvalue, interp])

            elif test == 'normality':
                # Jarque-Bera test
                jb_stat, jb_pvalue, skew, kurtosis = jarque_bera(resid_1d)
                interp = interpret_pvalue(jb_pvalue, 'normality')
                results.append(['Jarque-Bera', float(jb_stat), float(jb_pvalue), interp])

    except Exception as exc:
        return f"Error running diagnostic tests: {exc}"

    # Validate results
    for row in results[1:]:
        if not isinstance(row[1], (int, float)) or math.isnan(row[1]) or math.isinf(row[1]):
            return "Error: test produced invalid statistic."
        if not isinstance(row[2], (int, float)) or math.isnan(row[2]) or math.isinf(row[2]):
            return "Error: test produced invalid p-value."

    return results

Overview

Excel Usage

Examples

Python Code

Online Calculator