library(bbmle) par(mfrow=c(3,1)) set.seed(123456) n<-5; r<-1000 b0<-1; b1<-0.5; su<-2 b0ols5<-numeric(r) b1ols5<-numeric(r) for(j in 1:r) { x<-rnorm(n,4,1) u<-rnorm(n,0,su) y<-b0+b1*x+u bols5<-coefficients(lm(y~x)) b0ols5[j]<-bols5["(Intercept)"] b1ols5[j]<-bols5["x"] } n<-50; r<-1000 b0<-1; b1<-0.5; su<-2 b0ols50<-numeric(r) b1ols50<-numeric(r) for(j in 1:r) { x<-rnorm(n,4,1) u<-rnorm(n,0,su) y<-b0+b1*x+u bols50<-coefficients(lm(y~x)) b0ols50[j]<-bols50["(Intercept)"] b1ols50[j]<-bols50["x"] } n<-500; r<-1000 b0<-1; b1<-0.5; su<-2 b0ols500<-numeric(r) b1ols500<-numeric(r) for(j in 1:r) { x<-rnorm(n,4,1) u<-rnorm(n,0,su) y<-b0+b1*x+u bols500<-coefficients(lm(y~x)) b0ols500[j]<-bols500["(Intercept)"] b1ols500[j]<-bols500["x"] } hist(b1ols5, breaks=100, xlim=c(-2,2), xlab="n=5") hist(b1ols50, breaks=100, xlim=c(-2,2), xlab="n=50") hist(b1ols500, breaks=100, xlim=c(-2,2), xlab="n=500") n<-5; r<-1000 b0<-1; b1<-0.5; su<-2 b0mle5<-numeric(r) b1mle5<-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)) b0mle5[j]<-summary(res)@coef["beta0","Estimate"] b1mle5[j]<-summary(res)@coef["beta1","Estimate"] } n<-50; r<-1000 b0<-1; b1<-0.5; su<-2 b0mle50<-numeric(r) b1mle50<-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)) b0mle50[j]<-summary(res)@coef["beta0","Estimate"] b1mle50[j]<-summary(res)@coef["beta1","Estimate"] } n<-500; r<-1000 b0<-1; b1<-0.5; su<-2 b0mle500<-numeric(r) b1mle500<-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)) b0mle500[j]<-summary(res)@coef["beta0","Estimate"] b1mle500[j]<-summary(res)@coef["beta1","Estimate"] } hist(b1mle5, breaks=100, xlim=c(-2,2),xlab="n=5") hist(b1mle50, breaks=100, xlim=c(-2,2),xlab="n=50") hist(b1mle500, breaks=100, xlim=c(-2,2),xlab="n=500") mean(b0ols5) mean(b1ols5) mean(b0ols50) mean(b1ols50) mean(b0ols500) mean(b1ols500) var(b0ols5) var(b1ols5) var(b0ols50) var(b1ols50) var(b0ols500) var(b1ols500) mean(b0mle5) mean(b1mle5) mean(b0mle50) mean(b1mle50) mean(b0mle500) mean(b1mle500) var(b0mle5) var(b1mle5) var(b0mle50) var(b1mle50) var(b0mle500) var(b1mle500)