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))
    on.exit(par(oldpar))
    data(LakeHuron)
    plot(LakeHuron,ylab="Depth (in feet)", xlab = "Time (in years)")
    if(device!="")dev.off()
  }

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

g9.2b <-
function(device=""){
    if(device!="")hardcopy(width=3.25, height=2, device=device)
    oldpar <- par(mar=c(4.1,4.1,3.6,1.1), mfrow=c(1,1),
                  mgp=c(2.25, 0.75, 0))
    on.exit(par(oldpar))
    data(LakeHuron)
    acf(LakeHuron, main="")
    mtext(side=3, line=0.75, "B: Estimated autocorrelations", adj=0)
    if(device!="")dev.off()
  }

g9.3 <-
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))
    data(LakeHuron)
    acf(LakeHuron, type="partial", main="")
    if(device!="")dev.off()
  }

g9.4 <-
function(device="", lag.max=11){
    if(device!="")hardcopy(width=3.75, height=2.75,
                           device=device, pointsize=c(7,7), trellis=T)
    trellis.par.set(superpose.line=list(lwd=c(2,1), lty=c(1,1),
                      col=c("black","gray")))
    library(lattice)
    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)-0.095,
                                            rep(1:lag.max,4)+0.095),
                     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),
                 scale = list(cex = 0.75), groups=gp,
                 par.strip.text = list(cex = 0.65),
                 type = "h", ylab = "Autocorrelation"))
    if(device!="") dev.off()
  }

g9.5 <-
function(device=""){
    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))
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
    xbomsoi <- with(bomsoi, data.frame(SOI=SOI, cuberootRain=avrain^0.33))
    xbomsoi$trendSOI <- lowess(xbomsoi$SOI)$y
    xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))
            xbomsoi$detrendRain <-
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))    
            xbomsoi$detrendSOI <-
            with(xbomsoi, SOI - trendSOI + mean(trendSOI))
            xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y    
  }
    ## Plot time series avrain and SOI: ts object bomsoi (DAAG)
    plot(ts(xbomsoi[, c("cuberootRain","SOI")], start=1900), 
         panel=function(y,...)panel.smooth(1900:2004, y,...),
         xlab = "Year", main="", ylim=list(c(250, 800),c(-20,25)))
    if(device!="") dev.off()
    invisible(xbomsoi)
  }

g9.6 <-
function(device=""){
    if(device!="")hardcopy(width=2.4, height=2.4, device=device)
    oldpar <-  par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.75,0.75,0))
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
    xbomsoi <- with(bomsoi, data.frame(SOI=SOI, cuberootRain=avrain^0.33))
    xbomsoi$trendSOI <- lowess(xbomsoi$SOI)$y
    xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))
            xbomsoi$detrendRain <-
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))    
            xbomsoi$detrendSOI <-
            with(xbomsoi, SOI - trendSOI + mean(trendSOI))
            xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y    
  }
    rainpos <- pretty(bomsoi$avrain, 5)
    with(xbomsoi,
         {plot(cuberootRain ~ SOI, xlab = "SOI",
               ylab = "Rainfall (cube root scale)", yaxt="n")
          axis(2, at = rainpos^0.33, labels=paste(rainpos))
          ## Relative changes in the two trend curves
          lines(lowess(cuberootRain ~ SOI))
          lines(lowess(trendRain ~ trendSOI), lwd=1.75, col="gray40")
        })    
    if(device!="") dev.off()
    invisible(xbomsoi)
  }

g9.7 <-
function(device=""){
    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))
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
            xbomsoi <-
            with(bomsoi, data.frame(SOI=SOI, cuberootRain=avrain^0.33))
            xbomsoi$trendSOI <- lowess(xbomsoi$SOI)$y
            xbomsoi$detrendRain <-
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))
            xbomsoi$detrendSOI <-
            with(xbomsoi, SOI - trendSOI + mean(trendSOI))
            xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y
          }
    rainpos <- pretty(xbomsoi$cuberootRain^3, 6)
    plot(cuberootRain ~ SOI, data = xbomsoi, 
         ylab = "Rainfall (cube root scale)", yaxt="n")
    axis(2, at = rainpos^0.33, labels=paste(rainpos))
    mtext(side = 3, line = 0.25, "A", adj = 0)
    with(xbomsoi, lines(lowess(cuberootRain ~ SOI)))
    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)))
    mtext(side = 3, line = 0.25, "B", adj = 0)
    if(device!="") dev.off()
  }

g9.8 <-
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), mfrow=c(1,1), mgp=c(2.25, 0.75, 0))
    on.exit(par(oldpar))
    if(!exists("xbomsoi")){
    xbomsoi <- with(bomsoi, data.frame(SOI=SOI, cuberootRain=avrain^0.33))
    xbomsoi$trendSOI <- lowess(xbomsoi$SOI)$y
    xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))
            xbomsoi$detrendRain <-
            with(xbomsoi, cuberootRain - trendRain + mean(trendRain))    
            xbomsoi$detrendSOI <-
            with(xbomsoi, SOI - trendSOI + mean(trendSOI))
            xbomsoi$trendRain <- lowess(xbomsoi$cuberootRain)$y    
  }    
    acf(resid(lm(detrendRain ~ detrendSOI, data = xbomsoi)), main="")
    if(device!="")dev.off()
  }

g9.9 <-
function(device="", lag.max=11){
    if(device!="")hardcopy(width=3.75, height=2.75,
                           device=device, pointsize=c(7,7), trellis=T)
    trellis.par.set(superpose.line=list(lwd=c(2,1), lty=c(1,1),
                      col=c("black","gray")))
    on.exit(par(oldpar))
    library(lattice)
    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)-0.095,
                                            rep(1:lag.max,4)+0.095),
                     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),
                 scale = list(cex = 0.75), groups=gp,
                 par.strip.text = list(cex = 0.65),
                 type = "h", ylab = "Autocorrelation"))
    if(device!="") dev.off()
  }

gdump <-
function(fnam=NULL, prefix="~/r-book/ed2/figures/figs",
           splitchar="/ch"){
    if(is.null(fnam)){
      path <- getwd()
      pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]]
      fnam <- paste(prefix, pathtag[length(pathtag)],
                    ".R", sep="")
    }
    else fnam <- paste(prefix, fnam, sep="/")
    objnames <- c(objects(pattern="^g", envir=sys.frame(0)), "hardcopy")
    cat("\nDump to file:", fnam, "\n")
    print(objnames)
    dump(objnames, fnam)
  }

gsave <-
function(fnam=NULL, prefix="~/r-book/ed2/figures/figs",
           splitchar="/ch", xtras=c("hardcopy", "renum.fun","renum.files")){
    if(is.null(fnam)){
      path <- getwd()
      pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]]
      fnam <- paste(prefix, pathtag[length(pathtag)],
                    ".RData", sep="")
    }
    else fnam <- paste(prefix, fnam, sep="/")
    objnames <- c(objects(pattern="^g", envir=sys.frame(0)), xtras)
    cat("\nDump to file:", fnam, "\n")
    print(objnames)
    save(list=objnames, file=fnam)
  }

hardcopy <-
function(width=3.75, height=3.75, color=FALSE, trellis=FALSE,
           device=c("","pdf","ps"), path="~/r-book/ed2/Art/", file=NULL,
           format=c("nn-nn", "name"), split="\\.", pointsize=c(8,4),
           fonts=NULL,
           horiz=FALSE, ...){
    if(!trellis)pointsize <- pointsize[1]
    funtxt <- sys.call(1)
    nam <- strsplit(as.character(funtxt), "(", fixed=TRUE)[[1]][1]
    suffix <- switch(device, ps=".eps", pdf=".pdf")
    if(is.character(path) & nchar(path)>1 & substring(path, nchar(path))!="/")
      path <- paste(path, "/", sep="")
    if(is.null(file)) if(format[1]=="nn-nn"){
      if(!is.null(split))dotsplit <- strsplit(nam, split)[[1]] else
      dotsplit <- nam
      if(length(dotsplit)==1)dotsplit <- c("", dotsplit)
      nn2 <- paste(if(nchar(dotsplit[2])==1)"0" else "", dotsplit[2],
                   sep="")
      if(nchar(dotsplit[1])>0){
        numstart <- which(unlist(strsplit(dotsplit[1], "")) %in% paste(0:9))[1]
        nn1 <- substring(dotsplit[1], numstart)
        nn1 <- paste(if(nchar(nn1) == 1) "0" else "", nn1, "-", sep="")
      } else nn1 <- ""
      file <- paste(nn1, nn2, sep="")
    } else file <- nam
    if(nchar(file)>4 & substring(file, nchar(file)-nchar(suffix)+1)==suffix)
      suffix <- ""
    file <- paste(path, file, suffix, sep="")
    print(paste("Output will be directed to file:", file))
    dev.out <- device[1]
    dev.fun <- switch(dev.out, pdf=pdf, ps=postscript)
    if(trellis){
      library(lattice)
      if(device=="ps")
        trellis.device(file=file, device=dev.fun,
                       color = color, horiz=horiz, fonts=fonts,
                       width=width, height=height, ...) else
      trellis.device(file=file, device=dev.fun, fonts=fonts,
                     color = color, width=width, height=height, ...)

      trellis.par.set(list(fontsize=list(text=pointsize[1], points=pointsize[2])))
    } else
    if (dev.out!=""){
      print(c(width, height))
      if(device=="ps")
        dev.fun(file=file, paper="special",  horiz=horiz, fonts=fonts,
                width=width, height=height, pointsize=pointsize[1], ...) else
      dev.fun(file=file, paper="special", fonts=fonts,
              width=width, height=height, pointsize=pointsize[1], ...)      
    }
    if(trellis)trellis.par.set(list(fontsize=list(text=pointsize[1],
                                      points=pointsize[2])))
  }