<jupyter_start><jupyter_text>Simulation of variance of OLS parametersRecently, I've spend some time to understand what *variance* of parameters actually means. If I understand it correctly, the variance of an estimate refers to how much we expect our estimate ($\hat{\beta}$) to vary from the true population parameter ($\beta$). In OLS, two terms contribute to the variance estimate: the unexplained variance and the design variance. Unexplained variance, also referred to as $\hat{\sigma}^{2}$, is defined as:\begin{align}\hat{\sigma}^{2} = \frac{\sum_{i=1}^{N} (y_{i} - \hat{y}_{i})^{2}}{N - K}\end{align}The design variance, for any contrast $c$, is defined as:\begin{align}desvar = c(X'X)^{-1}c'\end{align}As such, the variance of any contrast of beta-parameter(s) is defined as:\begin{align}\mathrm{var}[c\beta] = \frac{\sum_{i=1}^{N} (y_{i} - \hat{y}_{i})^{2}}{N - K} c(X'X)^{-1}c'\end{align}While mathematically, I understand *how* variance is computed. But it's conceptual meaning -- i.e. the expected squared deviation from the true population parameter -- is still a bit muddy for me. As such, I decided to simulate the process, which often helps me to understand statistical issues.First, some imports.<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Now, let's define an arbitrary *true* model of some dependent variable $y$. Let's assume the true function is:\begin{align}y = 10 + 5x\end{align}<jupyter_code>x = np.arange(100)
y = 10 + 5 * x
plt.plot(x, y)
plt.title("True function of y")
plt.show()<jupyter_output><empty_output><jupyter_text>Let's now add some random gaussian noise ($\mu = 0, \sigma=1$) to our dependent variable and estimate the parameter of the intercept and $x$ using linear regression. We'll also calculate the variance of this estimate (using the contrast-vector `[1, 0]` for the intercept-parameter and `[0, 1]` for the $x$-parameter). We'll iterate this process 10,000 times and keep track of the parameter-estimates and their variance-estimates.<jupyter_code>iters = 50000
betas = np.zeros((iters, 2))
vars = np.zeros((iters, 2))
c = np.array([[1, 0], [0, 1]])
X = np.hstack((np.ones((100, 1)), x[:, np.newaxis]))
    
for i in range(iters):
    y2 = y + np.random.normal(0, 1, 100) * 100
    betas[i, :] = np.linalg.lstsq(X, y2)[0]
    y_hat = X.dot(betas[-1, :])
    sigma_hat = np.sum((y_hat - y2) ** 2) / (X.shape[0] - X.shape[1])
    desvar = c.dot(np.linalg.pinv(X.T.dot(X))).dot(c.T)
    vars[i, :] = np.diag(sigma_hat * desvar)<jupyter_output><empty_output><jupyter_text>Let's plot the histograms of the two parameter-estimates:<jupyter_code>plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.hist(betas[:, 0], bins=50)
plt.title(r"$\hat{\beta}_{0}$", fontsize=20)

plt.subplot(1, 2, 2)
plt.hist(betas[:, 1], bins=50)
plt.title(r"$\hat{\beta}_{x}$", fontsize=20)

plt.show()<jupyter_output><empty_output><jupyter_text>Now, the central limit theorem (CTL) predicts that the mean of this distribution should approximate the true parameters (I think). Let's check that:<jupyter_code>betas.mean(axis=0)<jupyter_output><empty_output><jupyter_text>Now, let's calculate the variance of our parameter-estimates across the simulations:<jupyter_code>betas.var(axis=0)<jupyter_output><empty_output><jupyter_text>This estimate from our simulation should then (I think) approximate the average variance estimate across simulations:<jupyter_code>vars.mean(axis=0)<jupyter_output><empty_output><jupyter_text>... which seems indeed the case!<jupyter_code>plt.figure(figsize=(15, 5))
plt.subplot(1, 2, 1)
plt.hist(vars[:, 0], bins=50)
plt.title(r"$\hat{\mathrm{var}[\beta}_{0}]$", fontsize=20)

plt.subplot(1, 2, 2)
plt.hist(vars[:, 1], bins=50)
plt.title(r"$\hat{\mathrm{var}[\beta}_{x}]$", fontsize=20)

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15759460860549507, 'fraction_numerical': 0.026438569206842923, 'mean_word_length': 5.688041594454073, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '191919', 'license': 'MIT', 'n_tokens_mistral': 1445, 'n_tokens_neox': 1327, 'n_words': 387, 'path': 'simulation_variance_OLS_parameters.ipynb', 'repo_name': 'lukassnoek/random_notebooks', 'size': '3858'}