#Decide on the sample size, and the number of Monte Carlo repitions we would like to use:
n = 50
rep = 10000

#Specify the population model:
beta1 = 0.5
beta2 = 2
x = runif(n)

#Create vectors of zeros, that will later be used to store OLS estimates for the intercept and slope:
b1 = b2 = rep(0,n)

#Start the Monte Carlo loop:
for(j in 1:rep){

#Create the disturbances vector (epsilon) by randomly generating numbers from the N(0,1) distribution.
eps = rnorm(n)

#Now "draw" a random sample of y values:
y = beta1 + beta2*x + eps

#Calculate and record OLS estimates:
b1[j] = lm(y~x)$coefficient[1]
b2[j] = lm(y~x)$coefficient[2]

#End loop:
}

#Now, we can examine the simulated sampling distribution. For example:
hist(b1)

#See if b2 appears to be unbiased:
mean(b2)