rm(list=ls())

#### Compute & plot Prediction and confidence bands for WHO data

D = read.csv("WHO_Life_expectancy.csv",header=TRUE)
D = D[D$Year == 2000,]
D = D[,c(1,4,7,8,17,22)]
D = na.omit(D)

Country = D[,1]
D = D[,-1]
attach(D)
n = nrow(D)
print(n)

## lm Log Covariates

LogPE = log(percentage.expenditure+1)
LogGDP = log(GDP)
LogSchooling = log(Schooling+1)
Dlog = data.frame(Life.expectancy=Life.expectancy,Alcohol=Alcohol,LogPE=LogPE,
                  LogGDP=LogGDP, LogSchooling=LogSchooling)
out = lm(Life.expectancy ~ Alcohol + LogPE +  LogGDP + LogSchooling,data=Dlog)
## print(summary(out))
## print(confint(out))

### PREDICTIONS 
### Prediction Band for hat(Y)

pdf("predictY.pdf")

m = 100  
newD = data.frame(Alcohol=rep(4.5,m),LogPE=rep(4.41,m),LogGDP=rep(6.9,m),
                  LogSchooling=seq(-6,10,length=m))
attach(newD)
tmp = predict(out,newdata=newD,interval="prediction",level=.95)
LS = seq(-6,10,length=m)
plot(LS,tmp[,1],lwd=3,ylim=c(min(tmp),max(tmp)),type="l",xlab="",
     ylab="",col="purple",cex.axis=1.5)    
polygon(c(rev(LS),LS),c(rev(tmp[,3]),tmp[,2]),density=NULL,angle=90,  
        border=NA,col="thistle")     
lines(LS,tmp[,1],col="purple",lty=4,lwd=4)
lines(LS,tmp[,2],col="purple",lwd=2)  
lines(LS,tmp[,3],col="purple",lwd=2)
title(ylab=expression(hat(Y)), line=1.9, cex.lab=2.1)
title(xlab="Log Schooling", line=2.5, cex.lab=2.2)
dev.off()

#### Confidence bands for hat(mu)

pdf("predictMU.pdf")

tmp2 = predict(out,newdata=newD,interval="confidence",level=.95)
plot(LS,tmp2[,1],lwd=4,col="purple",ylim=c(min(tmp2),max(tmp2)),type="l",
     xlab="",ylab="",cex.axis=1.5) 
polygon(c(rev(LS),LS),c(rev(tmp2[,3]),tmp2[,2]),density=NULL,angle=90, 
        border=NA,col="thistle")
lines(LS,tmp2[,1],col="purple",lty=4,lwd=3)
lines(LS,tmp2[,2],col="purple",lwd=1.5)
lines(LS,tmp2[,3],col="purple",lwd=1.5)
title(ylab=expression(mu(x)), line=2.1,cex.lab=1.9)
title(xlab="Log Schooling", line=2.5, cex.lab=2)
dev.off()