Code for ‘Data Analysis And Graphics Using R’, 3rd edn, CUP, 2010
John Maindonald and John Braun
Oct 3, 2014
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.5if(any(!z)){
notAvail <- paste(names(z)[!z], collapse=", ")
warning(paste("The following packages should be installed:", notAvail))
}
g9.1()
g9.2A()
g9.2BC()
g9.3()
g9.4()
g9.5()
g9.6()
g9.7()
g9.8()
g9.9()
g9.10()## Series: NA3rtRain
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 intercept SOI CO2
## 0.209 5.196 0.038 0.009
## s.e. 0.100 0.611 0.006 0.002
##
## sigma^2 estimated as 0.198: log likelihood=-68.95
## AIC=147.9 AICc=148.4 BIC=161.5