Run a linear regression under the classical OLS assumptions: well specification, homoscedacity and normal errors. Then draw a confidence box and a confidence ellipse.

classical_ols.py

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse, Rectangle
from matplotlib.transforms import Affine2D
from scipy.stats import t, f
from scipy.linalg import sqrtm
import statsmodels.api as sm



# generate some data that obeys a classical OLS assumptions
n_data= 100
conf_level=0.1
x = np.random.standard_normal(size=n_data)
X = np.c_[np.ones(n_data), x]
beta_circ = np.array([0.5, 0.2])
sigma2 = 0.7
e = np.random.standard_normal(size=n_data) * sigma2
y = X @ beta_circ + e
d = len(beta_circ)


# run the statsmodels regression for comparison
olsres=sm.OLS(y, X).fit()
print(olsres.summary(alpha=conf_level))

# run the regression under homoscedacity and well specification assumptions
beta_hat, _, _, _ = np.linalg.lstsq(X, y, rcond=-1)
e_hat = X @ beta_hat - y
s2 = e_hat @ e_hat / (n_data - d)  # unbiased estimate of the residual noise
var_hat = s2 * np.linalg.pinv(X.T @ X) # the variance of the estimator.
se_hat = np.sqrt(np.diag(var_hat))
tstats = beta_hat/ se_hat
fstat = beta_hat.T @ X.T @ X @ beta_hat / d / s2 

rss_restricted = (y-y.mean()) @ (y-y.mean())
rss_unrestricted = e_hat @ e_hat
n_restrictions = d-1
fstat2= (rss_restricted -rss_unrestricted) * (n_data-d) /n_restrictions / rss_unrestricted
print({'point estimates':beta_hat, 'standard error estiate':se_hat,'t stats':tstats,'f stat':fstat, 'f stat 2': fstat2})


# Plot the two confidence region style: box and ellipse
# both at level 10%
fig, ax = plt.subplots()
ax.scatter(beta_circ[0], beta_circ[1], label=r"$\beta_\circ$")
ax.scatter(beta_hat[0], beta_hat[1], label=r"$\hat \beta$", color="C1")
ax.set_aspect("equal")


# Make a confidence Box
n_std = t(df=n_data - d).ppf(q=1 - conf_level / 2)
rect = Rectangle((-1, -1), width=2, height=2, alpha=0.3,color="C1",label="Conf Box")
transf = (
    Affine2D().scale(se_hat[0] * n_std, se_hat[1] * n_std).translate(beta_hat[0], beta_hat[1])
)

rect.set_transform(transf + ax.transData)
ax.add_patch(rect)


# Make a confidence ellipse
tau = f(dfn=d, dfd=n_data - d).ppf(1 - conf_level)
Q = np.linalg.pinv(var_hat) / d / tau
ellipse = Ellipse(xy=(0, 0), width=2, height=2, alpha=0.3,color="C2",label="Conf Ell")
transf = Affine2D(matrix=np.block(
    [[sqrtm(np.linalg.pinv(Q)), beta_hat[:, np.newaxis]], [np.zeros((1, 2)), np.ones((1, 1))]]
))
ellipse.set_transform(transf + ax.transData)
ax.add_patch(ellipse)


ax.legend()
plt.show()