<jupyter_start><jupyter_text>Different ways to compute OLS Import all necessary modules<jupyter_code>import numpy as np
import matplotlib.pylab as plt
import seaborn as sns

from numpy.linalg import inv, lstsq

sns.set_context('notebook')<jupyter_output><empty_output><jupyter_text>If you want to embed plots inside IPython notebook, you need to turn on this option.<jupyter_code>%matplotlib inline<jupyter_output><empty_output><jupyter_text>Generate the data Define the number of observations in the sample, $N$, and the number of simulations, $S$.<jupyter_code>N, S = 100, 1000<jupyter_output><empty_output><jupyter_text>Parameters of the joint distribution of $X$ and $e$.<jupyter_code>mean = [0,0]
rho = .1
cov = [[1, rho], [rho, 1]]<jupyter_output><empty_output><jupyter_text>Regression level and slope.<jupyter_code>alpha, beta = 2, 3<jupyter_output><empty_output><jupyter_text>Simulate $X$ and $e$ from multivariate normal distribution.<jupyter_code>def simulate_data(mean, cov, alpha, beta, size):
    X, e = np.random.multivariate_normal(mean, cov, size).T
    Y = alpha + beta * X + e
    return Y, X

Y, X = simulate_data(mean, cov, alpha, beta, (N, S))<jupyter_output><empty_output><jupyter_text>OLS EstimatorDefine two functions which should return the same slope parameter. The difference is in the implementation. The first one is a bit more natural since it uses NumPy matrices. This allows to use `*` symbol to multiply matrices. the second uses NumPy arrays instead. This implementation requires the use of `dot` function.<jupyter_code>def ols_matrix(Y,X):
    Y = np.matrix(Y).T
    X = np.matrix(np.vstack((np.ones_like(X), X))).T
    beta = np.array(inv(X.T * X) * (X.T * Y))
    return float(beta[1])

def ols_array(Y,X):
    X = np.vstack((np.ones_like(X), X)).T
    Qxx = np.dot(X.T, X)
    Qxy = np.dot(X.T, Y)
    beta = np.dot(inv(Qxx), Qxy)
    return float(beta[1])

def ols_lstsq(Y,X):
    X = np.vstack((np.ones_like(X), X)).T
    beta = lstsq(X,Y)[0]
    return float(beta[1])<jupyter_output><empty_output><jupyter_text>Compute the estimates Now we need to initialize two containers for our estimates in each simulated sample.<jupyter_code>beta_hat1, beta_hat2, beta_hat3 = [], [], []<jupyter_output><empty_output><jupyter_text>For each simulated sample we compute the estimate and append it to the list. IPython "magic" `%%timeit` will show how much time was spent in the current block.<jupyter_code>%%time
for y, x in zip(Y, X):
    beta_hat1.append(ols_matrix(y, x))<jupyter_output>CPU times: user 484 ms, sys: 16 ms, total: 500 ms
Wall time: 254 ms<jupyter_text>Same for the second implimentation through arrays instead of matrices.<jupyter_code>%%time
for y, x in zip(Y, X):
    beta_hat2.append(ols_array(y, x))
%%time
for y, x in zip(Y, X):
    beta_hat3.append(ols_lstsq(y, x))<jupyter_output>CPU times: user 436 ms, sys: 8 ms, total: 444 ms
Wall time: 233 ms<jupyter_text>Analyze the results Are the results equal?<jupyter_code>print(np.array_equal(beta_hat1, beta_hat2))
print(np.array_equal(beta_hat1, beta_hat3))
print(np.allclose(beta_hat1, beta_hat3, atol=1e-20))<jupyter_output>True
False
True<jupyter_text>Plot the histogram.<jupyter_code>plt.figure(figsize = (12, 6))
plt.hist(beta_hat1, 50, histtype='stepfilled', normed=True, lw=0, alpha=.5, label='Density')
plt.axvline(beta, color='red', lw=5, label='True')
plt.axvline(np.mean(beta_hat1), color='black', lw=5, label='Mean Estimate')
plt.xlabel(r'$\hat{\beta}$')
plt.ylabel('%')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Monte Carlo experimentRun estimation for all sample sizes and all simulated data.<jupyter_code>nobs = [250, 500, 1000, 2000]

beta_hat = []
for N in nobs:
    Y, X = simulate_data(mean, cov, alpha, beta, (N, S))
    temp = []
    for y, x in zip(Y, X):
        temp.append(ols_array(y, x))
    beta_hat.append(temp)
    
beta_hat = np.array(beta_hat)<jupyter_output><empty_output><jupyter_text>Plot densities of the estimates.<jupyter_code>plt.figure(figsize=(8, 5))
for i in range(len(nobs)):
    sns.kdeplot(beta_hat[i], alpha=.4, lw=3, shade=True, label=nobs[i])
plt.axvline(beta, color='red', lw=3, label='True')
plt.xlabel(r'$\hat{\beta}$')
plt.ylabel('%')
plt.legend(title='N')
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1380762930025743, 'fraction_numerical': 0.019658319681722443, 'mean_word_length': 5.7095761381475665, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 54, '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': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '377075', 'license': 'MIT', 'n_tokens_mistral': 1687, 'n_tokens_neox': 1601, 'n_words': 414, 'path': 'notebooks/computing_ols.ipynb', 'repo_name': 'khrapovs/metrix', 'size': '4273'}