Regression Analysis for Financial Relationships › A Complete Regression Workflow

from scipy import stats


def complete_regression_analysis(y, X, variable_names=None):
    """
    Perform comprehensive regression analysis with diagnostics.

    Parameters:
    -----------
    y : array-like
        Dependent variable
    X : array-like
        Independent variables (without constant)
    variable_names : list, optional
        Names for the X variables

    Returns:
    --------
    dict : Dictionary containing model results and diagnostics
    """
    if variable_names is None:
        variable_names = [
            f"X{i}" for i in range(X.shape[1] if len(X.shape) > 1 else 1)
        ]

    # Convert to DataFrame if needed
    if isinstance(X, np.ndarray):
        X = pd.DataFrame(X, columns=variable_names)

    # Add constant
    X_const = sm.add_constant(X)

    # Fit models with different standard errors
    model_ols = sm.OLS(y, X_const).fit()
    model_hac = sm.OLS(y, X_const).fit(cov_type="HAC", cov_kwds={"maxlags": 4})

    # Diagnostics
    residuals = model_ols.resid

    # Durbin-Watson
    dw = durbin_watson(residuals)

    # Breusch-Pagan
    bp_stat, bp_pval, _, _ = het_breuschpagan(residuals, X_const)

    # Jarque-Bera (normality)
    jb_stat, jb_pval = stats.jarque_bera(residuals)

    results = {
        "model_ols": model_ols,
        "model_hac": model_hac,
        "durbin_watson": dw,
        "breusch_pagan": (bp_stat, bp_pval),
        "jarque_bera": (jb_stat, jb_pval),
        "residuals": residuals,
    }

    return results

# Apply to our three-factor data
analysis = complete_regression_analysis(
    ff_df["excess_return"], ff_df[["MKT", "SMB", "HML"]], ["MKT", "SMB", "HML"]
)