import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import scipy.stats as stats

# set the random seed:
np.random.seed(1234567)

# set sample size and number of simulations
n = 1000
r = 10000

# set true parameters:
beta0 = 1
beta1 = 0.5
su = 2
sx = 1
ex = 4

# initialize b0 and b1 to store results later:
b0 = np.empty(r)
b1 = np.empty(r)

# repeat r times:
for i in range(r):
 # draw a sample:
 x = stats.norm.rvs(ex, sx, size=n)
 u = stats.norm.rvs(0, su, size=n)
 y = beta0 + beta1 * x + u
 df = pd.DataFrame({'y': y, 'x': x})

 # estimate OLS:
 reg = smf.ols(formula='y ~ x', data=df)
 results = reg.fit()
 b0[i] = results.params['Intercept']
 b1[i] = results.params['x']