## Use simulation to assess the coverage probabilities of confidence
## intervals for regression coefficients in a bivariate regression.

## Sample size.
n = 200

## Number of simulation replications.
nrep = 1e4

## Correlation between X1 and X2
r = 0.3

## Keep track of the coverage for each slope.
Q = c(0,0,0)

for (k in 1:nrep) {

  ## The correlation between the two covariates is r.
  X1 = rnorm(n)
  X2 = r*X1 + sqrt(1-r^2)*rnorm(n)

  ## The design matrix.
  X = cbind(array(1,n), X1, X2)

  ## The (X'*X)^(-1) matrix.
  qrd = qr(X)
  R = qr.R(qrd)
  H = solve(t(R) %*% R)

  ## Simulate data from the model.
  beta = rnorm(3)
  Y = X %*% beta + rnorm(n)

  ## Fit the model.
  m = lm(Y ~ X+0)

  ## The parameter estimates.
  sigma2_hat = sum(m$residuals^2)/(n-3)
  beta_hat = m$coefficients

  ## The standard errors for the slopes.
  SE = sqrt(diag(H)*sigma2_hat)
  
  ## The CI.
  LB = beta_hat - 2*SE
  UB = beta_hat + 2*SE

  Q = Q + (beta > LB) * (beta < UB)
}

Q = Q/nrep