library(bbmle)

par(mfrow=c(2,1))

set.seed(123456)
n<-500; r<-1000
b0<-1; b1<-0.5; su<-2


b0ols<-numeric(r)
b1ols<-numeric(r)

for(j in 1:r) {
  x<-rnorm(n,4,1)
  u<-rnorm(n,0,su)
  y<-b0+b1*x+u
  
  bols<-coefficients(lm(y~x))
  b0ols[j]<-bols["(Intercept)"]
  b1ols[j]<-bols["x"]
  }

hist(b1ols, breaks=100, xlim=c(0,1), xlab="n=500")


b0mle<-numeric(r)
b1mle<-numeric(r)

for(j in 1:r) {

  x<-rnorm(n,4,1)
  u<-rnorm(n,0,su)
  y<-b0+b1*x+u
  
fn<-function(beta0, beta1, sigma) {
  (n/2)*log(sigma^2)+1/(2*sigma^2)*(sum((y-beta0-beta1*x)^2))
  }
res<-mle2(fn,start=list(beta0=0.8, beta1=0.3, sigma=1.5))
summary(res)
str(summary(res))

b0mle[j]<-summary(res)@coef["beta0","Estimate"]
b1mle[j]<-summary(res)@coef["beta1","Estimate"]
}

hist(b1mle, breaks=100, xlim=c(0,1),xlab="n=500")

mean(b0ols)
mean(b1ols)

var(b0ols)
var(b1ols)

mean(b0mle)
mean(b1mle)

var(b0mle)
var(b1mle)