---
title: "Prediction for the Rocket Propulsion Data"
output: pdf_document
---

```{r}
data.source<-"http://www.math.mcgill.ca/yyang/regression/data/2-1-RocketProp.csv"
RocketProp<-read.csv(file=data.source)
names(RocketProp)<-c('i','Strength','Age')
x<-RocketProp$Age
y<-RocketProp$Strength
n<-length(x)
plot(x,y,pch=19,xlab='Age',ylab='Shear Strength')
fit.RP<-lm(y~x)
summary(fit.RP)
abline(coef(fit.RP),col='red')
title('Prediction Rocket Propulsion Data')
```


Predictions on a grid of new x values using the \texttt{predict} function:



```{r}
xnew<-seq(0,25,by=0.1)

#Confidence interval
ynew.interval<-predict(fit.RP,newdata=data.frame(x=xnew),interval='confidence')
plot(x,y,pch=19,xlab='Age',ylab='Shear Strength')
abline(coef(fit.RP),col='red')
title('Prediction for Rocket Propulsion Data')
lines(xnew,ynew.interval[,2],lty=2,col='red')
lines(xnew,ynew.interval[,3],lty=2,col='red')

#Prediction interval
yonew.interval<-predict(fit.RP,newdata=data.frame(x=xnew),interval='prediction')
lines(xnew,yonew.interval[,2],lty=2,col='blue')
lines(xnew,yonew.interval[,3],lty=2,col='blue')
legend(20,2600,c('Prediction','Conf. Interv.','Pred. Interv.'),col=c('red','red','blue'),
       lty=c(1,2,2))
```



Using the matrix formulae:



```{r}
X<-cbind(rep(1,n),x)
(XtX<-t(X)%*%X)
Xty<-t(X) %*% y
beta.hat<-solve(XtX,Xty)
Hat.matrix<-X %*% solve(XtX) %*% t(X)

SS.Res<-(t(y)%*%(diag(1,n)-Hat.matrix)%*%y)[1,1]
MS.Res<-SS.Res/(n-2)
sigma.hat<-sqrt(MS.Res)

nnew<-length(xnew)
Xnew<-cbind(rep(1,nnew),xnew)

Ynew<-Xnew%*%beta.hat[,1]
Hnew<-diag(Xnew %*% solve(XtX) %*% t(Xnew))
Var.Ynew<-sigma.hat^2*Hnew
Var.Yonew<-sigma.hat^2*(1+Hnew)

#Look up the quantile
alpha<-0.05
q.val<-qt(1-alpha/2,n-2)

#Intervals
Conf.vals<-cbind(Ynew-q.val*sqrt(Var.Ynew),Ynew+q.val*sqrt(Var.Ynew))
Pred.vals<-cbind(Ynew-q.val*sqrt(Var.Yonew),Ynew+q.val*sqrt(Var.Yonew))

#Check computation - should be zero
sum(ynew.interval[,2:3]-Conf.vals)
sum(yonew.interval[,2:3]-Pred.vals)
```