I am trying to estimate 𝜷 = (𝛽􏰁, 𝛽􏰂)′ using OLS(matrix form) and store the values in a matrix of dimension 𝑟 × 2, with monte carlo simulation

Question

I do not quite understand how to store the values in a matrix of dimension 𝑟 × 2. This is how far i have come in R :

My regression : 𝑦=𝛽0+𝛽1𝑥+𝑢, Where B0= 1 and B1 = -1

set.seed(123)
n = 20
nreps=10000
beta_0 = -1 
beta_1 = 1 

ols = vector(mode="numeric", length=nreps)

##Start MC

for (r in 1:nreps) {
  u = rnorm(n, mean = 0, sd = 1) 
  x = rnorm(n, mean = 1, sd = 2^2) 
  y = beta_0 + beta_1*x + u
  
  head(x)

  beta_hat=solve(t(x) %*% x) %*% (t(x) %*% y)
  beta_hat 
  
  # print(paste("OLS estimate =", beta_hat))4f
  
  ols[r] = beta_hat
}
ols = cbind(ols)
ols_average=mean(ols[,1])
ols_sd=sd(ols[,1])
plot(density(ols)) 

Answer 1

To estimate the model's coefficients, use lsfit, the low-level function used by lm.

set.seed(307)
n <- 20
nreps <- 10000
beta_0 <- -1 
beta_1 <- 1 

beta_mat <- matrix(nrow = nreps, ncol = 2,
                   dimnames = list(NULL, c("Intercept", "x")))
for (r in seq.int(nreps)) {
  u <- rnorm(n, mean = 0, sd = 1) 
  x <- rnorm(n, mean = 1, sd = 2^2) 
  y <- beta_0 + beta_1*x + u
  fit <- lsfit(x, y)
  beta_mat[r, ] <- coef(fit)
}

beta_average <- colMeans(beta_mat)
beta_average
#Intercept         x 
#-1.005184  1.001375 


beta_sd <- apply(beta_mat, 2, sd)

dens <- density(beta_mat[,2])

hist(beta_mat[,2], freq = FALSE, ylim = c(0, 1.1*max(dens$y)),
     xlab = expression(beta[1]), main = expression("Histogram of" ~ beta[1]))
lines(dens)