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"Error: {name} must be a 2D list."
        if not data:
            return f"Error: {name} cannot be empty."
        for i, row in enumerate(data):
            if not isinstance(row, list):
                return f"Error: {name} row {i} must be a list."
            if not row:
                return f"Error: {name} row {i} cannot be empty."
            for j, val in enumerate(row):
                if not isinstance(val, (int, float)):
                    return f"Error: {name}[{i}][{j}] must be numeric."
                if math.isnan(val) or math.isinf(val):
                    return f"Error: {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'

    try:
        # 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"Error: diag_test must be one of {valid_tests}."

        # Convert to numpy arrays
        resid_array = np.array(residuals_2d, dtype=float)
        exog_array = np.array(exog_2d, dtype=float)

        # Flatten residuals to 1D array
        if resid_array.shape[1] != 1:
            return "Error: 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"Error: number of residuals ({n_obs}) must match number of observations in exog ({exog_array.shape[0]})."

        if n_obs < 3:
            return "Error: 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
        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])

        # 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
    except Exception as exc:
        return f"Error: {exc}"

Overview

Excel Usage

Examples

Python Code

Online Calculator