Chapter 9: Time Series Models
Packages required: “DAAG”, “lattice”, “MASS”, “forecast”, “grid”

The script that follows is designed to be executed as it stands. It sets up
functions that create Ch 9 graphs, then runs those functions. If executed by
clicking on the RStudio “Compile Notebook …” command, it will be processed
through R Markdown, generating a document that includes code, output and
graphs generated by executing the functions. If some needed packages are
missing, warning messages will appear. See further:
http://rmarkdown.rstudio.com/r_notebook_format.html

g9.1 <-
function(device=""){
    if(device!="")hardcopy(width=4, height=2, device=device)
    oldpar <- par(mar=par()$mar-c(.5,0,1.5,0), mfrow=c(1,1), tcl=-0.35)
    on.exit(par(oldpar))
    plot(LakeHuron,ylab="Depth (in feet)", xlab = "Time (in years)")
    if(device!="")dev.off()
  }


g9.2A <-
function(device=""){
    if(device!="")hardcopy(width=5.0, height=1.65, device=device)
    oldpar <- par(pty="s")
    on.exit(par(oldpar))
    lag.plot(LakeHuron, set.lags=1:4,do.lines=F, oma=c(0,1.5,3.1,1.5),
             layout=c(1,4), cex.lab=1.15)
    mtext(side=3, line=3, "A: Lag plots", adj=-0.06, cex=1.1)
    if(device!="")dev.off()
  }

g9.2BC <-
function(device="", lag.max=15, offset=0.125, col2="gray"){
    if(device!="")hardcopy(width=5.25, height=2.15, device=device,
                           pointsize=c(text=8,points=3),)
    if(!require(lattice))return("Package 'lattice' must be installed.")       
    ci95 <- 2/sqrt(length(LakeHuron))
    gph.key <-  list(space="top", text=list(c("Autocorrelation",
                                    "Partial autocorrelation")),
                     columns=2,
                     lines=list(lwd=c(1.25,2), col=c("black",col2)),
                     padding.text=1)
    gph0 <- xyplot(1 ~ 1, panel=function(x,y,...){}, xlab="", ylab="",
                   scales=list(draw=F), key=gph.key, aspect=1,
                   par.settings=list(axis.line=list(col="transparent"),
                     fontsize=list(text=8, points=3)))

    acfval <- acf(LakeHuron, main="", plot=FALSE, lag.max=lag.max)$acf
    pacfval <- pacf(LakeHuron, main="", plot=FALSE, lag.max=lag.max)$acf
    xy <- data.frame(acf=c(acfval,pacfval),
                     Lag=c(0,rep(1:lag.max,2))+c(0,rep(c(-offset,offset),
                       rep(lag.max,2))),
                     gp=rep(c("acf","pacf"), c(lag.max+1,lag.max)))
    if(!require(grid))return("Package 'grid' must be installed.")      
    obsgph <- xyplot(acf ~ Lag, data = xy, groups=gp, type=c("h"),
                     par.strip.text = list(cex = 0.85), 
                     ylab = "Autocorrelation", origin=0, ylim=c(-0.325, 1.04),
                     par.settings=simpleTheme(
                       col=c("black",col2), lty=rep(1,2),
                       lwd=c(1.25,2.5), pch=20),
                     scales=list(alternating=FALSE, tck=0.5),
                     main=textGrob("B: LakeHuron series",
                       x=unit(.075, "npc"), y = unit(0, "npc"), just="left", 
                       gp=gpar(fontsize=9)),           
                     panel=function(x,y,...){
                       panel.xyplot(x,y,...)
                       panel.abline(h=0, lwd=0.8)
                       panel.abline(h=ci95, lwd=0.8, lty=2)
                       panel.abline(h=-ci95, lwd=0.8, lty=2)                       
                     }
                     )
    print(obsgph, pos=c(0,0,0.5,0.9))
    acfval <- ARMAacf(ar=0.8, lag.max=15)
    pacfval <- ARMAacf(ar=0.8, lag.max=15, pacf=TRUE)
    xy <- data.frame(acf=c(acfval,pacfval),
                     Lag=c(0,rep(1:lag.max,2))+c(0,rep(c(-offset,offset),
                       rep(lag.max,2))),
                     gp=rep(c("acf","pacf"), c(lag.max+1,lag.max)))
    obsgph <- xyplot(acf ~ Lag, data = xy, groups=gp, type=c("h"),
                     par.strip.text = list(cex = 0.85), 
                     ylab = "Autocorrelation", origin=0, ylim=c(-0.325, 1.04),
                     par.settings=simpleTheme(
                       col=c("black",col2), lty=rep(1,2),
                       lwd=c(1.25,2.5), pch=20), 
                     scales=list(alternating=FALSE, tck=0.5),
                     main=textGrob("C: Theoretical AR1 process",
                       x=unit(.075, "npc"), y = unit(0, "npc"), just="left", 
                       gp=gpar(fontsize=9)),                                  
                     panel=function(x,y,...){
                       panel.xyplot(x,y,...)
                       panel.abline(h=0, lwd=0.8)
                     }
                     )
    print(obsgph, pos=c(0.5,0,1,0.9),newpage=FALSE)
    print(gph0, pos=c(0,0,1,1),newpage=FALSE)
    if(device!="")dev.off()
  }

g9.3 <-
function(device="", lag.max=11, col2="gray", offset=0.09){
    if(device!="")hardcopy(width=4.5, height=3,
                           device=device, pointsize=c(7,4), trellis=TRUE)
    trellis.par.set(superpose.line=list(lwd=c(1.25,2.25), lty=c(1,1),
                      col=c("black",col2)))
    if(!require(lattice))return("Package 'lattice' must be installed.")
    four <- c(ARMAacf(ma=c(0.5),lag.max=lag.max),
              ARMAacf(ma=c(0,0,0.5), lag.max=lag.max),
              ARMAacf(ma=c(0,0,0.5,0,0.5),lag.max=lag.max),
              ARMAacf(ma=rep(0.5,6),lag.max=lag.max))
    pfour <- c(ARMAacf(ma=c(0.5),lag.max=lag.max, pacf=T),
               ARMAacf(ma=c(0,0,0.5), lag.max=lag.max, pacf=T),
               ARMAacf(ma=c(0,0,0.5,0,0.5),lag.max=lag.max, pacf=T),
               ARMAacf(ma=rep(0.5,6),lag.max=lag.max, pacf=T))
    xy <- data.frame(acf = c(four,pfour), Lag = c(rep(0:lag.max, 4)-offset,
                                            rep(1:lag.max,4)+offset),
                     gp=rep(c("acf", "pacf"), c((lag.max+1)*4, lag.max*4)),
                     series =
                     c(rep(c("ma = 0.5",
                             "ma = c(0, 0, 0.5)",
                             "ma = c(0, 0, 0.5, 0, 0.5)",
                             "ma = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5)"),
                           rep(lag.max+1,4)),
                       rep(c("ma = 0.5",
                             "ma = c(0, 0, 0.5)",
                             "ma = c(0, 0, 0.5, 0, 0.5)",
                             "ma = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5)"),
                           rep(lag.max,4))))
    
    print(xyplot(acf ~ Lag|series, data = xy, layout = c(2,2), groups=gp,
                 par.strip.text = list(cex = 0.85), type=c("h"),
                 ylab = "Autocorrelation", origin=0,
                 as.table=TRUE, between=list(x=0.5, y=0.5),
                 par.settings=simpleTheme(
                   col=c("black",col2),
                   lwd=c(1.25,2), pch=20),
                 key=list(space="top", columns=2, text=list(c("Autocorrelation",
                                                    "Partial autocorrelation")),
                   lines=list(lwd=c(1.25,2), col=c("black",col2)),
                   points=list(pch=c(20,20), col=c("black",col2))),
                 scales=list(alternating=FALSE, tck=0.5),
                 panel=function(x,y,...){
                   panel.xyplot(x,y,...)
                   panel.abline(h=0, lwd=0.8)
                 }
                 ))
    if(device!="") dev.off()
  }

g9.4 <-
function( device=""){
    if(device!="")hardcopy(width=3.25, height=1.8,
                           device=device)
    oldpar <- par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.25, 0.75, 0), mfrow=c(1,1))
    on.exit(par(oldpar))
    if(!require(forecast))return("Package 'forecast' must be installed.")      
    LHfit <- auto.arima(LakeHuron)
    fcast <- forecast(LHfit)
    plot(fcast, shadecols=c("gray","gray60"), tcl=-0.35)
    if(device!="")dev.off()
  }

g9.5 <-
function(device="", lag.max=8, ma3=c(.125, .25, .375, .5), offset=0.085,
           n=98){
    if(device!="")hardcopy(width=4.5, height=3,
                           device=device, pointsize=c(text=7,points=4),
                           trellis=T)
    simfun <- function(ma, n, lag.max){
      z1 <- acf(arima.sim(model=list(ma=ma), n=n), lag.max=lag.max,
                plot=FALSE)$acf
      z4 <- ARMAacf(ma=ma,lag.max=lag.max)
      z6 <- acf(arima.sim(model=list(ma=ma), n=n), lag.max=lag.max,
                plot=FALSE)$acf
      c(z1[-1],z4,z6[-1])
    }
    if(!require(lattice))return("Package 'lattice' must be installed.")      
    five1 <- simfun(ma=c(0,0,ma3[1]), n=n, lag.max=lag.max)
    five2 <- simfun(ma=c(0,0,ma3[2]), n=n, lag.max=lag.max)
    five3 <- simfun(ma=c(0,0,ma3[3]), n=n, lag.max=lag.max)
    five4 <- simfun(ma=c(0,0,ma3[4]), n=n, lag.max=lag.max)
    xy <- data.frame(acf = c(five1,five2,five3,five4),
                     Lag=rep(c(1:lag.max, 0:lag.max,
                       1:lag.max),4),
                     ma3 = rep(paste("ma = c(0, 0, ", ma3,")",sep=""),
                       rep(lag.max*3+1,4)),
                     gp=rep(rep(1:3, c(lag.max, lag.max+1,
                       lag.max)), 4))
    xy$Lag <- xy$Lag+(xy$gp-1)*offset
    print(xyplot(acf~Lag|ma3, groups=gp, data = xy, layout = c(2,2),
                 type=c("h"),
                 par.strip.text = list(cex = 0.85),
                 ylab = "Autocorrelation", origin=0, as.table=TRUE,
                 between=list(x=0.5, y=0.5),
                 par.settings=simpleTheme(
                   col=rep(c("gray","black","gray"),c(1,1,1)),
                   lty=rep(1,3),
                   lwd=c(2,1.25,2), pch=20),
                 key=list(space="top", columns=2, text=list(c("Theoretical",
                                                    "Simulation")),
                   lines=list(lwd=c(1.25,2), col=c("black","gray"))),
                 scales=list(x=list(at=0:8), alternating=FALSE, tck=0.5),
                 panel=function(x,y,...){
                   panel.xyplot(x,y,...)
                   panel.abline(h=0, lwd=0.75)
                 }
                 ))
    if(device!="") dev.off()
  }

g9.6 <-
function(device="", df=bomregions2012, varnam="mdRain"){
    if(device!="")hardcopy(width=3.25, height=3.25, device=device)
    oldpar <-  par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.75, 0.75,0))
    if(!require(DAAG))return("Package 'DAAG' must be installed.")     
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
      avrain <- df[,"mdbRain"]
      xbomsoi <- with(df, data.frame(SOI=SOI, mdb3rtRain=avrain^0.33))
      xbomsoi$trendSOI <- lowess(xbomsoi$SOI, f=0.1)$y
      xbomsoi$trendRain <- lowess(xbomsoi$mdb3rtRain, f=0.1)$y
      xbomsoi$detrendRain <-
        with(xbomsoi, mdb3rtRain - trendRain + mean(trendRain))    
      xbomsoi$detrendSOI <-
        with(xbomsoi, SOI - trendSOI + mean(trendSOI))
    }
    ## Plot time series avrain and SOI: ts object xbomsoi
    tsmdb <- ts(xbomsoi[, c("mdb3rtRain","SOI")], start=1900)
    plot(tsmdb, 
         panel=function(y,...)panel.smooth(time(tsmdb), y,...),
         xlab = "Year", main="", ylim=list(c(250, 800),c(-20,25)), tcl=-0.35)
    if(device!="") dev.off()
    invisible(xbomsoi)
  }

g9.7 <-
function(device="", df=bomregions2012){
    if(device!="")hardcopy(width=4.25, height=2, device=device)
    oldpar <-  par(mar=c(4.1,4.6,2.1,1.1), mfrow = c(1,2), mgp=c(2.5,0.75,0),
                   oma = c(0, 0, 0, 0), pty="s")
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
      avrain <- df[,"mdbRain"]
      xbomsoi <- with(df, data.frame(SOI=SOI, mdb3rtRain=avrain^0.33))
      xbomsoi$trendSOI <- lowess(xbomsoi$SOI, f=0.1)$y
      xbomsoi$trendRain <- lowess(xbomsoi$mdb3rtRain, f=0.1)$y
      xbomsoi$detrendRain <-
        with(xbomsoi, mdb3rtRain - trendRain + mean(trendRain))    
      xbomsoi$detrendSOI <-
        with(xbomsoi, SOI - trendSOI + mean(trendSOI))
    }
    rainpos <- pretty(xbomsoi$mdb3rtRain^3, 6)
    plot(mdb3rtRain ~ SOI, data = xbomsoi, 
         ylab = "Rainfall (cube rt scale)", yaxt="n")
    axis(2, at = rainpos^0.33, labels=paste(rainpos))
    mtext(side = 3, line = 0.8, "A", adj = -0.025)
    with(xbomsoi, lines(lowess(mdb3rtRain ~ SOI, f=0.75)))
    plot(detrendRain ~ detrendSOI, data = xbomsoi,
         xlab="Detrended SOI", ylab = "Detrended rainfall", yaxt="n")
    axis(2, at = rainpos^0.33, labels=paste(rainpos)) 
    with(xbomsoi, lines(lowess(detrendRain ~ detrendSOI, f=0.75)))
    mtext(side = 3, line = 0.8, "B", adj = -0.025)
    if(device!="") dev.off()
  }

g9.8 <-
function(device="", df=bomregions2012){
    if(device!="")hardcopy(width=5.5, height=2, device=device)
    oldpar <- par(mar=c(4.1,4.1,2.6,1.1), mfrow=c(1,1), mgp=c(2.25, 0.75, 0),
                  pty="s")
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){  
      xbomsoi <- with(df, data.frame(SOI=SOI, mdb3rtRain=mdbRain^0.33))
      xbomsoi$trendSOI <- lowess(xbomsoi$SOI, f=0.1)$y
      xbomsoi$trendRain <- lowess(xbomsoi$mdb3rtRain, f=0.1)$y
      xbomsoi$detrendRain <-
        with(xbomsoi, mdb3rtRain - trendRain + mean(trendRain))    
      xbomsoi$detrendSOI <-
        with(xbomsoi, SOI - trendSOI + mean(trendSOI))
    }
    par(fig=c(0,0.5,0,1))
    acf(resid(lm(mdb3rtRain ~ SOI, data = xbomsoi)), main="")
    mtext(side = 3, line = 1, "A", adj = -0.025)
    par(fig=c(0.5,1,0,1), new=TRUE)
    pacf(resid(lm(mdb3rtRain ~ SOI, data = xbomsoi)), main="", ylim=c(-0.2,1))
    mtext(side = 3, line = 1, "B", adj = -0.025)
    if(device!="")dev.off()
  }

g9.9 <-
function(device="", df=bomregions2012){
    if(device!="")hardcopy(width=5.5, height=2, device=device)
    oldpar <- par(mar=c(4.1,4.1,2.1,1.1), mfrow=c(1,1), mgp=c(2.25, 0.75, 0),
                  pty="s")
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
      avrain <- df[,"mdbRain"]
      xbomsoi <- with(df, data.frame(SOI=SOI, mdb3rtRain=avrain^0.33))
      xbomsoi$trendSOI <- lowess(xbomsoi$SOI, f=0.1)$y
      xbomsoi$trendRain <- lowess(xbomsoi$mdb3rtRain, f=0.1)$y
      xbomsoi$detrendRain <-
        with(xbomsoi, mdb3rtRain - trendRain + mean(trendRain))    
      xbomsoi$detrendSOI <-
        with(xbomsoi, SOI - trendSOI + mean(trendSOI))
      xbomsoi$trendRain <- lowess(xbomsoi$mdb3rtRain)$y    
    }
    par(fig=c(0,0.5,0,1))
    md.arima <- with(xbomsoi, arima(mdb3rtRain,  xreg=SOI,
                                    order=c(0,1,1)))
    res <- as.vector(resid(md.arima))
    plot(res ~ SOI, ylab="Residuals", type="p", data=xbomsoi)
    mtext(side = 3, line = 1, "A: Residuals", adj = -0.025)
    par(fig=c(0.5,1,0,1), new=TRUE)
    acf(resid(md.arima), main="")
    mtext(side = 3, line = 1, "B: Autocorrelation plot; residuals",
          adj = -0.025)
    if(device!="")dev.off()
  }

g9.10 <-
function(device="", df=bomregions2012, robust=TRUE){
    if(!require(MASS))return("Package 'MASS' must be installed.")  
    if(!require(forecast))return("Package 'forecast' must be installed.")      
    if(device!="")hardcopy(width=5.5, height=4, device=device)
    oldpar <- par(mar=c(3.6,4.1,2.1,1.1), mfrow=c(1,1), mgp=c(2.25, 0.75, 0),
                  cex.lab=0.9, cex.axis=0.9)
    on.exit(par(oldpar))
    ## A: N Aus rainfall vs year
    par(fig=c(0,0.5,0.5,1))
    NA3rtRain <- df$northRain^(1/3)
    with(df, plot(NA3rtRain ~ Year, yaxt="n",
                  xlab="", ylab="northRain (mm)"))
    with(df, lines(lowess(NA3rtRain ~ Year, f=0.75)))
    yloc <- pretty(df$northRain)
    axis(2, at=yloc^(1/3), labels=paste(yloc))
    mtext(side = 3, line = 0.5, "A: N. Aust. rainfall vs Year",
          adj = -0.025)    
    ## B: NA3rtRain vs SOI
    par(fig=c(0.5,1, 0.5, 1), new=TRUE)
    plot(NA3rtRain ~ SOI, data=df)
    lines(lowess(NA3rtRain ~ df$SOI, f=0.75))
    mtext(side = 3, line = 0.5,
          "B: NA3rtRain vs SOI",
          adj = -0.025)
    ## C: Partial residual plot; include CO2?
    par(fig=c(0,0.5,0,0.5), new=TRUE)
    NA3rtRainSOI.res <- resid(lqs(NA3rtRain ~ SOI, data=df))
    CO2SOI.res <- resid(lqs(CO2 ~ SOI, data=df))
    plot(NA3rtRainSOI.res ~ CO2SOI.res)
    lines(lowess(NA3rtRainSOI.res ~ CO2SOI.res, f=0.75))
    mtext(side = 3, line = 0.5,
          "C: Partial residual plot; include CO2?",
          adj = -0.025)
    ## D: ARIMA(0,0,1), xreg=cbind(SOI, CO2))
    par(fig=c(0.5,1, 0,0.5), new=TRUE)
    north.arima <- with(df, auto.arima(NA3rtRain,
                                       xreg=cbind(SOI, CO2)))
    print(north.arima)
    plot(unclass(resid(north.arima)) ~ Year, type="p", data= df,
         xlab="", ylab="Resids, ARIMA(0,0,1)",
         cex=0.8, cex.axis=0.9, cex.lab=0.9)
    mtext(side = 3, line = 0.5, adj = -0.025,
          "D: ARIMA(0,0,1), xreg=cbind(SOI, CO2))")    
    if(device!="")dev.off()
  }

pkgs <- c("DAAG", "lattice", "MASS", "forecast", "grid")
z <- sapply(pkgs, require, character.only=TRUE, warn.conflicts=FALSE, quietly=TRUE)
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## This is forecast 5.5
if(any(!z)){
  notAvail <- paste(names(z)[!z], collapse=", ")
  warning(paste("The following packages should be installed:", notAvail))
}

g9.1()

plot of chunk unnamed-chunk-1

g9.2A()

plot of chunk unnamed-chunk-1

g9.2BC()

plot of chunk unnamed-chunk-1

g9.3()

plot of chunk unnamed-chunk-1

g9.4()