###  Assignment #3 R code


setwd("/Users/bob/courses/stat910/assign/")

### load data and SS functions

load("tsa3.rda")



##############################################################
#  3.9, S&S
##############################################################

length(cmort)

time <- 1:508  # weekly smoothed averages during 1970s


# --- sequence plot

plot(time,cmort, type="o")

# --- correlations; partial correlations: signature AR(2)

acf2(cmort, max.lag=60)

# --- fit suggested AR(2) using ols and yule walker
#     similar unless you raise the order considerably

order <- 2;

ar2.ols <- ar(cmort, aic=FALSE, order.max=order, demean=TRUE, intercept=TRUE, method="ols")

ar2.yw  <- ar(cmort, aic=FALSE, order.max=order, demean=TRUE, intercept=TRUE, method="yw") # need demean

summary(ar2.ols) # not very helpful

names(ar2.ols)
ar2.ols; ar2.ols$asy.se.coef; ar2.ols$var.pred

#       1       2  
#  0.4286  0.4418  
#  0.0398  0.0398

# --- Yule Walker is very similar when order =2
ar2.yw; ar2.yw$asy.se.coef; ar2.yw$var.pred

plot(ar2.ols$ar, ar2.yw$ar); abline(a=0,b=1,col="red")  # plot them when order is high

# --- prediction intervals, using AR(2)

pred <- predict(ar2.ols, cmort, n.ahead=20, se.fit=TRUE); names(pred)

# [1] 87.59986 86.76349 87.33714 87.21350
# [1]  5.68485  6.18497  7.13423  7.59335

# --- draw the predictions with interval at end of data

n<-508
show <- 450:n
plot(time[show],cmort[show], xlim=c(show[1],520))
lines(time[show],cmort[show]-ar2$resid[show], col="red")

nf <- length(pred$pred);
lines(n+1:nf, pred$pred, lty=2, col="red")
lines(n+1:nf, pred$pred+1.96*pred$se, col="red")
lines(n+1:nf, pred$pred-1.96*pred$se, col="red")

# --- residual correlations show problems

ar2 <- arima(cmort, order=c(2,0,0), include.mean=TRUE, method="ML"); print(ar2)
#        ar1     ar2  intercept
#       0.4301  0.4424    88.8416
# s.e.  0.0397  0.0398     1.9415

# --- need to go to a long-enough lag to see the annual cycle

tsdiag(ar2, gof.lag=70)


##############################################################
#  3.31, S&S
##############################################################


Temp <- read.table("http://www.stat.pitt.edu/stoffer/tsa2/data/globtemp2.dat")

# --- means are deviations from "normal"

colnames(Temp) <- c("Year", "Mean", "5-Year Mean")

time <- 1880:2004
xt   <- Temp[,2]

# --- begin with sequence plot - always

plot(time, xt, type="o")

# --- evidently not stationary, so difference it

dtime <- time[-1]
dxt <- diff(xt)

plot(dtime, dxt, type="o", xlab="Change in Temperature", ylab="Year")

# --- correlations; partial correlations

acf2(dxt);

# --- fit several low-order arma models

m300 <- arima(dxt, order=c(3,0,0), include.mean=TRUE, method="ML"); print(m300)
tsdiag(m300)

m003<- arima(dxt, order=c(0,0,3), include.mean=TRUE, method="ML"); print(m003)
tsdiag(m003)

m303<- arima(dxt, order=c(3,0,3), include.mean=TRUE, method="ML"); print(m303)
tsdiag(m303)

# --- put the AIC into a matrix  (indices off by 1)
maxp <- maxq <- 6
aic.table <- matrix(NA,nrow=1+maxp, ncol=1+maxq)
rownames(aic.table)<-paste("p",0:maxp,sep="=")     # nice labels help!
colnames(aic.table)<-paste("q",0:maxq,sep="=")

for (p in 0:maxp) for (q in 0:maxq) {
	aic.table[1+p,1+q] <- arima(dxt, order=c(p,0,q), include.mean=TRUE, method="ML")$aic
}

aic.table

pp <- 0:(1+maxp); qq <- 0:(1+maxq)  # one more than nrow or ncol
image(pp,qq,aic.table, col=topo.colors(50))

# better image without q=0 which is really poor fit
image(pp,qq[-1],aic.table[,-1], col=topo.colors(50))

# --- MA(4) seems to work as well as any, at least over these choices of fits 
m004<- arima(dxt, order=c(0,0,4), include.mean=TRUE, method="ML"); print(m004)
tsdiag(m004)

# --- show predictions...  let ARIMA do the predictions.  Oddly, estimates differ a bit

m014 <- arima(xt, order=c(0,1,4), include.mean=TRUE, method="ML"); print(m014); print(m004) 
tsdiag(m014)

n<-length(xt)
show <- 100:n

plot(time[show],xt[show], ylim=c(0,0.85),xlim=c(1979,2010))
lines(time[show],xt[show]-residuals(m014)[show], col="red")

pred <- predict(m014, n.ahead=5)

nf <- length(pred$pred);
lines(2004+1:nf, pred$pred, lty=2, col="red")
lines(2004+1:nf, pred$pred+1.96*pred$se, col="red")
lines(2004+1:nf, pred$pred-1.96*pred$se, col="red")