sim.R

require('mvtnorm')

set.seed(1001)

rho = 0.25

sigma = matrix(0, 500, 500)
for (i in 1:500) {
  for (j in 1:500) {
    sigma[i, j] = rho^(abs(i - j))
  }
}

X = rmvnorm(n = 150, mean = rep(0, 500), sigma = sigma)

c1 = 3.5  # non-zero coefficients
beta1 = matrix(c(rep(c1, 5), rep(0, 495)))

sd1 = 1
eps = matrix(rnorm(150, 0, sd1))

y = (X %*% beta1) + eps