setwd("D://R Scripts/Linear Models/Data Files")

library(nlme)
library(lmtest)

par(mfrow = c(2., 2.), pty = "s", oma = c(0., 0., 4., 0.))
confood<-read.table("Confood.dat",header=T)
attach(confood)
Promotion=as.factor(Promotion)
plot(log(Price),log(Sales),type="n")
text(log(Price),log(Sales),c("+","*")[Promotion])
plot(Week,log(Sales),type="l")
text(Week,log(Sales),c("+","*")[Promotion],col=c("red","blue")[Promotion])
plot(log(SalesLag1),log(Sales),xlab=expression(log(Sales)[t-1]),ylab=expression(log(Sales)[t]),type="n")
text(log(Sales)[-1],log(Sales),c("+","*")[Promotion])
#acf(log(Sales))

#------------------------------------------------------------------------
#            LINEAR MODEL (OLS)
#-----------------------------------------------------------------------


par(mfrow=c(2,2))
sales.lm=lm(log(Sales)~log(Price)+Week+Promotion)
print(summary(sales.lm))
res.st=rstandard(sales.lm)
n=length(res.st)
plot(sales.lm,which=c(1,2))
#boxcox(sales.lm)
print(dwtest(sales.lm,alternative=c("two.sided")))
plot(res.st[-1],res.st[-n],xlab=expression(res.st[t-1]),ylab=expression(res.st[t]))
acf(res.st,type="partial")

#---------------------------------------------------------------------------
#           LINEAR MODEL WITH AR(1) ERRORS (GLS)
#---------------------------------------------------------------------------

sales.gls=gls(log(Sales)~log(Price)+Week+Promotion,correlation=corAR1(form=~Week),method="ML")
par(mfrow=c(1,1))
print(summary(sales.gls),Correlation=F)
plot(Week,log(Sales),type="n")
lines(Week,fitted(sales.gls),lty=1,col="blue")
lines(Week,fitted(sales.lm),lty=1,col="red")
text(Week,log(Sales),c("+","*")[Promotion],col=c("red","blue")[Promotion])

var=diag(model.matrix(sales.lm)%*%sales.gls$varBeta%*%t(model.matrix(sales.lm)))
lines(Week,fitted(sales.gls)+qnorm(.95)*sqrt(var),lty=2,col="blue")
lines(Week,fitted(sales.gls)-qnorm(.95)*sqrt(var),lty=2,col="blue")
lines(Week,predict(sales.lm,interval="confidence")[,2],lty=2,col="red")
lines(Week,predict(sales.lm,interval="confidence")[,3],lty=2,col="red")

legend(17,8.8,legend=c("gls fitted","lm fitted", "gls CI", "lmI"),lty=c(1,1,2,2),col=c("blue","red","blue","red"))