# STAT3530 -- the sampling distribution of the least squares estimator
# Simulate many datasets from a model we choose, fit each one, and see how the
# estimates behave. Run a line at a time.

# --- the model we simulate from -------------------------------------------
set.seed(30925)

n <- 25
beta <- c(10, 2)       # the true intercept and slope
sigma <- 8             # the true error standard deviation

x <- seq(1, 20, length = n)   # predictors: we choose these, they are not random
# x <- x - mean(x)              # center the predictor (leave uncommented until end)
X <- cbind(1, x)              # predictor matrix

# simulate one dataset
set.seed(83026)
y <- X %*% beta + rnorm(n, mean = 0, sd = sigma)
plot(x, y, ylab = 'y')
abline(coef = beta, lty = 2)
abline(lm(y ~ x), col = 'blue')

# --- simulate many datasets, and store the estimates from each -------------
nsim <- 2000
B <- matrix(NA, nsim, 2)
colnames(B) <- c('intercept', 'slope')

for (i in 1:nsim) {
  y <- X %*% beta + rnorm(n, mean = 0, sd = sigma)   # a new dataset
  B[i, ] <- coef(lm(y ~ x))                          # keep its estimates
}

head(B)                # one row per simulated dataset

# --- plot some example estimates ------------------------------------------

# what happens as you plot more lines?
n_fits <- 10

# plot
for(j in 1:n_fits){
  if(j == 1){
    plot(x, X %*% (B[j, ]), type = 'l', 
         ylab = 'y', col = rgb(0, 0, 1, 0.25))
  }else{
    lines(x, X %*% (B[j, ]), col = rgb(0, 0, 1, 0.25))
  }
  lines(x, X %*% beta, type = 'l', lwd = 2, col = 'red', ylab = 'y')
}
  
# --- summarise the simulated estimates ------------------------------------
colMeans(B)                    # one mean per coordinate
apply(B, 2, var)               # one variance per coordinate
cov(B[, 1], B[, 2])            # covariance between coordinates
cor(B[, 1], B[, 2])            # correlation between coordinates

# --- what the theory says --------------------------------------------------
# var(betahat) = sigma^2 (X'X)^{-1}
V <- sigma^2 * solve(t(X) %*% X)
V

# --- compare ---------------------------------------------------------------
rbind(simulated = colMeans(B),      theory = beta)
rbind(simulated = apply(B, 2, var), theory = diag(V))

c(simulated = cov(B[, 1], B[, 2]),  theory = V[1, 2])
c(simulated = cor(B[, 1], B[, 2]),  theory = V[1, 2] / sqrt(V[1, 1] * V[2, 2]))

# --- what does the sampling distribution look like? ------------------------
plot(B, xlab = 'intercept estimate', ylab = 'slope estimate')
points(beta[1], beta[2], col = 'red', pch = 19, cex = 2)   # the truth

# --- your turn -------------------------------------------------------------
# 1. The two estimates are strongly negatively correlated. Look back at the
#    scatterplot and explain why, in terms of how a line can tilt.
#
# 2. Re-run everything with x replaced by x - mean(x) by uncommenting line 13. 
#    What happens to the correlation, and can you see why?