stats.py

#!/usr/bin/env python3
"""
ols_tstats.py

Compute OLS estimates, standard errors, t-statistics and p-values
for a linear model Y = H * Theta + noise.

Inputs:
    Y : numpy array (n x 1) observed values
    H : numpy array (n x k) design matrix
"""

import numpy as np
from scipy.stats import t as student_t

def ols_tstats(Y: np.ndarray, H: np.ndarray, theta_hat, XtX_inv):
    n, k = H.shape
    # Residuals and variance estimate
    residuals = Y - H @ theta_hat
    sigma2_hat = (residuals.T @ residuals) / (n - k)

    # Standard errors
    var_theta = sigma2_hat * XtX_inv
    se_theta = np.sqrt(np.diag(var_theta))

    # t-statistics
    t_stats = theta_hat.flatten() / se_theta

    # p-values (two-sided)
    p_values = 2 * (1 - student_t.cdf(np.abs(t_stats), df=n - k))

    return se_theta, t_stats, p_values

if __name__ == "__main__":
    # Example usage
    np.random.seed(0)
    n = 30
    # Example: intercept + 2 predictors
    X1 = np.random.rand(n)
    X2 = np.random.rand(n)
    H = np.column_stack([np.ones(n), X1, X2])
    Y = 1.0 + 2.0*X1 - 0.5*X2 + 0.1*np.random.randn(n)

    theta_hat, se_theta, t_stats, p_values = ols_tstats(Y, H)

    print("Coefficients:", theta_hat)
    print("Std Errors :", se_theta)
    print("t-stats    :", t_stats)
    print("p-values   :", p_values)