data <- read.csv("hubble.csv", header = TRUE)     # replace by your own path here

Y <- data[,2]
x <- data[,1]
n <- length(Y)

lm1 <- lm(Y ~ x) 
summary(lm1)
 
plot(x,Y) 
  lines(x, predict(lm1))

## now let's do it "by hand"
 Xmat <- cbind(rep(1,n), x)    # "design matrix" 
 XTX1 <- solve(t(Xmat) %*% Xmat)
 betahat <- XTX1  %*% t(Xmat) %*% Y
 Yhat <- Xmat %*% betahat      # predicted values 

 shat <- sqrt( sum((Y - Yhat)^2/(n-2)) )           # estimate of the residual st.dev. 
 SE <- shat*sqrt(diag(XTX1))



 # examples of prediction and conf.int in R

x <- rnorm(15)
y <- x + rnorm(15)
predict(lm(y ~ x))
new <- data.frame(x = seq(-3, 3, 0.5))
predict(lm(y ~ x), new, se.fit = TRUE)
pred.w.plim <- predict(lm(y ~ x), new, interval="prediction")
pred.w.clim <- predict(lm(y ~ x), new, interval="confidence")
matplot(new$x,cbind(pred.w.clim, pred.w.plim[,-1]),
        lty=c(1,2,2,3,3), type="l", ylab="predicted y")