# Time Series Analysis and Its Applications: With R Examples

## Second Edition

This is the site for the second edition of the text and is no longer maintained.

Follow this link if you're looking for the site of the 4th edition
.

# R Code Used in the Text  (Chapters 1-5)

Everything you see in a box below is R code. You can copy-and-paste any line (or multiple lines) into R. Make sure the data files are in the mydata directory (or change the code accordingly). R code for Chapter 6 has its own page and we also have some code for Chapter 7. We have improved, corrected or streamlined some of the code below, so there are a few cases where the code here differs slightly from what is in the text.

## CHAPTER 1

Example 1.1
``` jj = ts(scan("/mydata/jj.dat"), start=1960, frequency=4)
plot(jj, ylab="Quarterly Earnings per Share")
```

Example 1.2
``` globtemp = ts(scan("/mydata/globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900)     # use data from 1990
plot(gtemp, type="o", ylab="Global Temperature Deviations")
```

Example 1.3
``` plot(speech <- ts(scan("/mydata/speech.dat")), ylab="speech")
```

Example 1.4
``` plot(nyse <- ts(scan("/mydata/nyse.dat")), ylab="NYSE Returns")
```

Example 1.5
``` soi = ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(2,1), mar=c(3, 4, 3, 2))   # sets up the graphics - ?par for info
plot(soi, ylab="", xlab="", main="Southern Oscillation Index")
plot(rec, ylab="", xlab="", main="Recruitment")
```

Example 1.6
``` fmri = read.table("/mydata/fmri.dat")
par(mfrow=c(2,1), mar=c(5, 4, 3, 2))
ts.plot(fmri[,2:5],lty=1:4,ylab="BOLD",xlab="Time (1 pt = 2 sec)",main="Cortex",col=1:4)
ts.plot(fmri[,6:9],lty=1:4,ylab="BOLD",xlab="",main="Thalamus & Cerebellum",col=1:4)

# col=1:4 adds color but if you print monochrome you may want to remove it
```

Example 1.7
``` x = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
par(mfrow=c(2,1))
plot.ts(x[,1], main="Earthquake", ylab="EQ5")
plot.ts(x[,2], main="Explosion", ylab="EXP6")
```

Example 1.9
``` w = rnorm(500,0,1)                  # 500 N(0,1) variates
v = filter(w, sides=2, rep(1,3)/3)  # moving average
par(mfrow=c(2,1))
plot.ts(w, main="white noise")
plot.ts(v, ylim=c(-3,3), main="moving average")
# now try this:
ts.plot(w, v, lty=2:1, col=1:2, lwd=1:2)
```

Example 1.10
``` w = rnorm(550,0,1) #  50 extra to avoid startup problems - the [-(1:50)] on the next line...
x = filter(w, filter=c(1,-.9), method="recursive")[-(1:50)]  #... removes the first 50 values
plot.ts(x)
```

Example 1.11
``` set.seed(154)
w = rnorm(200,0,1)
x = cumsum(w)
wd = w +.2
xd = cumsum(wd)
plot.ts(xd, ylim=c(-5,55))
lines(x)
lines(.2*(1:200), lty="dashed")
```

Example 1.12
``` c = 2*cos(2*pi*(1:500)/50 + .6*pi)
w = rnorm(500,0,1)
par(mfrow=c(3,1), mar=c(3,2,2,1), cex.main=1.5)
plot.ts(c, main=expression(x[t]==2*cos(2*pi*t/50+.6*pi)))
plot.ts(c+w, main=expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,1)))
plot.ts(c+5*w, main=expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,25)))
```

Example 1.24
``` speech = scan("/mydata/speech.dat")
acf(speech,250)
```

Example 1.25*
``` soi = ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(3,1))
acf(soi, 50)
acf(rec, 50)
ccf(soi, rec, 50)
```
* You don't have to read the data in again if you already did it for Example 1.5 and you've saved your workspace. This goes for the rest of the examples. Type objects() to see a list of your objects.

## CHAPTER 2

Example 2.1
``` globtemp = ts(scan("/mydata/globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900)           # use years 1990 to 1997
fit = lm(gtemp~time(gtemp), na.action=NULL)    # regress gtemp on time
#-- note: na.action=NULL is used to retain time series attributes - ?lm for info
summary(fit)                # view the results
plot(gtemp, type="o", ylab="Global Temperature Deviation")   #  plot the series
abline(fit)                 # add regression line to the plot
```

Example 2.2
``` mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
temp = ts(scan("/mydata/temp.dat"), start=1970, frequency=52)
part = ts(scan("/mydata/part.dat"), start=1970, frequency=52)
par(mfrow=c(3,1), mar=c(3,4,3,2))
plot(mort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(temp, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")
x11()                          #  open another graphic device
pairs(cbind(mort, temp, part)) #  scatterplot matrix
temp = temp-mean(temp)
temp2 = temp^2
trend = time(mort)       #  time
fit = lm(mort~ trend + temp + temp2 + part, na.action=NULL)
summary(fit)             #  regression results
anova(lm(mort~cbind(trend, temp, temp2, part))) # pretty ANOVA table
num = length(mort)       #  sample size
AIC(fit)/num - log(2*pi)              #  AIC (corresponds to def 2.1)
AIC(fit, k=log(num))/num - log(2*pi)  #  BIC (corresponds to def 2.3)
(AICc = log(sum(resid(fit)^2)/num)+(num+5)/(num-5-2)) #  AICc  (as in def 2.2)
```

Examples 2.3 and 2.4
``` globtemp = ts(scan("/mydata/globtemp.dat"), start=1856)
gtemp = window(globtemp, start=1900)              # use years 1990 to 1997
fit = lm(gtemp~time(gtemp), na.action=NULL)       # regress gtemp on time
par(mfrow=c(2,1))
plot(resid(fit), type="o", main="detrended")
plot(diff(gtemp), type="o", main="first difference")  #  this and below will do example 2.4
x11()                #  open a new graphic device
par(mfrow=c(3,1))    #  plot ACFs
acf(gtemp, 48)
acf(resid(fit), 48)
acf(diff(gtemp), 48)
```

Example 2.6
You can download a function called lag.plot2 that will do the second part of this example over here.
``` soi = ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
lag.plot(soi, lags=12, layout=c(3,4), diag=FALSE)

par(mfrow=c(3,3), mar=c(2.5, 4, 4, 1)) #  set up plot area
for(h in 0:8){                         #  loop through lags 0-8
plot(lag(soi,-h),rec, main=paste("soi(t-",h,")",sep=""), ylab="rec(t)",xlab="")
}
```

Example 2.8
``` x = 2*cos(2*pi*(1:500)/50 + .6*pi) + rnorm(500,0,5)
I = abs(fft(x)/sqrt(500))^2   #  the periodogram
P = (4/500)*I                 #  the scaled periodogram
f = 0:250/500
plot(f, P[1:251], type="l", xlab="frequency", ylab=" ")
abline(v=seq(0,.5,.02), lty="dotted")
# Note: the dominant period is 50 points or 1/50 = .02 cycles per point
```

Example 2.10
``` mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
ma5 = filter(mort, sides=2, rep(1,5)/5)
ma53 = filter(mort, sides=2, rep(1,53)/53)
plot(mort, type="p", xlab="week", ylab="mortality")
lines(ma5)
lines(ma53)
```

Example 2.11
``` wk = time(mort) - mean(time(mort))   # week (centered) - wk is basically t/52 ...
wk2 = wk^2                           # ... because frequency=52 for mort
wk3 = wk^3
c = cos(2*pi*wk)
s = sin(2*pi*wk)
reg1 = lm(mort~wk + wk2 + wk3, na.action=NULL)
reg2 = lm(mort~wk + wk2 + wk3 + c + s, na.action=NULL)
plot(mort, type="p", ylab="mortality")
lines(fitted(reg1))
lines(fitted(reg2))

# here's the original code where mort isn't a ts object
mort = scan("/mydata/cmort.dat")
t = 1:length(mort)
t2 = t^2
t3 = t^3
c = cos(2*pi*t/52)
s = sin(2*pi*t/52)
fit1 = lm(mort~t + t2 + t3)
fit2 = lm(mort~t + t2 + t3 + c + s)
plot(t, mort)
lines(fit1\$fit)
lines(fit2\$fit)```

Example 2.12
``` plot(mort, type="p", ylab="mortality")
lines(ksmooth(time(mort), mort, "normal", bandwidth=10/52))  #  or try bandwidth=5/52
lines(ksmooth(time(mort), mort, "normal", bandwidth=2))

# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(ksmooth(t, mort, "normal", bandwidth=10))
lines(ksmooth(t, mort, "normal", bandwidth=104))

# Red indicates a change; Fig. 2.15 shows a bandwidth of 10.
```

Example 2.13
``` par(mfrow=c(2,1))
plot(mort,type="p", ylab="mortality", main="nearest neighbor")
lines(supsmu(time(mort), mort, span=.5))
lines(supsmu(time(mort), mort, span=.01))
plot(mort, type="p", ylab="mortality", main="lowess")
lines(lowess(mort, f=.02))
lines(lowess(mort, f=2/3))
```

Example 2.14
``` plot(mort, type="p", ylab="mortality")
lines(smooth.spline(time(mort), mort))
lines(smooth.spline(time(mort), mort, spar=1))

# - original code when mort wasn't a time series with frequency=52
plot(t, mort)
lines(smooth.spline(t, mort, spar=.0000001))
lines(smooth.spline(t, mort, spar=.1))
```

Example 2.15
``` par(mfrow=c(2,1))
plot(temp, mort, main="lowess")
lines(lowess(temp,mort))
plot(temp, mort, main="smoothing splines")
lines(smooth.spline(temp,mort))
```

## CHAPTER 3

Example 3.1
``` par(mfrow=c(2,1))
plot(arima.sim(list(order=c(1,0,0), ar=.9), n=100),
ylab="x",main=(expression(AR(1)~~~phi==+.9)))      # ~ is a space and == is equal
plot(arima.sim(list(order=c(1,0,0), ar=-.9), n=100),
ylab="x",main=(expression(AR(1)~~~phi==-.9)))
```

Example 3.3
``` par(mfrow=c(2,1))
plot(arima.sim(list(order=c(0,0,1), ma=.5), n=100),
ylab="x",main=(expression(MA(1)~~~theta==+.5)))
plot(arima.sim(list(order=c(0,0,1), ma=-.5), n=100),
ylab="x",main=(expression(MA(1)~~~theta==-.5)))
```

Example 3.9
``` set.seed(5)
ar2 = arima.sim(list(order=c(2,0,0), ar=c(1.5,-.75)), n = 144)
t = (1:144)/12
plot(t, ar2, type="l", xlab="Time (one unit = 12 points)")
abline(v=0:12, lty="dotted")
```

Example 3.10
``` ARMAtoMA(ar=.9, ma=.5, 50)       #  for a list
plot(ARMAtoMA(ar=.9, ma=.5, 50)) #  for a graph
```

Example 3.14
``` ar2.acf = ARMAacf(ar=c(1.5,-.75), ma=0, 24)
ar2.pacf = ARMAacf(ar=c(1.5,-.75), ma=0, 24, pacf=TRUE)
par(mfrow=c(1,2))
plot(ar2.acf, type="h", xlab="lag")
abline(h=0)
plot(ar2.pacf, type="h", xlab="lag")
abline(h=0)
```

Example 3.16
``` rec = ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(2,1)); acf(rec, 48); pacf(rec, 48)
(fit = ar.ols(rec,order=2,demean=FALSE,intercept=TRUE)) #  estimates
fit\$asy.se    #  standard errors
```

Example 3.26
``` rec.yw = ar.yw(rec, order=2)
rec.yw\$x.mean
rec.yw\$ar
sqrt(diag(rec.yw\$asy.var.coef))
rec.yw\$var.pred
U = rec.pr\$pred + rec.pr\$se
L = rec.pr\$pred - rec.pr\$se
minx = min(rec,L)
maxx = max(rec,U)
ts.plot(rec, fore\$pred, xlim=c(1980,1990), ylim=c(minx,maxx))
lines(rec.pr\$pred, col="red", type="o")
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")
```

Example 3.29
``` rec.mle = ar.mle(rec, order=2)
rec.mle\$x.mean
rec.mle\$ar
sqrt(diag(rec.mle\$asy.var.coef))
rec.mle\$var.pred
```

Example 3.33
``` x = scan("/mydata/ar1boot.dat")  #  see note below
m = mean(x)                      #  estimate of mu
fit = ar.yw(x, order=1)
phi = fit\$ar                     #  estimate of phi
nboot = 200              #  number of bootstrap replicates
resids = fit\$resid[-1]   #  the first resid is NA (so remove it)
x.star = x               #  initialize x.star
phi.star = matrix(0, nboot, 1)
for (i in 1:nboot) {
resid.star = sample(resids, replace=TRUE)
for (t in 1:99){
x.star[t+1] = m + phi*(x.star[t]-m) + resid.star[t]
}
phi.star[i] = ar.yw(x.star, order=1)\$ar
}
hist(phi.star, prob=TRUE)   # histogram of bootstrapped φs

#  If you want to simulate your own data, use the following commands:
e = rexp(150, rate = .5)
u = runif(150,-1,1)
de = e*sign(u)
x = 50 + arima.sim(n = 100, list(ar =.95), innov = de, n.start = 50)
```

Example 3.35
``` gnp96 = read.table("/mydata/gnp96.dat")
gnp = ts(gnp96[,2], start=1947, frequency=4)
plot(gnp)
acf(gnp, 50)
gnpgr = diff(log(gnp))   # growth rate
plot.ts(gnpgr)
par(mfrow=c(2,1))
acf(gnpgr, 24)
pacf(gnpgr, 24)
# Fit ARIMA and view results:
(gnpgr.ar = arima(gnpgr, order = c(1, 0, 0)))  # potential problem here (see R Issue 1)
(gnpgr.ma = arima(gnpgr, order = c(0, 0, 2)))
ARMAtoMA(ar=.35, ma=0, 10) # prints psi-weights
```

Example 3.36
``` tsdiag(gnpgr.ma, gof.lag=20)
```

Example 3.37
``` varve = scan("/mydata/varve.dat")
(varve.ma = arima(log(varve), order = c(0, 1, 1)))
tsdiag(varve.ma)
(varve.arma = arima(log(varve), order = c(1, 1, 1)))
tsdiag(varve.arma, gof.lag=20)
```

Example 3.39
``` n = length(gnpgr)             # sample size
kma = length(gnpgr.ma\$coef)   # number of parameters in ma model
sma=gnpgr.ma\$sigma2           # mle of sigma^2
kar = length(gnpgr.ar\$coef)   # number of parameters in ar model
sar=gnpgr.ar\$sigma2           # mle of sigma^2

# AIC                         Returned Value
log(sma) + (n+2*kma)/n        # MA2 -8.298
log(sar) + (n+2*kar)/n        # AR1 -8.294
# AICc
log(sma) + (n+kma)/(n-kma-2)  # MA2 -8.288
log(sar) + (n+kar)/(n-kar-2)  # AR1 -8.285
# BIC
log(sma) + kma*log(n)/n       # MA2 -9.252
log(sar) + kar*log(n)/n       # AR1 -9.264

##  Note: You can use sarima.R to get the above
##  (remember to source the code first):
sarima(gnpgr,0,0,2)  #  for the MA(2)
sarima(gnpgr,1,0,0)  #  for the AR(1)
```

Example 3.41
``` phi = c(rep(0,11),.8)
acf = ARMAacf(ar=phi, ma=-.5, 50)
pacf = ARMAacf(ar=phi, ma=-.5, 50, pacf=TRUE)
par(mfrow=c(1,2))
plot(acf, type="h", xlab="lag")
abline(h=0)
plot(pacf, type="h", xlab="lag")
abline(h=0)
```

Example 3.43
``` prod=ts(scan("/mydata/prod.dat"), start=1948, frequency=12)
par(mfrow=c(2,1)) #  (P)ACF of data
acf(prod, 48)
pacf(prod, 48)
par(mfrow=c(2,1)) #  (P)ACF of d1 data
acf(diff(prod), 48)
pacf(diff(prod), 48)
par(mfrow=c(2,1)) #  (P)ACF of d1-d12 data
acf(diff(diff(prod),12), 48)
pacf(diff(diff(prod),12), 48)
###  fit model (iii)
prod.fit3 = arima(prod, order=c(1,1,1), seasonal=list(order=c(2,1,1), period=12))
prod.fit3 #  to view the results
tsdiag(prod.fit3, gof.lag=48) #  diagnostics
###  forecasts for the final model
U = prod.pr\$pred + 2*prod.pr\$se
L = prod.pr\$pred - 2*prod.pr\$se
ts.plot(prod,prodpr\$pred, col=1:2, type="o", ylim=c(105,175), xlim=c(1975,1980))
lines(U, col="blue", lty="dashed")
lines(L, col="blue", lty="dashed")

##  Note: This can be done using acf2.R, sarima.R
##        and sarima.for.R (over here) as follows
##        (remember to source the code first):
prod=scan("/mydata/prod.dat")
acf2(prod,48)
acf2(diff(prod), 48)
acf2(diff(diff(prod),12), 48)
###  fit model (iii)
sarima(prod,1,1,1,2,1,1,12)
###  forecasts for the final model
sarima.for(prod,12,1,1,1,2,1,1,12)
```

## CHAPTER 4

Example 4.1
``` t = 1:100
x1 = 2*cos(2*pi*t*6/100) + 3*sin(2*pi*t*6/100)
x2 = 4*cos(2*pi*t*10/100) + 5*sin(2*pi*t*10/100)
x3 = 6*cos(2*pi*t*40/100) + 7*sin(2*pi*t*40/100)
x = x1 + x2 + x3
par(mfrow=c(2,2))
plot.ts(x1, ylim=c(-10,10), main = expression(omega==6/100~~~A^2==13))
plot.ts(x2, ylim=c(-10,10), main = expression(omega==10/100~~~A^2==41))
plot.ts(x3, ylim=c(-10,10), main = expression(omega==40/100~~~A^2==85))
plot.ts(x, ylim=c(-16,16),main="sum")
```

Example 4.2
``` P = abs(2*fft(x)/100)^2
f = 0:50/100
plot(f, P[1:51], type="o", xlab="frequency", ylab="periodogram")
```

Example 4.7
``` x = c(1,2,3,2,1)
c1 = cos(2*pi*1:5*1/5)
s1 = sin(2*pi*1:5*1/5)
c2 = cos(2*pi*1:5*2/5)
s2 = sin(2*pi*1:5*2/5)
omega1 = cbind(c1, s1)
omega2 = cbind(c2, s2)
anova(lm(x~omega1+omega2))    # ANOVA Table
abs(fft(x))^2/5               # the periodogram (as a comparison)
```

Example 4.9
``` soi = scan("/mydata/soi.dat")
rec = scan("/mydata/recruit.dat")
par(mfrow=c(2,1))
soi.per = spec.pgram(soi, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=1/48, lty="dotted")
rec.per = spec.pgram(rec, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=1/48, lty="dotted")
soi.per\$spec[40] #  soi pgram at freq 1/12 = 40/480
soi.per\$spec[10] #  soi pgram at freq 1/48 = 10/480
#  -- conf intervals -- #  returned value:
U = qchisq(.025,2)     #  0.05063562
L = qchisq(.975,2)     #  7.377759
2*soi.per\$spec[10]/L   #  0.1747835
2*soi.per\$spec[10]/U   #  25.46648
2*soi.per\$spec[40]/L   #  3.162688
2*soi.per\$spec[40]/U   #  460.813
# -- replace soi with rec above to get recruit values

#######################################################################################
###- Repeat example using "frequency=12" for monthly data and note the differences -###
###-  basically, "everything" is in multiples of Δ = 1/12                          -###
#######################################################################################
soi=ts(scan("/mydata/soi.dat"), start=1950, frequency=12)
rec=ts(scan("/mydata/recruit.dat"), start=1950, frequency=12)
par(mfrow=c(2,1))
soi.per = spec.pgram(soi, taper=0, log="no")
abline(v=1, lty="dotted")
abline(v=1/4, lty="dotted")
rec.per = spec.pgram(rec, taper=0, log="no")
abline(v=1, lty="dotted")
abline(v=1/4, lty="dotted")
soi.per\$spec[40] # soi pgram at freq 1/12 = 40/480
soi.per\$spec[10] # soi pgram at freq 1/48 = 10/480
# -- conf intervals -- # returned value:
U = qchisq(.025,2)     # 0.05063562
L = qchisq(.975,2)     # 7.377759
2*soi.per\$spec[10]/L   # 0.01456530
2*soi.per\$spec[10]/U   # 2.122207
2*soi.per\$spec[40]/L   # 0.2635573
2*soi.per\$spec[40]/U   # 38.40108
```

Example 4.10
``` par(mfrow=c(2,1))
k = kernel("daniell",4)
soi.ave = spec.pgram(soi, k, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=2/12, lty="dotted")
abline(v=3/12, lty="dotted")
abline(v=1/48, lty="dotted")
# -- Repeat 5 lines above using rec in place of soi
soi.ave\$bandwidth           #  reported bandwidth
soi.ave\$bandwidth*sqrt(12)  #  Bw

df = soi.ave\$df     #  df = 16.9875 (returned values)
U = qchisq(.025,df) #  U = 7.555916
L = qchisq(.975,df) #  L = 30.17425
soi.ave\$spec[10]    #  0.5942431
soi.ave\$spec[40]    #  1.428959
#  -- intervals --
df*soi.ave\$spec[10]/L #  0.334547
df*soi.ave\$spec[10]/U #  1.336000
df*soi.ave\$spec[40]/L #  0.8044755
df*soi.ave\$spec[40]/U #  3.212641
# -- repeat above commands with soi replaced by rec
```

Example 4.11
``` kernel("modified.daniell", c(3,3))          # for a list
plot(kernel("modified.daniell", c(3,3)))    # for a graph
par(mfrow=c(2,1))
k = kernel("modified.daniell", c(3,3))
soi.smo = spec.pgram(soi, k, taper=0, log="no")
abline(v=1/12, lty="dotted")
abline(v=1/48, lty="dotted")
# -- Repeat above 3 lines with rec replacing soi
df = soi.smo\$df           #  df=17.42618
Lh = 1/sum(k[-k\$m:k\$m]^2) #  Lh=9.232413
Bw = Lh/480               #  Bw=0.01923419

#  An easier way to obtain soi.smo (you can add log="no"):
soi.smo = spectrum(soi, spans=c(7,7), taper=0)
```

Example 4.12
``` par(mfrow=c(3,1))
spectrum(soi, spans=c(7,7), taper=0, main="No Taper")
abline(v=1/12,lty="dashed")
abline(v=1/48,lty="dashed")
spectrum(soi, spans=c(7,7), main="10% Taper")
abline(v=1/12,lty="dashed")
abline(v=1/48,lty="dashed")
spectrum(soi, spans=c(7,7), taper=.5, main="50% Taper")
abline(v=1/12,lty="dashed")
abline(v=1/48,lty="dashed")
```

Example 4.13
``` x = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
eqP = x[1:1024, 1]; eqS = x[1025:2048, 1]
exP = x[1:1024, 2]; exS = x[1025:2048, 2]
par(mfrow=c(2,2))
eqPs=spectrum(eqP, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.04))
eqSs=spectrum(eqS, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.4))
exPs=spectrum(exP, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.04))
exSs=spectrum(exS, spans=c(21,21), log="no", xlim=c(0,.25), ylim=c(0,.4))
exSs\$df
```

Example 4.16
``` x = ts(cbind(soi,rec))
s = spec.pgram(x, kernel("daniell",9), taper=0)
s\$df                     # df = 35.8625
f = qf(.999, 2, s\$df-2)  # f = 8.529792
c = f/(18+f)             # c = 0.3188779
plot(s, plot.type = "coh", ci.lty = 2)
abline(h = c)
```

Example 4.17
``` par(mfrow=c(3,1))
plot.ts(soi)         # the data
plot.ts(diff(soi))   # first difference
k = kernel("modified.daniell", 6)  #-- 12 month filter
soif = kernapply(soi,k)
plot.ts(soif)
x11()                #  open new graphics device
spectrum(soif, spans=9, log="no")  #-- spectral analysis
abline(v=1/52, lty="dotted")
windows()
w = seq(0,.5, length=1000)         #-- frequency response
FR = abs(1-exp(2i*pi*w))^2
plot(w, FR, type="l")
```

Example 4.19
``` spec.ar(soi, log="no")           # plot min AIC spectrum
abline(v=1/52, lty="dotted")     # locate El Niño period
abline(v=1/12, lty="dotted")     # locate yearly period
soi.ar = ar(soi, order.max=30)   # obtain AICs **(see below)
plot(0:30, soi.ar\$aic, type="l") # plot AICs
soi.ar                           # view results

#** Here's a nicer display with the three information criteria:
n = length(soi)
AIC = matrix(0,30,1) -> AICc -> BIC -> sig2
for (k in 1:30){
sigma2=ar(soi, order=k, aic=FALSE, method="burg", var.method=2)\$var.pred
#-- ?ar details suggest the above is the mle of σ2
BIC[k]=log(sigma2)+(k*log(n)/n)
AICc[k]=log(sigma2)+((n+k)/(n-k-2))
AIC[k]=log(sigma2)+((n+2*k)/n)
sig2[k]=sigma2
}
IC=cbind(AIC, AICc, BIC+1)       # BIC+1 reduces white space in the plot
ts.plot(IC,type="o",xlab="order",ylab="AIC / BIC",main="Information Criteria",col=c(1,2,4))
text(12.5,-1.55,"AIC", col=1)
text(14.7,-1.48,"AICc", col=2)
text(15.5,-1.39,"BIC", col=4)
abline(v=1:30, lty=3)```

Example 4.20
``` eqexp = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
ex = eqexp[,2]                   # the explosion series
##  -- dynamic spectral analysis -- ##
nobs = length(ex)                # number of observations
wsize = 256                      # window size
overlap = 128                    # overlap
ovr = wsize-overlap
nseg = floor(nobs/ovr)-1;        # number of segments
krnl = kernel("daniell", c(1,1)) # kernel
ex.spec = matrix(0,wsize/2,nseg)
for (k in 1:nseg) {
a = ovr*(k-1)+1
b = wsize+ovr*(k-1)
ex.spec[,k]=spectrum(ex[a:b],krnl,taper=.5,plot=FALSE)\$spec
}
##  --plot results -- ##
x = seq(0, .5, len = nrow(ex.spec))
y = seq(0, ovr*nseg, len = ncol(ex.spec))
persp(x, y, ex.spec, zlab="Power", xlab="frequency",
ylab="time", ticktype = "detailed", theta=25, d=2)
```

Example 4.21
```# -- uses the waveslim package --#
library(waveslim)
eqexp = matrix(scan("/mydata/eq5exp6.dat"), ncol=2)
eq = eqexp[,1] #  the earthquake series
eq.dwt = dwt(eq, n.levels=6)
#  -- plot the dwt and calculate TP -- #
TP = matrix(0,7,1)
par(mfcol=c(7,1), pty="m", mar=c(3,4,2,2))
for(i in 1:6){
plot.ts(up.sample(eq.dwt[[i]], 2^i), type="h", axes=FALSE,
ylab=names(eq.dwt)[i])
abline(h=0)
axis(side=2)
TP[i]=sum(eq.dwt[[i]]^2)
}
plot.ts(up.sample(eq.dwt[[7]], 2^6), type="h", axes=FALSE,
ylab=names(eq.dwt)[7])
abline(h=0)
axis(side=2)
axis(side=1)
TP[7]=sum(eq.dwt[[7]]^2)
TP/sum(eq^2) #  the energy distribution
#  To plot the ordinates of the wavelet transforms on the
#  same scale, include a command like ylim=c(-1.5,1.5) in
#  each plot.ts() command.
```

Example 4.22
```# -- uses the waveslim package --#
# -- assumes eq is retained from previous example --#
library(waveslim)
eq.dwt = dwt(eq, n.levels=6)
eq.trsh = universal.thresh(eq.dwt, hard=FALSE)
eq.smo = idwt(eq.trsh)
par(mfrow=c(3,1))
plot.ts(eq, ylab="Earthquake", ylim=c(-.5,.5))
plot.ts(eq.smo,ylab="Smoothed Earthquake ",ylim=c(-.5,.5))
plot.ts(eq-eq.smo, ylab="Noise", ylim=c(-.5,.5))
```

Examples 4.23 and 4.24   Doug Wiens, University of Alberta, kindly contributed R code for these examples.
```Example 4.23
Example 4.24
```

## CHAPTER 5

Example 5.1
```# -- uses the fracdiff package --#
library(fracdiff)
lvarve = log(varve)-mean(log(varve))
varve.fd = fracdiff(lvarve, nar=0, nma=0, M=30)
varve.fd\$d
varve.fd\$stderror.dpq
```

Example 5.3
```# -- uses the tseries package --#
library(tseries)
gnpr = diff(log(gnp96[,2]))
gnpr.ar = ar.mle(gnpr, order=1) #  recall phi1 =.347
y = gnpr.ar\$resid[-1]           #  first resid is NA
arch.y = garch(y, order=c(0,1))
summary(arch.y)
```
NB: The code for Examples 5.3 and 5.4 are slightly different than what is in text because the tseries package changed. Also, it may be possible to fit simultaneous [ARMA mean and GARCH variance] models using garchFit() from the package fSeries. Unfortunately, I couldn't get it to work with this example (not that I tried very hard).

Example 5.4
```# -- uses the tseries package --#
library(tseries)
nyse = scan("/mydata/nyse.dat")
nyse.g = garch(nyse, order=c(1,1))
summary(nyse.g)
u = predict(nyse.g)
plot(800:1000, nyse[800:1000], type="l", xlab="Time", ylab="NYSE Returns")
lines(u[,1], col="blue", lty="dashed")
lines(u[,2], col="blue", lty="dashed")
```

Example 5.7
This code does not appear in the text where the results were obtained using ASTSA. Here we use gls() from the nlme package.
``` library(nlme)
mort = scan("/mydata/cmort.dat")
temp = scan("/mydata/temp.dat")
part = scan("/mydata/part.dat")
temp = temp - mean(temp)
temp2 = temp^2
trend = 1:length(mort)
# -- lm() fit in Example 2.2 indicates AR(2) for the residuals
fit.gls = gls(mort~trend + temp + temp2 + part, correlation=corARMA(p=2), method="ML")
# -- the results (abbreviated)
summary(fit.gls)
Correlation Structure: ARMA(2,0)
Parameter estimate(s):
Phi1      Phi2
0.3848530 0.4326282
Coefficients:   Value Std.Error  t-value p-value
(Intercept)  87.63747 2.7327956 32.06880  0.0000
trend        -0.02915 0.0081705 -3.56781  0.0004
temp         -0.01881 0.0430264 -0.43715  0.6622
temp2         0.01542 0.0020289  7.60244  0.0000
part          0.15437 0.0248071  6.22297  0.0000
Residual standard error: 7.699336
Degrees of freedom: 508 total; 503 residual
# -- check noise for whiteness (* see note below)
w = filter(residuals(fit.gls), filter=c(1,-.3848530, -.4326282), sides=1)
w = w[3:508]    #  first two are NA
Box.test(w, 20, type="Ljung")  #  Ljung-Box-Pierce Statistic
X-squared = 26.793, df = 20, p-value = 0.1412
pchisq(26.793, 18, lower=FALSE)  #  the p-value (they are resids from an ar2 fit)
0.08296072
```
* The errors from the gls fit, say e(t), are estimated to be of the form e(t) = .3848530 e(t-1) + .4326282 e(t-2) + w(t). The filter statement gets at the w(t), which should be white.

Examples 5.9 and 5.10
```mort = ts(scan("/mydata/cmort.dat"),start=1970, frequency=52)
temp = ts(scan("/mydata/temp.dat"),start=1970, frequency=52)
part = ts(scan("/mydata/part.dat"),start=1970, frequency=52)
dmort = resid(lm(mort~time(mort)))  # detrend series first
dtemp = resid(lm(temp~time(temp)))
dpart = resid(lm(part~time(part)))

lag.max = 10
BIC = matrix(NA,lag.max,1)
fit = vector("list",lag.max)
n = length(mort)
k = 3      # number of series
for (p in 1:lag.max){
fit[[p]] = ar(cbind(dmort,dtemp,dpart), order=p)
BIC[p] = log(det(fit[[p]]\$var.pred)) + (k^2*p*log(n)/n)
}

BIC              # Example 5:10 - view the BICs, note p=2 is min, and
fit[[2]]         # view the parameter estimates for the min BIC model- VAR(2)

#-- for Example 5.9, look at fit[[1]]
#-- if you look at
fit[[lag.max]]\$aic
#-- you will see all the AICs and you'll notice the min is at p=8
```