# This program reproduces the figures and output of 
# Cod Catch data of Example 6.6 (Table 6.4) in Bowerman, O'Connell and Koehler
# Created by C. Andy Tsao @ NDHU, 060529

# Change the working directory to where the data is stored
# Load packages 
#      foreign for reading Splus data
#      car for DW test
library(foreign);library(car);

ds<-read.S("t6-4 hotel.splus")
dt<-data.frame(y=ds[,1],t=as.integer(row.names(ds)))
dt; attach(dt)

# Time series plot
# Transform ys<- sqrt(y)
ys<-sqrt(y); yqr<-y^(0.25); logy<-log(y);
par(mfrow=c(4,1));
plot(t,y,main="Monthly Hotel Room Average",xlab="Time",ylab="y");
plot(t,ys,main="Monthly Hotel Room Average",xlab="Time",ylab=expression(sqrt(y)));
plot(t,yqr,main="Monthly Hotel Room Average",xlab="Time",ylab=expression(y^(0.25)));
plot(t,logy,main="Monthly Hotel Room Average",xlab="Time",ylab=expression(log(y)));
# Follow the recomendation in BOK, we use log(y) as the new response hereafter

# Creating sin(2Pi*t/L), cos(2Pi*t/L) (conti) variables
L<-12;
s2t<-sin(2*pi*t/L); c2t<-cos(2*pi*t/L);
s4t<-sin(4*pi*t/L); c4t<-cos(4*pi*t/L);

# Creating Month, Season/Quater (categorical) variables 
mconvert <- function(t) # month converter 
{  v<-mat.or.vec(NROW(t),2)

tt<-t%%12; v[,1]<-tt;

tt<-tt+(tt==0)*12; tt<-tt-.01; tt<-floor(tt/3)
v[,2]<-(tt+1)%%4  
list(month=v[,1],qtr=v[,2]); # returing month, quarters
}

mct<-mconvert(t);
M<-as.factor(mct$month); Q<-as.factor(mct$qtr); 

# Polynomial regression of 1st and 2nd order
m.month<-lm(logy~t+M); anova(m.month);summary(m.month);
durbin.watson(m.month);

m.trig2<-lm(logy~t+s2t+c2t);anova(m.trig2);summary(m.trig2);
durbin.watson(m.trig2);

m.trig4<-lm(logy~t+s2t+c2t+s4t+c4t);anova(m.trig4);summary(m.trig4);
durbin.watson(m.trig4);

m.season<-lm(logy~t+Q);anova(m.season);summary(m.season);
durbin.watson(m.trig4);

nt<-169:180; 
ns2t<-sin(2*pi*nt/L); nc2t<-cos(2*pi*nt/L);
ns4t<-sin(4*pi*nt/L); nc4t<-cos(4*pi*nt/L);

new<- data.frame(t=nt,s2t=ns2t,c2t=nc2t,s4t=ns4t,c4t=nc4t)
p.mt4<-predict(m.trig4, new, se.fit = TRUE);
pred.value<-p.mt4$fit; se.mpred<-p.mt4$se
pred.w.plim <- predict(m.trig4, new, interval="prediction")
pred.w.clim <- predict(m.trig4, new, interval="confidence")
cbind(nt,pred.value,se.mpred,pred.w.clim[,-1],pred.w.plim[,-1])

# Overlay CI and PI at the given t values
windows()
matplot(new$t,cbind(pred.w.clim, pred.w.plim[,-1]),
        lty=c(1,2,2,3,3), type="l", ylab="predicted y")