Last active
March 21, 2022 11:06
-
-
Save el-hult/ff429527771d649d95fd0e629db585bb to your computer and use it in GitHub Desktop.
Run a linear regression under the classical OLS assumptions: well specification, homoscedacity and normal errors. Then draw a confidence box and a confidence ellipse.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment