# This program reproduces the figures and output of
# Cod Catch data of Example 9.1 (Table 9.1) in Bowerman, O'Connell and Koehler
# for stational and non-stationary ts
# Created by C. Andy Tsao @ NDHU, 060604

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

ds<-read.spss("t9-1 towel.sav",to.data.frame=TRUE) #Input SPSS data
dt<-data.frame(y=ds$Y,t=as.integer(row.names(ds)))
dt; attach(dt)

# Time series plot
yt<-as.ts(y);        # Create yt as a ts object
zt1<-diff(yt,differences=1);
zt2<-diff(yt,differences=2);

#postscript("test.eps",paper="a4") # Creating a ps file 
par(mfrow=c(3,1));
plot(yt);plot(zt1);plot(zt2)
#dev.off();                      # Turing off ps device

# Create lagged ts
yt1<-lag(yt,k=1); yt2<-lag({t,k=2); yt3<-lag(yt,k=3);
windows()
plot(yt,yt3,xy.labels=FALSE,xy.lines=FALSE) # Fig 9.2

# Calculate and plot sac (acf) and spac (pacf) functions
acfyt<-acf(yt); acfzt1<-acf(zt1); acf(zt2);
pacfyt<-pacf(yt);pacf(zt1); pacf(zt2);
par(mfrow=c(2,1));
acfyt<-acf(yt); pacfzt1<-pacf*zt1);
acfyt$acf;  # sacf of Fig 9.3
acfzt1$acf; # sacf of Fig 9.7

acf(zt1); pacf(zt1);  # sac and psac plots to tentatively id a model
# Refer to P422-423. Example 9.6-part1: MA(1) fits zt better contrasting to AR(1)

library(tseries)
z.arma<- arma(zt1,include.intercept = FALSE,order = c(0, 1))
summary(z.arma);plot(z.arma);  # Example 9.6-part2
z.arma$coef; z.arma$fit;

z.arima<-arima0(zt1,include.mean = FALSE,order = c(0,0,1), method = "CSS")
yt.arima<-arima0(yt,include.mean = FALSE,order = c(0,1,1), method = "CSS")
fyt<-yt-yt.arima$res  # Fitted value of yt  (in the past)
fyt; ts.plot(yt,fyt);
py121<-predict(yt.arima,1);
cbind(py121$pred,py121$pred-1.96*py121$se,py121$pred+1.96*py121$se)