Thanks to atireto's comment, I realized that I haven't taken into account the covariances between the estimated parameters. I was assuming that a and b are uncorrelated random variables, but they are not. So I corrected the sampling of a and b by converting them to correlated random variables using Cholesky decomposition of variance-covariance matrix (see here).
Here's the corrected code:
rm(list = ls())
n = 20 # number of points
N = 10000 # simulation rounds
x = 1:n # x
s = 200 # standard deviation of error
set.seed(12)
e = rnorm(n, 0, s)
a = 5; b = 5 # slope, intercept
y = a * x + b + e
## lm model
mod1 = lm(y ~ x)
summ_mod1 = summary(mod1)
## std errors
coefs = summ_mod1$coefficients[, 1:2]
sigma = summ_mod1$sigma
## predict new data using simulation
new.x = (min(x) - 2 * max(x)):(max(x) + 2 * max(x))
upper = rep(0, length(new.x))
lower = rep(0, length(new.x))
tmp = rep(0, N)
for (i in 1:length(new.x)) {
tmp = rnorm(n = N, mean = coefs[1, 1], sd = coefs[1, 2]) +
rnorm(n = N, mean = coefs[2, 1], sd = coefs[2, 2]) * new.x[i] +
rnorm(n = N, mean = 0, sd = sigma)
upper[i] = mean(tmp) + 1.96 * sd(tmp)
lower[i] = mean(tmp) - 1.96 * sd(tmp)
}
plot(x, y, type = "p",
xlim = c(min(new.x), max(new.x)),
ylim = c(min(lower), max(upper)))
lines(new.x, lower, col = "blue")
lines(new.x, upper, col = "blue")
## R's prediction interval
pred.int = predict(object = mod1,
newdata = data.frame(x = new.x),
interval = "predict",
level = 0.95)
pred.lower = pred.int[,2]
pred.upper = pred.int[,3]
lines(new.x, pred.lower, col = "red")
lines(new.x, pred.upper, col = "red")
## taking the variance-covariance matrix into account
L = chol(vcov(mod1))
for (i in 1:length(new.x)) {
beta = matrix(c(rnorm(n = N),
rnorm(n = N)),
nrow = 2)
beta_cor = t(L) %*% beta
tmp = coefs[1, 1] + beta_cor[1, ] +
(coefs[2, 1] + beta_cor[2, ]) * new.x[i] +
rnorm(n = N, mean = 0, sd = sigma)
upper[i] = mean(tmp) + 1.96 * sd(tmp)
lower[i] = mean(tmp) - 1.96 * sd(tmp)
}
lines(new.x, lower, col = "green")
lines(new.x, upper, col = "green")
legend(x = -40, y = 1600,
legend = c("sim - without cor", "sim - with co", "predict.lm"),
lty = c(1, 1, 1),
col = c("blue", "green", "red"))
It still doesn't perfectly match the R's predict.lm result, but it's very close. In fact, if I change the 1.96 coefficient to something like 2.1, it matches it perfectly. Predict.lm is supposed to be 95% prediction interval, so I'm not sure this is happening.