## Chapter 9: Time Series Models
## Sec 9.1: Time Series -- Some Basic Ideas
## ss 9.1.1: Preliminary graphical explorations
## Plot depth measurements: ts object LakeHuron (datasets)
plot(LakeHuron, ylab="depth (in feet)", xlab = "Time (in years)")
opar <- par(mfrow=c(1,3), pty="s")
lag.plot(LakeHuron, lags=3, do.lines=FALSE)
par(opar)
## ss 9.1.2: The autocorrelation and partial autocorrelation function
acf(LakeHuron)
## pacf(LakeHuron) gives the plot of partial autocorelations
## ss 9.1.3: Autoregressive (AR) models
## The { AR(1)} model
## Footnote Code
## Yule-Walker autocorrelation estimate
LH.yw <- ar(x = LakeHuron, order.max = 1, method = "yw")
# autocorrelation estimate
# order.max = 1 for the AR(1) model
LH.yw$ar # autocorrelation estimate of alpha
## [1] 0.8319
## Maximum likelihood estimate
LH.mle <- ar(x = LakeHuron, order.max = 1, method = "mle")
LH.mle$ar # maximum likelihood estimate of alpha
## [1] 0.8375
LH.mle$x.mean # estimated series mean
## [1] 579.1
LH.mle$var.pred # estimated innovation variance
## [1] 0.5093
## The general { AR(}p{ )} model
ar(LakeHuron, method="mle")
##
## Call:
## ar(x = LakeHuron, method = "mle")
##
## Coefficients:
## 1 2
## 1.044 -0.250
##
## Order selected 2 sigma^2 estimated as 0.479
## ss 9.1.4: *~Autoregressive moving average (ARMA) models -- theory
## Footnote Code
## series.diff <- diff(series)
## ss 9.1.5: Automatic model selection?
library(forecast)
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: timeDate
## This is forecast 5.5
auto.arima(LakeHuron)
## Series: LakeHuron
## ARIMA(0,1,4) with drift
##
## Coefficients:
## ma1 ma2 ma3 ma4 drift
## 0.058 -0.316 -0.303 -0.235 -0.021
## s.e. 0.103 0.106 0.110 0.130 0.017
##
## sigma^2 estimated as 0.481: log likelihood=-101.1
## AIC=214.2 AICc=215.1 BIC=229.6
## Check that model removes most of the correlation structure
acf(resid(arima(LakeHuron, order=c(p=1, d=1, q=2))))
## The following achieves the same effect, for these data
acf(resid(auto.arima(LakeHuron)))
## ss 9.1.6: A time series forecast
LH.arima <- auto.arima(LakeHuron)
fcast <- forecast(LH.arima)
plot(fcast)
## Use of simulation as a check
oldpar <- par(mfrow=c(2,2), mar=c(3.1,4.6,2.6, 1.1))
for(i in 1:4){
ma3 <- 0.125*i
simts <- arima.sim(model=list(order=c(0,0,3), ma=c(0,0,ma3)), n=98)
acf(simts, main="", xlab="")
mtext(side=3, line=0.5, paste("ma3 =", ma3), adj=0)
}
par(oldpar)
set.seed(29) # Ensure that results are reproducible
estMAord <- matrix(0, nrow=4, ncol=20)
for(i in 1:4){
ma3 <- 0.125*i
for(j in 1:20){
simts <- arima.sim(n=98, model=list(ma=c(0,0,ma3)))
estMAord[i,j] <- auto.arima(simts, start.q=3)$arma[2]
}
}
detectedAs <- table(row(estMAord), estMAord)
dimnames(detectedAs) <- list(ma3=paste(0.125*(1:4)),
Order=paste(0:(dim(detectedAs)[2]-1)))
print(detectedAs)
## Order
## ma3 0 1 2 3 4
## 0.125 15 2 3 0 0
## 0.25 9 4 1 6 0
## 0.375 4 1 2 11 2
## 0.5 1 0 1 18 0
## Sec 9.2: *~Regression modeling with ARIMA errors
library(DAAG)
## Loading required package: lattice
## Plot time series mdbRain and SOI: ts object bomregions (DAAG)
plot(ts(bomregions[, c("mdbRain","SOI")], start=1900),
panel=function(y,...)panel.smooth(bomregions$Year, y,...))
xbomsoi <-
with(bomregions, data.frame(SOI=SOI, mdbRain=mdbRain,
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))
## Footnote Code
## Plot cube root of Murray-Darling basin rainfall vs SOI
oldpar <- par(mfrow=c(1,2), pty="s")
plot(mdb3rtRain ~ SOI, data = xbomsoi, ylab = "Rainfall (cube root scale)", yaxt="n")
rainpos <- pretty(xbomsoi$mdbRain, 6)
axis(2, at = rainpos^0.33, labels=paste(rainpos))
with(xbomsoi, lines(lowess(mdb3rtRain ~ SOI, f=0.1)))
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.1)))
par(oldpar)
par(mfrow=c(1,2), pty="s")
acf(resid(lm(mdb3rtRain ~ SOI, data = xbomsoi)), main="")
pacf(resid(lm(mdb3rtRain ~ SOI, data = xbomsoi)), main="")
par(mfrow=c(1,1), pty="m")
## Use auto.arima() to fit model
(mdb.arima <- with(xbomsoi, auto.arima(mdb3rtRain, xreg=SOI)))
## Series: mdb3rtRain
## ARIMA(0,1,2)
##
## Coefficients:
## ma1 ma2 SOI
## -0.909 -0.027 0.042
## s.e. 0.103 0.104 0.007
##
## sigma^2 estimated as 0.276: log likelihood=-83.55
## AIC=175.1 AICc=175.5 BIC=185.8
## Fit undifferenced model
(mdb0.arima <- with(xbomsoi,
auto.arima(mdb3rtRain, xreg=SOI, d=0)))
## Series: mdb3rtRain
## ARIMA(0,0,0) with non-zero mean
##
## Coefficients:
## intercept SOI
## 7.603 0.040
## s.e. 0.050 0.007
##
## sigma^2 estimated as 0.275: log likelihood=-84.26
## AIC=174.5 AICc=174.7 BIC=182.6
## Footnote Code
## Increase in mdbRain, given a unit increase in SOI
with(xbomsoi, (range(mdbRain)^(1/3) + 0.0427)^3 - range(mdbRain))
## [1] 5.197 11.283
## *The {northRain} series
par(mfrow=c(2,2))
## Panel A
NA3rtRain <- bomregions$northRain^(1/3)
with(bomregions, plot(NA3rtRain ~ Year, yaxt="n",
xlab="", ylab="northRain (mm)"))
with(bomregions, lines(lowess(NA3rtRain ~ Year, f=0.75)))
yloc <- pretty(bomregions$northRain)
axis(2, at=yloc^(1/3), labels=paste(yloc))
## Panel B
plot(NA3rtRain ~ SOI, data=bomregions)
lines(lowess(NA3rtRain ~ bomregions$SOI, f=0.75))
## Panel C
NA3rtRainSOI.res <- resid(lm(NA3rtRain ~ SOI, data=bomregions))
CO2SOI.res <- resid(lm(CO2 ~ SOI, data=bomregions))
plot(NA3rtRainSOI.res ~ CO2SOI.res)
lines(lowess(NA3rtRainSOI.res ~ CO2SOI.res, f=0.75))
(north.arima <- with(bomregions,
auto.arima(NA3rtRain, xreg=cbind(SOI, CO2))))
## Series: NA3rtRain
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 intercept SOI CO2
## 0.204 5.400 0.035 0.008
## s.e. 0.102 0.694 0.007 0.002
##
## sigma^2 estimated as 0.206: log likelihood=-68.46
## AIC=146.9 AICc=147.5 BIC=160.4
## Now plot the residuals, as in Panel D
plot(unclass(resid(north.arima)) ~ Year, type="p",
data= bomregions)
## Fit a constant mean model with the same correlation structure
with(bomregions, arima(x = NA3rtRain, order = c(0, 0, 1), xreg=CO2))
## Series: NA3rtRain
## ARIMA(0,0,1) with non-zero mean
##
## Coefficients:
## ma1 intercept CO2
## 0.223 5.753 0.007
## s.e. 0.108 0.786 0.002
##
## sigma^2 estimated as 0.258: log likelihood=-80.84
## AIC=169.7 AICc=170.1 BIC=180.4
(.258-.206)/.258
## [1] 0.2016
## Other checks and computations
north.arima <- with(bomregions,
auto.arima(NA3rtRain, xreg=cbind(SOI, CO2)))
Box.test(resid(north.arima), lag=20, type="Ljung")
##
## Box-Ljung test
##
## data: resid(north.arima)
## X-squared = 31.52, df = 20, p-value = 0.04866
## Examine normality of estimates of "residuals" ("innovations")
par(pty="s")
qqnorm(resid(north.arima, type="normalized"))
par(pty="m")
(mdb2.arima <- with(xbomsoi, auto.arima(mdb3rtRain,
xreg=poly(SOI,2))))
## Series: mdb3rtRain
## ARIMA(0,1,1)
##
## Coefficients:
## ma1 1 2
## -0.938 2.937 0.851
## s.e. 0.040 0.512 0.509
##
## sigma^2 estimated as 0.269: log likelihood=-82.21
## AIC=172.4 AICc=172.8 BIC=183.1
## Sec 9.3: *Nonlinear time series
x <- numeric(999) # x will contain the ARCH(1) errors
x0 <- rnorm(1)
for (i in 1:999){
x0 <- rnorm(1, sd=sqrt(.25 + .95*x0^2))
x[i] <- x0
}
library(tseries)
garch(x, order = c(0, 1), trace=FALSE)
##
## Call:
## garch(x = x, order = c(0, 1), trace = FALSE)
##
## Coefficient(s):
## a0 a1
## 0.280 0.963
## Sec 9.4: Further Reading
## Spatial statistics
## Other time series models and packages
Box.test(resid(lm(detrendRain ~ detrendSOI, data = xbomsoi)),
type="Ljung-Box", lag=20)
##
## Box-Ljung test
##
## data: resid(lm(detrendRain ~ detrendSOI, data = xbomsoi))
## X-squared = 29.68, df = 20, p-value = 0.07512
## xx <- matrix(x, ncol=1000)
library(tseries)
data(ice.river)
river1 <- diff(log(ice.river[, 1]))
