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
 rec.pr = predict(rec.yw, n.ahead=24)
 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             
 prod.pr = predict(prod.fit3, n.ahead=12)
 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)
 gnp96 = read.table("/mydata/gnp96.dat")
 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