MC - Simple Linear Regression with Intercept

Model Building

$y={\beta }_0+{\beta }_{1}x+u$

Variable Declaration

y i = β 0 + β 1 x i + u i       { ( x i , y i ) , i = 1 , ⋯   , n } y_i=\beta_0+\beta_1x_i+u_i\ \ \ \ \ \{(x_i,y_i), i=1,\cdots,n\} yi​=β0​+β1​xi​+ui​     {(xi​,yi​),i=1,⋯,n} >

Variable (Parameter)	Meaning	Set Value (Distribution)
$x$	explanatory variable	$x\sim Normal(5,50)$
$u$	error term	$u\sim Normal(0,200)$
${\beta }_0$	intercept parameter	25
${\beta }_1$	slope parameter	3
$n$	simple size	10
$N$	frequency of simpling	10000

• Note:

1. There is almost no limit to the value of $X$ in OLS estimates, and in fact, we allow both $X$ and $Y$ to be $R.V.$
2. $R.V.$ $Y$ is generated by the following formula: $Y={\beta }_0+{\beta }_{1}X+u$

Results

${\hat{y}}_i={\hat{\beta }}_0+{\hat{\beta }}_1{x}_i,i=1,\cdots ,n$

Estimator	Mean Error	Standard Deviation of Error
${\hat{\beta }}_0$	0.68150	81.85889
${\hat{\beta }}_1$	-0.00912	1.31418

• Mean Abs. Error: 232.94118
• Mean Sq. Error: 40849.77162

Graphs

1. Graph1: Distribution of ${\hat{\beta }}_0$
2. Graph2: DIstribution of ${\hat{\beta }}_1$ (hint: There are some problems with the vertical data but they don’t affect the distribution.)

from IPython.display import Image
Image(filename='Graphs1.png')

from IPython.display import Image
Image(filename='Graphs2.png')

Coding

import numpy as np
import pandas as pd
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import statsmodels.api as sma
from scipy.stats import norm

# Params Setting
b0 = 25
b1 = 3
n = 10
N = 10000
x = np.random.normal(5, 50, n)
b0_set = np.zeros(N)
b1_set = np.zeros(N)
yab_set = np.zeros(N)
ysq_set = np.zeros(N)

# MC Simulation
i = 0
while i < N:
    u = np.random.normal(0, 200, n)
    y = b0 + b1 * x + u
    model = sma.OLS(y,sma.add_constant(x))
    results = model.fit()
    b0hat = results.params[0]
    b1hat = results.params[1]
    b0_set[i] = b0hat
    b1_set[i] = b1hat
    yhat = b0hat + b1hat * x
    yab_set[i] = np.abs(y,yhat).mean()
    ysq_set[i] = np.sum((y-yhat)**2) / n
    i +=1

# Results
dict_result = {}
dict_result['Mean Error'] = [b0_set.mean()-b0, b1_set.mean()-b1]
dict_result['Std. Error'] = [b0_set.std(), b1_set.std()]
df = pd.DataFrame(dict_result, index=['beta0_hat', 'beta1_hat'])
print(df)
print("Mean Abs. Error: {}".format(yab_set.mean()))
print("Mean Sq. Error: {}".format(ysq_set.mean()))

# Graphs
plt.figure(figsize=(25,15))
plt.style.use('seaborn')

# Graph1: Distribution of estimator beta0_hat
plt.subplot(221)
plt.title('Graph1: Distribution of estimator beta0_hat',fontsize=25, c='black')
plt.grid(True)
plt.xlabel('Value of beta0_hat',fontsize=15,c='black')
plt.ylabel('Probability',fontsize=15,c='black')
m0,bins0,patches0 = plt.hist(b0_set, bins=50, density=True, color='blue', edgecolor='black', linewidth=1, alpha=0.7)
B0_set = norm.pdf(bins0, b0_set.mean(), b0_set.std())
plt.plot(bins0,B0_set,'r--')

# Graph2: Distribution of estimator beta1_hat
plt.subplot(222)
plt.title('Graph2: Distribution of estimator beta1_hat',fontsize=25, c='black')
plt.grid(True)
plt.xlabel('Value of beta1_hat',fontsize=15,c='black')
plt.ylabel('Probability',fontsize=15,c='black')
m1,bins1,patches1 = plt.hist(b1_set, bins=40, density=True, color='r', edgecolor='black', linewidth=1, alpha=0.5)
B1_set = norm.pdf(bins1, b1_set.mean(), b1_set.std())
plt.plot(bins1,B1_set,'b-.')
plt.savefig('Graphs1&2.png', bbox_inches='tight')
plt.show()

原创文章，转载请注意出处。