OLS from scratch

ols.py

import numpy as np

# Core functions
def ols(X, y):
    ret = np.dot(np.dot(np.linalg.inv(np.dot(X.transpose(), X)), X.transpose()), y)
    return ret.tolist()

def se(arr):
    return np.std(arr) / arr.size

def boot(X, y):
    n = X.shape[0]
    idx = np.arange(0, n-1)
    idx = np.random.choice(idx, size = n, replace = True)
    return X[idx,:],  y[idx]

def boot_summary(boots):
    ret = {}
    ret["coef"] = np.apply_along_axis(np.mean, 0, boots)
    ret["se"] = np.apply_along_axis(se, 0, boots)
    return ret
    
def ols_boot(X, y, n_boots=1000):
    boots = np.zeros(shape=(n_boots, X.shape[1]))
    for idx_boot in range(0, n_boots):
        x_boot, y_boot = boot(X, y)
        boots[idx_boot,:]  = ols(x_boot, y_boot)
    return boot_summary(boots)

# Example
N = 1000
x = np.random.normal([1, 100, 1000], [1,2,3], (N,3))
x = np.column_stack([np.ones(x.shape[0]), x])
y = np.matmul(x, [10,  10.5, 20.5, 40.5])
ols_boot(x, y, 1000)
#> {'coef': array([ 9.99999926, 10.5       , 20.5       , 40.5       ]),
#>  'se': array([1.47045649e-09, 1.70255501e-13, 3.86609646e-13, 1.36142024e-12])}