data.source<-"http://www.math.mcgill.ca/dstephens/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')
summary(fit.RP)
abline(coef(fit.RP),col='red')
title('Prediction Rocket Propulsion Data')


#Predictions on a grid of x* values.

xstar<-seq(0,25,by=0.1)

#Confidence interval for prediction
ystar.interval<-predict(fit.RP,newdata=data.frame(x=xstar),interval='confidence')
lines(xstar,ystar.interval[,2],lty=2,col='red')
lines(xstar,ystar.interval[,3],lty=2,col='red')

#Prediction interval for prediction
yostar.interval<-predict(fit.RP,newdata=data.frame(x=xstar),interval='prediction')
lines(xstar,yostar.interval[,2],lty=3,col='blue')
lines(xstar,yostar.interval[,3],lty=3,col='blue')

legend(15,2600,c('Prediction','Conf. Interv.','Pred. Interv.'),col=c('red','red','blue'),lty=c(1,2,3))

dev.copy2pdf(file='RocketPropPredict.pdf',paper='USr',width=11,height=9)

############################################
#Using the matrix formulae

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)

nstar<-length(xstar)
Xstar<-cbind(rep(1,nstar),xstar)

Ystar<-Xstar%*%beta.hat[,1]
hstar<-diag(Xstar %*% solve(XtX) %*% t(Xstar))
Var.Ystar<-sigma.hat^2*hstar
Var.Yostar<-sigma.hat^2*(1+hstar)

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

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

#Check
ystar.interval[,2:3]-Conf.vals
yostar.interval[,2:3]-Pred.vals