compare_ols.py

import statsmodels.api as sm
import numpy as np
from sklearn.linear_model import LinearRegression

# get original NumPy implementation
β = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y).flatten()

# get sklearn's LinearRegression implementation weights
lr = LinearRegression()
sklearn_coefficients = lr.fit(X,y).coef_.flatten()

# get statsmodel implementation weights
result = sm.OLS(y, X).fit()
sm_coefficients = result.params.values

# zip up all the results into one consolidated dataframe to inspect
coefficients_df = pd.DataFrame(list(zip(features, sklearn_coefficients, sm_coefficients, β)), 
                               columns=["feature", "sklearn", "statsmodel", "numpy"])