Time Series Analysis and Its Applications: With R Examples

Third Edition

• |
• Home
• |
 

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

Click this line to go to the site for the 4th edition.

Code Used in the Third Edition Examples

Below is the code used for each numerical example in the text. This stuff won't work unless you have loaded astsa and the data files at the start of the session. If this is your first time here, you might want to read the astsa package notes page for further information.

[+] Chapter 1

Example 1.1

plot(jj, type="o", ylab="Quarterly Earnings per Share")


Example 1.2

plot(gtemp, type="o", ylab="Global Temperature Deviations")


Example 1.3

plot(speech)  


Example 1.4

plot(nyse, ylab="NYSE Returns") 


Example 1.5

par(mfrow = c(2,1))  # set up the graphics
plot(soi, ylab="", xlab="", main="Southern Oscillation Index")
plot(rec, ylab="", xlab="", main="Recruitment") 


Example 1.6

par(mfrow=c(2,1), mar=c(3,2,1,0)+.5, mgp=c(1.6,.6,0))  
ts.plot(fmri1[,2:5], lty=c(1,2,4,5), ylab="BOLD", xlab="", main="Cortex")
ts.plot(fmri1[,6:9], lty=c(1,2,4,5), ylab="BOLD", xlab="", main="Thalamus & Cerebellum")
mtext("Time (1 pt = 2 sec)", side=1, line=2) 


Example 1.7

par(mfrow=c(2,1))
plot(EQ5, main="Earthquake")
plot(EXP6, main="Explosion")


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 (not in text):  
dev.new()  # open a new graphic device
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
x = filter(w, filter=c(1,-.9), method="recursive")[-(1:50)]
plot.ts(x, main="autoregression")


Example 1.11

set.seed(154)  # so you can reproduce the results
w = rnorm(200,0,1)
x = cumsum(w)    
wd = w +.2
xd = cumsum(wd)      
plot.ts(xd, ylim=c(-5,55), main="random walk")
lines(x)
lines(.2*1:200, lty="dashed")


Example 1.12

cs = 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)   # help(par) for info
plot.ts(cs, main = expression(x[t]==2*cos(2*pi*t/50+.6*pi)))
plot.ts(cs + w, main = expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,1)))
plot.ts(cs + 5*w, main = expression(x[t]==2*cos(2*pi*t/50+.6*pi)+N(0,25)))


Example 1.24

acf(speech, 250)


Example 1.25

par(mfrow=c(3,1))
acf(soi, 48, main="Southern Oscillation Index")
acf(rec, 48, main="Recruitment")
ccf(soi, rec, 48, main="SOI vs Recruitment", ylab="CCF")


Example 1.26

persp(1:64, 1:36, soiltemp, phi=30, theta=30, scale=FALSE, expand=4, ticktype="detailed", 
          xlab="rows", ylab="cols", zlab="temperature")
dev.new()          
plot.ts(rowMeans(soiltemp), xlab="row", ylab="Average Temperature")


Example 1.27

fs = abs(fft(soiltemp-mean(soiltemp)))^2/(64*36) # see Ch 4 for info on FFT
cs = Re(fft(fs, inverse=TRUE)/sqrt(64*36))  # ACovF
rs = cs/cs[1,1]                             # ACF

rs2 = cbind(rs[1:41,21:2], rs[1:41,1:21])   #  these lines are just to center
rs3 = rbind(rs2[41:2,], rs2)                #  the 0 lag  

par(mar = c(1,2.5,0,0)+.1)
persp(-40:40, -20:20, rs3, phi=30, theta=30, expand=30, scale="FALSE", ticktype="detailed",
       xlab="row lags", ylab="column lags", zlab="ACF")

[+] Chapter 2

Example 2.1

summary(fit <- lm(gtemp~time(gtemp)))  # regress gtemp on time
plot(gtemp, type="o", ylab="Global Temperature Deviation")
abline(fit)  # add regression line to the plot               


Example 2.2

par(mfrow=c(3,1))
plot(cmort, main="Cardiovascular Mortality", xlab="", ylab="")
plot(tempr, main="Temperature", xlab="", ylab="")
plot(part, main="Particulates", xlab="", ylab="")

dev.new()  
pairs(cbind(Mortality=cmort, Temperature=tempr, Particulates=part))

#  Regression
temp  = tempr-mean(tempr)  # center temperature    
temp2 = temp^2             # square it  
trend = time(cmort)        # time

fit = lm(cmort~ trend + temp + temp2 + part, na.action=NULL)
            
summary(fit)       # regression results
summary(aov(fit))  # ANOVA table   (compare to next line)
summary(aov(lm(cmort~cbind(trend, temp, temp2, part)))) # Table 2.1

num = length(cmort)                                     # sample size
AIC(fit)/num - log(2*pi)                                # AIC 
BIC(fit)/num - log(2*pi)                                # BIC 
# AIC(fit, k=log(num))/num - log(2*pi)                  # BIC (alt method)    
(AICc = log(sum(resid(fit)^2)/num) + (num+5)/(num-5-2)) # AICc


Examples 2.3

fish = ts.intersect(rec, soiL6=lag(soi,-6), dframe=TRUE)   
summary(fit <- lm(rec~soiL6, data=fish, na.action=NULL))
plot(fish$rec)  # plot the data and the fitted values (not shown in text) 
lines(fitted(fit), col=2)


Examples 2.4 and 2.5

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")

dev.new()
par(mfrow=c(3,1))     # plot ACFs
acf(gtemp, 48, main="gtemp")
acf(resid(fit), 48, main="detrended")
acf(diff(gtemp), 48, main="first difference")


Example 2.6

par(mfrow=c(2,1))
plot(varve, main="varve", ylab="")
plot(log(varve), main="log(varve)", ylab="" )


Example 2.7

lag1.plot(soi, 12)
dev.new()
lag2.plot(soi, rec, 8)


Example 2.8

set.seed(1000)  # so you can reproduce these results
x = 2*cos(2*pi*1:500/50 + .6*pi) + rnorm(500,0,5)
z1 = cos(2*pi*1:500/50)  
z2 = sin(2*pi*1:500/50)

summary(fit <- lm(x~0+z1+z2))  # zero to exclude the intercept
plot.ts(x, lty="dashed")
lines(fitted(fit), lwd=2)


Example 2.9

I = abs(fft(x))^2/500  # the periodogram
P = (4/500)*I[1:250]   # the scaled periodogram
f = 0:249/500          # frequencies
plot(f, P, type="l", xlab="Frequency", ylab="Scaled Periodogram")


Example 2.10

ma5 = filter(cmort, sides=2, rep(1,5)/5)
ma53 = filter(cmort, sides=2, rep(1,53)/53)
plot(cmort, type="p", ylab="mortality")
lines(ma5)
lines(ma53)


Example 2.11

wk = time(cmort) - mean(time(cmort))  # wk is essentially t/52 centered at zero  
wk2 = wk^2 
wk3 = wk^3
cs = cos(2*pi*wk)  
sn = sin(2*pi*wk)
reg1 = lm(cmort~wk + wk2 + wk3, na.action=NULL)
reg2 = lm(cmort~wk + wk2 + wk3 + cs + sn, na.action=NULL)
plot(cmort, type="p", ylab="mortality")
lines(fitted(reg1))
lines(fitted(reg2))


Example 2.12

plot(cmort, type="p", ylab="mortality")
lines(ksmooth(time(cmort), cmort, "normal", bandwidth=5/52))
lines(ksmooth(time(cmort), cmort, "normal", bandwidth=2))


Example 2.13

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


Example 2.14

plot(cmort, type="p", ylab="mortality")
lines(smooth.spline(time(cmort), cmort))          
lines(smooth.spline(time(cmort), cmort, spar=1))  


Example 2.15

par(mfrow=c(2,1), mar=c(3,2,1,0)+.5, mgp=c(1.6,.6,0))
plot(tempr, cmort, main="lowess", xlab="Temperature", ylab="Mortality")
lines(lowess(tempr,cmort))
plot(tempr, cmort, main="smoothing splines", xlab="Temperature", ylab="Mortality")
lines(smooth.spline(tempr, cmort))

[+] Chapter 3

Example 3.1

par(mfrow=c(2,1))                         
# in the expressions below, ~ 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))) 
plot(arima.sim(list(order=c(1,0,0), ar=-.9), n=100), ylab="x", main=(expression(AR(1)~~~phi==-.9))) 


Example 3.4

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.10

z = c(1,-1.5,.75)    # coefficients of the polynomial
(a = polyroot(z)[1]) # = 1+0.57735i,  print one root which is 1 + i 1/sqrt(3)
arg = Arg(a)/(2*pi)  # arg in cycles/pt  
1/arg                # = 12,  the period

set.seed(8675309)    # Jenny, I got your number
ar2 = arima.sim(list(order=c(2,0,0), ar=c(1.5,-.75)), n = 144)
plot(ar2, axes=FALSE, xlab="Time")
axis(2); axis(1, at=seq(0,144,by=12)); box()  # work the plot machine
abline(v=seq(0,144,by=12), lty=2)

ACF = ARMAacf(ar=c(1.5,-.75), ma=0, 50)
plot(ACF, type="h", xlab="lag")
abline(h=0)


Example 3.11

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


Example 3.15

ar2.acf = ARMAacf(ar=c(1.5,-.75), ma=0, 24)[-1]
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.17

acf2(rec, 48)     # will produce values and a graphic 
(regr = ar.ols(rec, order=2, demean=F, intercept=TRUE))  # regression
regr$asy.se.coef  # standard errors                             


Example 3.24

regr = ar.ols(rec, order=2, demean=FALSE, intercept=TRUE)
fore = predict(regr, n.ahead=24)
ts.plot(rec, fore$pred, col=1:2, xlim=c(1980,1990)) 
lines(fore$pred, type="p", col=2)
lines(fore$pred + fore$se, lty="dashed", col=4)
lines(fore$pred - fore$se, lty="dashed", col=4)


Example 3.27

rec.yw = ar.yw(rec, order=2)
rec.yw$x.mean  # = 62.26278 (mean estimate)
rec.yw$ar      # = 1.3315874, -.4445447  (parameter estimates)
sqrt(diag(rec.yw$asy.var.coef))  # = .04222637, .04222637  (standard errors)
rec.yw$var.pred  # = 94.79912 (error variance estimate)

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, rec.pr$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.28

set.seed(2)
ma1 = arima.sim(list(order = c(0,0,1), ma = 0.9), n = 50)
acf(ma1, plot=FALSE)[1]  # = .507 (lag 1 sample ACF)


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.35

set.seed(111)

# simulate from model
phi.yw = rep(NA, 1000)
for (i in 1:1000){
 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)
 phi.yw[i] = ar.yw(x, order=1)$ar 
}

# plot simulated distribution of phi 
hist(phi.yw, prob=TRUE, main="")
lines(density(phi.yw, bw=.015))

# bootstrap 
x = ar1boot 
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 
x.star = x    # initialize x.star
phi.star.yw = rep(NA, nboot)
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.yw[i] = ar.yw(x.star, order=1)$ar 
}

# bootstrap distribution of phi
dev.new()
hist(phi.star.yw, 10, main="", prob=TRUE, ylim=c(0,14), xlim=c(.75,1.05))
lines(density(phi.star.yw, bw=.02))
u = seq(.75, 1.05, by=.001) 
lines(u, dnorm(u,mean=.96,sd=.03), lty="dashed", lwd=2)


Example 3.37

set.seed(666)   # not that 666
x = arima.sim(list(order = c(0,1,1), ma = -0.8), n = 100)
(x.ima = HoltWinters(x, beta=FALSE, gamma=FALSE))  # α is 1-λ here
plot(x.ima)


Example 3.38

plot(gnp)
acf2(gnp, 50)           
gnpgr = diff(log(gnp))  # growth rate
plot(gnpgr)
acf2(gnpgr, 24)
sarima(gnpgr, 1, 0, 0)  # AR(1)
sarima(gnpgr, 0, 0, 2)  # MA(2) 
ARMAtoMA(ar=.35, ma=0, 10)  # prints psi-weights


Example 3.40

sarima(log(varve), 0, 1, 1, no.constant=TRUE)   # ARIMA(0,1,1)
dev.new()
sarima(log(varve), 1, 1, 1, no.constant=TRUE)   # ARIMA(1,1,1)


Example 3.44

phi = c(rep(0,11),.8)
ACF = ARMAacf(ar=phi, ma=-.5, 50)[-1]  
PACF = ARMAacf(ar=phi, ma=-.5, 50, pacf=TRUE) 
par(mfrow=c(1,2))
plot(ACF,  type="h", xlab="lag", ylim=c(-.4,.8))
abline(h=0)
plot(PACF, type="h", xlab="lag", ylim=c(-.4,.8))
abline(h=0)


Example 3.46

acf2(prodn,48)
dev.new()
acf2(diff(prodn), 48)
dev.new()
acf2(diff(diff(prodn),12), 48)
dev.new()
sarima(prodn,2,1,1,0,1,3,12)        # fit model (ii)
dev.new()
sarima.for(prodn,12,2,1,1,0,1,3,12) # 12 month ahead forecasts 

[+] Chapter 4

Example 4.1

x1 = 2*cos(2*pi*1:100*6/100)  + 3*sin(2*pi*1:100*6/100)
x2 = 4*cos(2*pi*1:100*10/100) + 5*sin(2*pi*1:100*10/100)
x3 = 6*cos(2*pi*1:100*40/100) + 7*sin(2*pi*1:100*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
Fr = 0:99/100                    
plot(Fr, P, type="o", xlab="frequency", ylab="periodogram")


Example 4.6

par(mfrow=c(3,1))       
arma.spec(log="no", main="White Noise")
arma.spec(ma=.5, log="no", main="Moving Average")  
arma.spec(ar=c(1,-.9), log="no", main="Autoregression")


Example 4.8

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 check)


Example 4.10

par(mfrow=c(2,1))      
soi.per = mvspec(soi, log="no")                 # using astsa package
# soi.per = spec.pgram(soi, taper=0, log="no")  # using stats package
abline(v=1/4, lty="dotted")
rec.per = mvspec(rec, log="no")
# rec.per = spec.pgram(rec, taper=0, log="no") 
abline(v=1/4, lty="dotted")

soi.per$spec[40]  # 0.97223;  soi pgram at freq 1/12 = 40/480
soi.per$spec[10]  # 0.05372;  soi pgram at freq 1/48 = 10/480

# conf intervals -  returned value:
U = qchisq(.025,2)    # 0.05063  
L = qchisq(.975,2)    # 7.37775
2*soi.per$spec[10]/L  # 0.01456
2*soi.per$spec[10]/U  # 2.12220
2*soi.per$spec[40]/L  # 0.26355
2*soi.per$spec[40]/U  # 38.40108

# Repeat lines above using rec in place of soi


Example 4.11

par(mfrow=c(2,1))
k = kernel("daniell",4)
soi.ave = mvspec(soi, k, log="no")                  # using astsa package
# soi.ave = spec.pgram(soi, k, taper=0, log="no")   # using stats package
abline(v=c(.25,1,2,3), lty=2)

soi.ave$bandwidth                 # 0.225; bandwidth, time unit = year
soi.ave$bandwidth/frequency(soi)  # 0.01875; bandwidth, time unit = month

df = soi.ave$df         # df = 16.9875
U = qchisq(.025, df)    # U = 7.555916
L = qchisq(.975, df)    # L = 30.17425
soi.ave$spec[10]        # 0.04952026
soi.ave$spec[40]        # 0.1190800

# intervals 
df*soi.ave$spec[10]/L   # 0.02787891
df*soi.ave$spec[10]/U   # 0.1113333
df*soi.ave$spec[40]/L   # 0.06703963
df*soi.ave$spec[40]/U   # 0.2677201

# Repeat above commands with soi replaced by rec, for example:
rec.ave = mvspec(rec, k, log="no")
abline(v=c(.25,1,2,3), lty=2)
# and so on.


Example 4.12

t = seq(0, 1, by=1/200)  # WARNING: using t is bad pRactice because it's reserved- but let's be bad
amps = c(1, .5, .4, .3, .2, .1)
x = matrix(0, 201, 6)
for (j in 1:6) x[,j] = amps[j]*sin(2*pi*t*2*j)
x = ts(cbind(x, rowSums(x)), start=0, deltat=1/200)               
ts.plot(x, lty=c(1:6, 1), lwd=c(rep(1,6), 2), ylab="Sinusoids")
names = c("Fundamental","2nd Harmonic","3rd Harmonic","4th Harmonic","5th Harmonic", 
          "6th Harmonic","Formed Signal")
legend("topright", names, lty=c(1:6, 1), lwd=c(rep(1,6), 2))
rm(t)                    # Redemption


Example 4.13

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 = mvspec(soi, k, log="no")                     # using astsa package
# soi.smo = spec.pgram(soi, k, taper=0, log="no")      # using stats package
  abline(v=1, lty="dotted")    
  abline(v=1/4, lty="dotted")
rec.smo = mvspec(rec, k, log="no")  
# rec.smo = spec.pgram(rec, k, taper=0, log="no")
  abline(v=1, lty="dotted")    
  abline(v=1/4, lty="dotted")

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 units = month
soi.smo$bandwidth         # Bw = 0.2308103, units = year (value returned by mvspec)         
soi.smo$bandwidth/12      # Bw = 0.01923419 units = month (value returned by mvspec)  

# An easier way to obtain soi.smo:
soi.smo = mvspec(soi, spans=c(7,7), log="no")                  # using astsa package
# soi.smo = spectrum(soi, spans=c(7,7), taper=0, log="no")     # using stats package


Example 4.14

s0 = mvspec(soi, spans=c(7,7), plot=FALSE)                # using astsa package
s50 = mvspec(soi, spans=c(7,7), taper=.5, plot=FALSE)     
# s0 = spectrum(soi, spans=c(7,7), taper=0, plot=FALSE)   # using stats package
# s50 = spectrum(soi, spans=c(7,7), taper=.5, plot=FALSE)
plot(s0$freq, s0$spec, log="y", type="l", lty=2, ylab="spectrum", xlab="frequency")  # dashed line
lines(s50$freq, s50$spec)   # solid line
title("Tapering")
legend("topright", c("with you", "without you"), lty=1:2)
arrows(1.5,.025,.7,.012, col="purple")
text(1.7,.027, "leakage", col="purple")


Example 4.15

# this part uses Yule-Walker for the estimates, so there's no likelihood involved in this
# I guess you could call it pseudo-AIC

spaic = spec.ar(soi, log="no", ylim=c(0,.3))             # min AIC spec, order = 15
text(frequency(soi)*1/52, .07, substitute(omega==1/52))  # El Nino Cycle
text(frequency(soi)*1/12, .27, substitute(omega==1/12))  # Yearly Cycle
(soi.ar = ar(soi, order.max=30))                         # estimates and AICs
dev.new() 
plot(1:30, soi.ar$aic[-1], type="o")                     # plot AICs


# Better comparison of pseudo-ICs 
n = length(soi)
AIC = rep(0,30) -> AICc -> BIC
for (k in 1:30){
  fit = ar(soi, order=k, aic=FALSE)
# sigma2 = var(fit$resid, na.rm=TRUE)  # this doesn't seem to work right - use next line 
  sigma2 = fit$var.pred                # with this, all IC pick the order 15 model
  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) }
IC = cbind(AIC, BIC+1)

dev.new()
ts.plot(IC, type="o", xlab="p", ylab="AIC / BIC")
grid()
text(15.2, -1.48, "AIC")
text(15, -1.35, "BIC")


Example 4.18

sr = mvspec(cbind(soi,rec), kernel("daniell",9), plot.type="coh")
# sr = spec.pgram(cbind(soi,rec), kernel("daniell",9), taper=0, plot.type="coh")
sr$df                     # df = 35.8625
f = qf(.999, 2, sr$df-2)  # f = 8.529792
C = f/(18+f)              # C = 0.3188779
abline(h = C)


Example 4.19

par(mfrow=c(3,1))
plot(soi)        # plot data
plot(diff(soi))  # plot first difference
k = kernel("modified.daniell", 6)  # filter weights
plot(soif <- kernapply(soi, k))    # plot 12 month filter 
  
dev.new()                          # open new graphics device
mvspec(soif, spans=9, log="no")    # spectral analysis
# spectrum(soif, spans=9, log="no")  
abline(v=12/52, lty="dashed")

dev.new()
w = seq(0, .5, length=500)         # frequency response
FR = abs(1-exp(2i*pi*w))^2
plot(w, FR, type="l", xlab=expression(omega))


Example 4.21

nobs = length(EXP6)  # 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] = mvspec(EXP6[a:b], krnl, taper=.5, plot=FALSE)$spec 
  }
x = seq(0, 10,  len = nrow(ex.spec)/2)
y = seq(0, ovr*nseg, len = ncol(ex.spec))
z = ex.spec[1:(nrow(ex.spec)/2),]

# below is text version
filled.contour(x,y,log(z),ylab="time",xlab="frequency (Hz)",nlevels=12,col=gray(11:0/11),main="Explosion")

dev.new()  # a nicer version with color
filled.contour(x, y, log(z), ylab="time", xlab="frequency (Hz)", main="Explosion")

dev.new()  # below not shown in text
persp(x,y,z,zlab="Power",xlab="frequency (Hz)",ylab="time",ticktype="detailed",theta=25,d=2,main="Explosion")  


Example 4.22

library(wavethresh)
eq = scale(EQ5)  # standardize the series
ex = scale(EXP6) 

eq.dwt = wd(eq, filter.number=8) 
ex.dwt = wd(ex, filter.number=8)

# plot the wavelet transforms 
par(mfrow = c(1,2))
plot(eq.dwt, main="Earthquake")   
plot(ex.dwt, main="Explosion")  

# total power 
TPe = rep(NA,11)  # for the earthquake series
 for (i in 0:10){TPe[i+1] = sum(accessD(eq.dwt, level=i)^2)}
TotEq = sum(TPe)  # check with sum(eq^2) 
TPx = rep(NA,11)  # for the explosion series
 for (i in 0:10){TPx[i+1] = sum(accessD(ex.dwt, level=i)^2)}
TotEx = sum(TPx)  # check with sum(ex^2) 

# make a nice table
Power = round(cbind(11:1, 0:10, 100*TPe/TotEq, 100*TPx/TotEx), digits=3)
colnames(Power) = c("Level", "Resolution", "EQ(%)", "EXP(%)")
Power


Example 4.23

library(wavethresh)
eq = scale(EQ5)
par(mfrow=c(3,1))
eq.dwt = wd(eq, filter.number=8) 
eq.smo = wr(threshold(eq.dwt, levels=5:10))
ts.plot(eq, main="Earthquake", ylab="Data")
ts.plot(eq.smo, ylab="Signal")
ts.plot(eq-eq.smo, ylab="Resid")


Example 4.24

LagReg(soi, rec, L=15, M=32, threshold=6)
LagReg(rec, soi, L=15, M=32, inverse=TRUE, threshold=.01)


Example 4.25

SigExtract(soi, L=9, M=64, max.freq=.05) 


Example 4.26

per = abs(fft(soiltemp-mean(soiltemp))/sqrt(64*36))^2       
per2 = cbind(per[1:32,18:2], per[1:32,1:18])   # this and line below is just rearranging
per3 = rbind(per2[32:2,], per2)                # results to get 0 frequency in the middle

par(mar=c(1,2.5,0,0)+.1)
persp(-31:31/64, -17:17/36, per3, phi=30, theta=30, expand=.6, ticktype="detailed", xlab="cycles/row", 
       ylab="cycles/column", zlab="Periodogram Ordinate")

[+] Chapter 5

Example 5.1

library(fracdiff)
lvarve = log(varve) - mean(log(varve))
varve.fd = fracdiff(lvarve, nar=0, nma=0, M=30)
varve.fd$d  # = 0.3841688
varve.fd$stderror.dpq  # = 4.589514e-06 (If you believe this, I have a bridge for sale.)

p = rep(1,31)
for (k in 1:30){ p[k+1] = (k-varve.fd$d)*p[k]/(k+1) }
plot(1:30, p[-1], ylab=expression(pi(d)), xlab="Index", type="h", lwd=2)
res.fd = diffseries(log(varve), varve.fd$d)       # frac diff resids            
res.arima = resid(arima(log(varve), order=c(1,1,1))) # arima resids

dev.new()
par(mfrow=c(2,1))  
acf(res.arima, 100, xlim=c(4,97), ylim=c(-.2,.2), main="arima resids")
acf(res.fd, 100, xlim=c(4,97), ylim=c(-.2,.2), main="frac diff resids")


Example 5.2

series = log(varve)  # specify series to be analyzed
d0 = .1              # initial value of d
n.per = nextn(length(series))
m = (n.per)/2  - 1
per = abs(fft(series-mean(series))[-1])^2  # remove 0 freq
per = per/n.per      # R doesn't scale fft by sqrt(n)
g = 4*(sin(pi*((1:m)/n.per))^2)

# Function to calculate -log.likelihood
whit.like = function(d){  
 g.d=g^d
 sig2 = (sum(g.d*per[1:m])/m)
 log.like = m*log(sig2) - d*sum(log(g)) + m
 return(log.like)   
}

# Estimation (?optim for details - output not shown)
(est = optim(d0, whit.like, gr=NULL, method="L-BFGS-B", hessian=TRUE, lower=-.5, upper=.5, 
             control=list(trace=1,REPORT=1)))

# Results  [d.hat = .380, se(dhat) = .028]  
cat("d.hat =", est$par, "se(dhat) = ",1/sqrt(est$hessian),"\n")  
g.dhat = g^est$par
sig2 = sum(g.dhat*per[1:m])/m
cat("sig2hat =",sig2,"\n")  # sig2hat = .229 

u = spec.ar(log(varve), plot=FALSE)  # produces AR(8)       
g = 4*(sin(pi*((1:500)/2000))^2)
fhat = sig2*g^{-est$par}             # long memory spectral estimate          
plot(1:500/2000, log(fhat), type="l", ylab="log(spectrum)", xlab="frequency")
lines(u$freq[1:250], log(u$spec[1:250]), lty="dashed") 
ar.mle(log(varve))                   # to get AR(8) estimates 


Example 5.3

library(tseries)
adf.test(log(varve), k=0)  # DF test
adf.test(log(varve))       # ADF test
pp.test(log(varve))        # PP test


Example 5.4

gnpgr = diff(log(gnp))  # get the returns
sarima(gnpgr, 1, 0, 0)  # fit an AR(1)
acf2(innov^2, 24)       # get (p)acf of the squared residuals
 
library(fGarch)
summary(garchFit(~arma(1,0)+garch(1,0), gnpgr))


Example 5.5

library(fGarch)
nyse = astsa::nyse   
summary(nyse.g <- garchFit(~arma(1,0)+garch(1,1), nyse))
u = nyse.g@sigma.t
plot(window(nyse, start=900, end=1000), ylim=c(-.22,.2), ylab="NYSE Returns")
lines(window(nyse-2*u, start=900, end=1000), lty=2, col=4)
lines(window(nyse+2*u, start=900, end=1000), lty=2, col=4)


Example 5.6

plot(flu, type="c")  # first 3 lines plot the data with months as symbols
Months = c("J","F","M","A","M","J","J","A","S","O","N","D")
points(flu, pch=Months, cex=.8, font=2) 

require(tsDyn)              # load package - install it if you don't have it
# vignette("tsDyn")         # for package details (it's quirky, so you'll need this)

dflu = diff(flu)
lag1.plot(dflu, 1, corr=FALSE)     # see the nonlinearity here
(u = setar(dflu, m=4, thDelay=0))  # fit model and view results (thDelay=0 is lag 1 delay)
BIC(u); AIC(u)                     # if you want to try other models ... m=3 works well too
plot(u)                            # graphics -  ?plot.setar for information


Example 5.7

trend = time(cmort)
temp = tempr - mean(tempr)
temp2 = temp^2

fit = lm(cmort~trend + temp + temp2 + part, na.action=NULL)
acf2(resid(fit), 52)   # implies AR2
sarima(cmort, 2,0,0, xreg=cbind(trend,temp,temp2,part))   # Note change with astsa version 1.3+   


Example 5.8

acf2(soi)
(fit = arima(soi, xreg=time(soi), order=c(1, 0, 0)))       
ar1 = as.numeric(fit$coef[1])    # = 0.5875387
soi.pw = resid(fit)  

rec.d = resid(lm(rec~time(rec), na.action=NULL))  
rec.fil = filter(rec.d, filter=c(1, -ar1), method="conv", sides=1)
ccf(soi.pw, rec.fil, main="", ylab="CCF", na.action=na.omit)


Example 5.9

rec.d = resid(lm(rec~time(rec), na.action=NULL))  
soi.d = resid(lm(soi~time(soi), na.action=NULL))

fish = ts.intersect(rec.d, rec.d1=lag(rec.d,-1), soi.d5=lag(soi,-5), dframe=TRUE)

summary(fish.fit <- lm(rec.d~0+rec.d1+soi.d5, data=fish))

# corrected version below
eta.hat = filter(resid(fish.fit), filter=-om1, method="recur", sides=1)
acf2(eta.hat)  # implies an AR(2)
(eta.fit <- arima(eta.hat, order=c(2,0,0)))


Example 5.10

library(vars)
x = cbind(cmort, tempr, part)
summary(VAR(x, p=1, type="both"))  # "both" fits constant + trend


Example 5.11

# continued from 5.10
VARselect(x, lag.max=10, type="both")
summary(fit <- VAR(x, p=2, type="both"))

acf(resid(fit), 52)
serial.test(fit, lags.pt=12, type="PT.adjusted")

(fit.pr = predict(fit, n.ahead = 24, ci = 0.95))  # 4 weeks ahead

dev.new()
fanchart(fit.pr)  # plot prediction + error

[+] Chapter 6

Example 6.1

blood = cbind(WBC, PLT, HCT)           # skip these first two lines if using astsa
blood = replace(blood, blood==0, NA)   # because blood is already there for you

plot(blood, type="o", pch=19, xlab="day", main="")  


Example 6.2

ts.plot(gtemp, gtemp2, lty=1:2, ylab="Temperature Deviations")


Example 6.5

# generate data 
set.seed(1)  
num = 50
w = rnorm(num+1,0,1)
v = rnorm(num,0,1)
                               
mu = cumsum(w)  # states:  mu[0], mu[1], . . ., mu[50] 
y = mu[-1] + v  # obs:  y[1], . . ., y[50]

# filter and smooth (Ksmooth0 does both)
mu0 = 0; sigma0 = 1;  phi = 1; cQ = 1; cR = 1   
ks = Ksmooth0(num, y, 1, mu0, sigma0, phi, cQ, cR)   

# pictures 
Time = 1:num; par(mfrow=c(3,1))

plot(Time, mu[-1], main="Prediction", ylim=c(-5,10))      
  lines(ks$xp)
  lines(ks$xp+2*sqrt(ks$Pp), lty="dashed", col="blue")
  lines(ks$xp-2*sqrt(ks$Pp), lty="dashed", col="blue")

plot(Time, mu[-1], main="Filter", ylim=c(-5,10))
  lines(ks$xf)
  lines(ks$xf+2*sqrt(ks$Pf), lty="dashed", col="blue")
  lines(ks$xf-2*sqrt(ks$Pf), lty="dashed", col="blue")

plot(Time, mu[-1],  main="Smoother", ylim=c(-5,10))
  lines(ks$xs)
  lines(ks$xs+2*sqrt(ks$Ps), lty="dashed", col="blue")
  lines(ks$xs-2*sqrt(ks$Ps), lty="dashed", col="blue") 

mu[1]; ks$x0n; sqrt(ks$P0n)   # initial value info

# In case you can't see the differences in the 3 figures for this example, try this... 
# Predictor, Filter, Smoother on Same Plot (not shown)
dev.new()
plot(Time, mu[-1])
lines(ks$xp, col=4)
lines(ks$xf, col=3)  
lines(ks$xs, col=2)
abline(v=Time, lty=3, col="lightblue")
names = c("predictor","filter","smoother")
legend("bottomright", names, col=4:2, lty=1, bg="white")


Example 6.6

# Generate Data
set.seed(999)
num = 100
N = num+1
x = arima.sim(n=N, list(ar = .8, sd=1))
y = ts(x[-1] + rnorm(num,0,1))     

# Initial Estimates 
u = ts.intersect(y, lag(y,-1), lag(y,-2)) 
varu = var(u)
coru = cor(u) 
phi = coru[1,3]/coru[1,2]             
q = (1-phi^2)*varu[1,2]/phi   
r = varu[1,1] - q/(1-phi^2) 
(init.par = c(phi, sqrt(q), sqrt(r))) 

# Function to evaluate the likelihood 
Linn=function(para){
  phi = para[1]; sigw = para[2]; sigv = para[3]   
  Sigma0 = (sigw^2)/(1-phi^2); Sigma0[Sigma0<0]=0   
  kf = Kfilter0(num,y,1,mu0=0,Sigma0,phi,sigw,sigv)
  return(kf$like)   
  }

# Estimation  
(est = optim(init.par, Linn, gr=NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1)))      
SE = sqrt(diag(solve(est$hessian)))
cbind(estimate=c(phi=est$par[1],sigw=est$par[2],sigv=est$par[3]), SE)


Example 6.7

# Setup 
y = cbind(gtemp, gtemp2)
num = nrow(y)
input = rep(1,num)
A = array(rep(1,2), dim=c(2,1,num))
mu0 = -.26; Sigma0 = .01;  Phi = 1

# Function to Calculate Likelihood 
Linn=function(para){
 cQ = para[1]      # sigma_w
  cR1 = para[2]    # 11 element of chol(R)
  cR2 = para[3]    # 22 element of chol(R)
  cR12 = para[4]   # 12 element of chol(R)
 cR = matrix(c(cR1,0,cR12,cR2),2)  # put the matrix together
 drift = para[5]
 kf = Kfilter1(num,y,A,mu0,Sigma0,Phi,drift,0,cQ,cR,input)
 return(kf$like) 
 }

# Estimation
init.par = c(.1,.1,.1,0,.05)  # initial values of parameters
(est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1))) 
SE = sqrt(diag(solve(est$hessian))) 

# Summary of estimation  
estimate = est$par; u = cbind(estimate, SE)
rownames(u)=c("sigw","cR11", "cR22", "cR12", "drift"); u  

# Smooth (first set parameters to their final estimates)
cQ=est$par[1]  
 cR1=est$par[2]   
 cR2=est$par[3]   
 cR12=est$par[4]  
cR = matrix(c(cR1,0,cR12,cR2), 2)
(R = t(cR)%*%cR)    #  to view the estimated R matrix
drift = est$par[5]  
ks = Ksmooth1(num,y,A,mu0,Sigma0,Phi,drift,0,cQ,cR,input)  

# Plot 
xsmooth = ts(as.vector(ks$xs), start=1880)
plot(xsmooth, lwd=2, ylim=c(-.5,.8), ylab="Temp Deviations")  
lines(gtemp, col="blue", lty=2)     # color helps here
lines(gtemp2, col="red", lty=2)


Example 6.8

library(nlme)   # loads package nlme 

# Generate data (same as Example 6.6)
set.seed(999); num = 100; N = num+1
x = arima.sim(n=N, list(ar = .8, sd=1))
y = ts(x[-1] + rnorm(num,0,1))     

# Initial Estimates 
u = ts.intersect(y,lag(y,-1),lag(y,-2)) 
varu = var(u); coru = cor(u) 
phi = coru[1,3]/coru[1,2]             
q = (1-phi^2)*varu[1,2]/phi   
r = varu[1,1] - q/(1-phi^2) 
cr = sqrt(r); cq = sqrt(q); mu0 = 0; Sigma0 = 2.8
(em = EM0(num, y, 1, mu0, Sigma0, phi, cq, cr, 75, .00001))   

# Standard Errors  (this uses nlme)
phi = em$Phi; cq = chol(em$Q); cr = chol(em$R)
mu0 = em$mu0; Sigma0 = em$Sigma0
para = c(phi, cq, cr)

# Evaluate likelihood at estimates 
Linn=function(para){
  kf = Kfilter0(num, y, 1, mu0, Sigma0, para[1], para[2], para[3])
  return(kf$like) 
  }
emhess = fdHess(para, function(para) Linn(para))
SE = sqrt(diag(solve(emhess$Hessian)))  

# Display summary of estimation 
estimate = c(para, em$mu0, em$Sigma0); SE = c(SE,NA,NA)
u = cbind(estimate, SE)
rownames(u) = c("phi","sigw","sigv","mu0","Sigma0")
u 


Example 6.9

y = cbind(WBC, PLT, HCT)
num = nrow(y)       
A = array(0, dim=c(3,3,num))  # creates num 3x3 zero matrices
for(k in 1:num) if (y[k,1] > 0) A[,,k]= diag(1,3) 

# Initial values 
mu0 = matrix(0,3,1) 
Sigma0 = diag(c(.1,.1,1),3)
Phi = diag(1,3)
cQ = diag(c(.1,.1,1),3)
cR = diag(c(.1,.1,1),3)  
(em = EM1(num,y,A,mu0,Sigma0,Phi,cQ,cR,100,.001))   # note change with astsa version 1.3+

# Graph smoother
ks = Ksmooth1(num, y, A, em$mu0, em$Sigma0, em$Phi, 0, 0, chol(em$Q), chol(em$R), 0) 
y1s = ks$xs[1,,]
y2s = ks$xs[2,,]
y3s = ks$xs[3,,]
p1 = 2*sqrt(ks$Ps[1,1,])
p2 = 2*sqrt(ks$Ps[2,2,])
p3 = 2*sqrt(ks$Ps[3,3,])

par(mfrow=c(3,1), mar=c(4,4,1,1)+.2)
plot(WBC, type="p", pch=19, ylim=c(1,5), xlab="day")
 lines(y1s); lines(y1s+p1, lty=2); lines(y1s-p1, lty=2)
plot(PLT, type="p",  ylim=c(3,6), pch=19, xlab="day")
 lines(y2s); lines(y2s+p2, lty=2); lines(y2s-p2, lty=2)
plot(HCT, type="p", pch=19, ylim=c(20,40), xlab="day")
 lines(y3s); lines(y3s+p3, lty=2); lines(y3s-p3, lty=2)


Example 6.10

num = length(jj)
A = cbind(1,1,0,0)                                  

# Function to Calculate Likelihood 
Linn=function(para){
 Phi = diag(0,4); Phi[1,1] = para[1] 
 Phi[2,]=c(0,-1,-1,-1); Phi[3,]=c(0,1,0,0); Phi[4,]=c(0,0,1,0)
 cQ1 = para[2]; cQ2 = para[3]     # sqrt q11 and sqrt q22
 cQ=diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2
 cR = para[4]       # sqrt r11
 kf = Kfilter0(num,jj,A,mu0,Sigma0,Phi,cQ,cR)
 return(kf$like)  
 }

# Initial Parameters 
mu0 = c(.7,0,0,0);  Sigma0 = diag(.04,4)  
init.par = c(1.03,.1,.1,.5)   # Phi[1,1], the 2 Qs and R

# Estimation
est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1))
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimate=est$par,SE)
rownames(u)=c("Phi11","sigw1","sigw2","sigv"); u     

# Smooth
Phi = diag(0,4); Phi[1,1] = est$par[1] 
Phi[2,]=c(0,-1,-1,-1); Phi[3,]=c(0,1,0,0); Phi[4,]=c(0,0,1,0)
cQ1 = est$par[2]; cQ2 = est$par[3]      
cQ = diag(0,4); cQ[1,1]=cQ1; cQ[2,2]=cQ2   
cR = est$par[4]   
ks = Ksmooth0(num,jj,A,mu0,Sigma0,Phi,cQ,cR)   

# Plot 
Tsm = ts(ks$xs[1,,], start=1960, freq=4)
Ssm = ts(ks$xs[2,,], start=1960, freq=4)
p1 = 2*sqrt(ks$Ps[1,1,]); p2 = 2*sqrt(ks$Ps[2,2,])
par(mfrow=c(3,1))
plot(Tsm, main="Trend Component", ylab="Trend")
 lines(Tsm+p1, lty=2, col=4); lines(Tsm-p1,lty=2, col=4)
plot(Ssm, main="Seasonal Component", ylim=c(-5,4), ylab="Season")
 lines(Ssm+p2,lty=2, col=4); lines(Ssm-p2,lty=2, col=4)
plot(jj, type="p", main="Data (points) and Trend+Season (line)") 
 lines(Tsm+Ssm)

# Forecast 
n.ahead=12; y = ts(append(jj, rep(0,n.ahead)), start=1960, freq=4)
rmspe = rep(0,n.ahead); x00 = ks$xf[,,num]; P00 = ks$Pf[,,num]
Q=t(cQ)%*%cQ;  R=t(cR)%*%(cR)
for (m in 1:n.ahead){ 
 xp = Phi%*%x00; Pp = Phi%*%P00%*%t(Phi)+Q
 sig = A%*%Pp%*%t(A)+R; K = Pp%*%t(A)%*%(1/sig)
 x00 = xp; P00 = Pp-K%*%A%*%Pp
 y[num+m] = A%*%xp; rmspe[m] = sqrt(sig)  
}   

dev.new() 
plot(y, type="o", main="", ylab="", ylim=c(5,30), xlim=c(1975,1984))
 upp = ts(y[(num+1):(num+n.ahead)]+2*rmspe, start=1981, freq=4)
 low = ts(y[(num+1):(num+n.ahead)]-2*rmspe, start=1981, freq=4)
lines(upp, lty=2);  lines(low, lty=2);  abline(v=1980.75, lty=3) 


Example 6.12

dm = cmort-mean(cmort)   # center mortality
trend = time(cmort) - mean(time(cmort))

ded = ts.intersect(dm, u1=lag(tempr,-1), u2=part, u3=lag(part,-4), u4=trend, dframe=TRUE)

y = ded$dm; input = cbind(ded$u1, ded$u2, ded$u3, ded$u4)
num = length(y); A = array(c(1,0), dim = c(1,2,num))

# Function to Calculate Likelihood 
Linn=function(para){
  phi1=para[1]; phi2=para[2]; cR=para[3] 
  b1=para[4]; b2=para[5]; b3=para[6]; b4=para[7]
  mu0 = matrix(c(0,0),2,1); Sigma0=diag(100,2)
  Phi = matrix(c(phi1, phi2, 1, 0), 2)
  Theta = matrix(c(phi1, phi2), 2)
  Ups = matrix(c(b1,0,b2,0,b3,0,0,0),2,4)
  Gam = matrix(c(0,0,0,b4),1,4)
  cQ = cR;  S = cR^2
  kf = Kfilter2(num,y,A,mu0,Sigma0,Phi,Ups,Gam,Theta,cQ,cR,S,input)
  return(kf$like) 
}
# Estimation
phi1=.4;  phi2=.4;  cR=5; b1=-.1; b2=.1; b3=.1; b4=-1.5
init.par = c(phi1,phi2,cR,b1,b2,b3,b4)  # initial parameter values
U = c(1,1,10,1,1,1,5)  #  upper search bound for L-BFGS-B; lower is -U 
est = optim(init.par,Linn,NULL,method="L-BFGS-B",hessian=TRUE,
                     lower=-U,upper=U,control=list(trace=1,REPORT=1))
SE = sqrt(diag(solve(est$hessian)))

# Results    
u = cbind(estimate=est$par, SE)
rownames(u)=c("phi1","phi2","sigv","TL1","P","PL4","trnd"); u  


Example 6.13

################ 
# NOTE: Change lines below to tol=.01 or tol=.001 and nboot=100 or nboot=200 
#            if this takes a long time to run - depends on your machine
            
tol = sqrt(.Machine$double.eps)  # determines convergence of optimizer     
nboot = 500                      # number of bootstrap replicates     
################# 

if (!'plyr' %in% installed.packages()[,1]) install.packages(plyr)   # check if plyr installed
library(plyr)
 
y = window(qinfl, c(1953,1), c(1965,2))  # inflation   
z = window(qintr, c(1953,1), c(1965,2))  # interest   
num = length(y) 
A = array(z, dim=c(1,1,num))
input = matrix(1,num,1) 
    
# Function to Calculate Likelihood   
Linn=function(para){
  phi=para[1]; alpha=para[2]; b=para[3]; Ups=(1-phi)*b
  cQ=para[4]; cR=para[5]  
  kf=Kfilter2(num,y,A,mu0,Sigma0,phi,Ups,alpha,1,cQ,cR,0,input)
  return(kf$like)    
  }

# Parameter Estimation   
mu0=1; Sigma0=.01; phi=.84; alpha=-.77; b=.85; cQ=.12; cR=1.1
init.par = c(phi,alpha,b,cQ,cR)    # initial parameters   
est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1,reltol=tol))  
SE = sqrt(diag(solve(est$hessian)))                     
phi = est$par[1]; alpha=est$par[2]; b=est$par[3]; Ups=(1-phi)*b         
cQ=est$par[4]; cR=est$par[5] 
cbind(estimate=est$par, SE)  

# BEGIN BOOTSTRAP   
# Likelihood for the bootstrapped data   
Linn2=function(para){
  phi=para[1]; alpha=para[2]; b=para[3]; Ups=(1-phi)*b
  cQ=para[4]; cR=para[5]  
  kf=Kfilter2(num,y.star,A,mu0,Sigma0,phi,Ups,alpha,1,cQ,cR,0,input)
  return(kf$like)  }

# Run the filter at the estimates   
kf = Kfilter2(num,y,A,mu0,Sigma0,phi,Ups,alpha,1,cQ,cR,0,input)          

# Pull out necessary values from the filter and initialize  
xp=kf$xp; innov=kf$innov; sig=kf$sig; K=kf$K; e=innov/sqrt(sig)
e.star=e; y.star=y; xp.star=xp; k=4:50 
para.star = matrix(0, nboot, 5)  # to store estimates
init.par=c(.84,-.77,.85,.12,1.1)    

pr <- progress_text()            # displays progress
pr$init(nboot)

for (i in 1:nboot){
 # cat("iteration:", i, "\n")  # old counter
 pr$step()                     # new counter
 e.star[k] = sample(e[k], replace=TRUE)   
 for (j in k){
  xp.star[j] = phi*xp.star[j-1]+Ups+K[j]*sqrt(sig[j])*e.star[j] 
  }   
 y.star[k] = z[k]*xp.star[k]+alpha+sqrt(sig[k])*e.star[k]  
 est.star = optim(init.par, Linn2, NULL, method="BFGS", control=list(reltol=tol))     
 para.star[i,] = cbind(est.star$par[1], est.star$par[2], est.star$par[3], abs(est.star$par[4]), 
                      abs(est.star$par[5]))
} 
                      
# Some summary statistics  
rmse = rep(NA,5)  # SEs from the bootstrap
for(i in 1:5){rmse[i]=sqrt(sum((para.star[,i]-est$par[i])^2)/nboot)
              cat(i, rmse[i],"\n") 
}
              
# phi and sigw  
phi = para.star[,1]; sigw = abs(para.star[,4]) 
phi = ifelse(phi<0, NA, phi)  # any φ < 0 not plotted

panel.hist <- function(x, ...){    # scatterplot with histograms
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks; nB <- length(breaks)
  y <- h$counts; y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, ...)
}

u = cbind(phi, sigw)
colnames(u) = c("f","s")
pairs(u, cex=1.5, pch=1, diag.panel=panel.hist, cex.labels=1.5, font.labels=5)


Example 6.17

y = as.matrix(flu)
num = length(y)
nstate = 4;

M1 = as.matrix(cbind(1,0,0,1))  # obs matrix normal
M2 = as.matrix(cbind(1,0,1,1))  # obs matrix flu epi

prob = matrix(0,num,1); yp = y  # to store pi2(t|t-1) & y(t|t-1)
xfilter = array(0, dim=c(nstate,1,num)) # to store x(t|t)

# Function to Calculate Likelihood 
Linn = function(para){
  alpha1=para[1]; alpha2=para[2]; beta0=para[3]      
  sQ1=para[4];  sQ2=para[5];  like=0
  xf=matrix(0, nstate, 1)  # x filter
  xp=matrix(0, nstate, 1)  # x predict
  Pf=diag(.1, nstate)      # filter covar
  Pp=diag(.1, nstate)      # predict covar
  pi11 <- .75 -> pi22;  pi12 <- .25 -> pi21; pif1 <- .5 -> pif2            
  phi=matrix(0,nstate,nstate)
  phi[1,1]=alpha1; phi[1,2]=alpha2; phi[2,1]=1; phi[4,4]=1 
  Ups = as.matrix(rbind(0,0,beta0,0))
  Q = matrix(0,nstate,nstate)
  Q[1,1]=sQ1^2; Q[3,3]=sQ2^2; R=0  # R=0 in final model
  # begin filtering 
    for(i in 1:num){
    xp = phi%*%xf + Ups; Pp = phi%*%Pf%*%t(phi) + Q
    sig1 = as.numeric(M1%*%Pp%*%t(M1) + R)
    sig2 = as.numeric(M2%*%Pp%*%t(M2) + R)
    k1 = Pp%*%t(M1)/sig1; k2 = Pp%*%t(M2)/sig2 
    e1 = y[i]-M1%*%xp; e2 = y[i]-M2%*%xp
    pip1 = pif1*pi11 + pif2*pi21; pip2 = pif1*pi12 + pif2*pi22;
    den1 = (1/sqrt(sig1))*exp(-.5*e1^2/sig1); 
    den2 = (1/sqrt(sig2))*exp(-.5*e2^2/sig2);
    denom = pip1*den1 + pip2*den2;
    pif1 = pip1*den1/denom; pif2 = pip2*den2/denom;
    pif1=as.numeric(pif1); pif2=as.numeric(pif2)
    e1=as.numeric(e1); e2=as.numeric(e2)
    xf = xp + pif1*k1*e1 + pif2*k2*e2
    eye = diag(1, nstate)
    Pf = pif1*(eye-k1%*%M1)%*%Pp + pif2*(eye-k2%*%M2)%*%Pp 
    like = like - log(pip1*den1 + pip2*den2)
    prob[i]<<-pip2; xfilter[,,i]<<-xf; innov.sig<<-c(sig1,sig2)
    yp[i]<<-ifelse(pip1 > pip2, M1%*%xp, M2%*%xp)  
    }    
 return(like)   
 }
 
# Estimation
alpha1=1.4; alpha2=-.5; beta0=.3; sQ1=.1; sQ2=.1
init.par = c(alpha1, alpha2, beta0, sQ1, sQ2)

(est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1)))
SE = sqrt(diag(solve(est$hessian)))    
u = cbind(estimate=est$par, SE)
rownames(u)=c("alpha1","alpha2","beta0","sQ1","sQ2"); u

# Graphics 
predepi = ifelse(prob<.5,0,1); k = 6:length(y)       
Time = time(flu)[k]

par(mfrow=c(3,1), mar=c(2,3,1,1)+.1, cex=.9)     
plot(Time, y[k], type="o", ylim=c(0,1),ylab="")
 lines(Time, predepi[k],  lty="dashed", lwd=1.2)
 text(1979,.95,"(a)") 
plot(Time, xfilter[1,,k], type="l", ylim=c(-.1,.4), ylab="")
 lines(Time, xfilter[3,,k]); lines(Time, xfilter[4,,k])
 text(1979,.35,"(b)")
plot(Time, y[k], type="p", pch=1, ylim=c(.1,.9),ylab="")
 prde1 = 2*sqrt(innov.sig[1]); prde2 = 2*sqrt(innov.sig[2])
 prde = ifelse(predepi[k]<.5, prde1,prde2)
 lines(Time, yp[k]+prde, lty=2, lwd=1.5)
 lines(Time, yp[k]-prde, lty=2, lwd=1.5)
 text(1979,.85,"(c)")


Example 6.18

y = log(nyse^2); num=length(y)

# Initial Parameters
phi0=0; phi1=.95; sQ=.2; alpha=mean(y); sR0=1; mu1=-3; sR1=2
init.par = c(phi0,phi1,sQ,alpha,sR0,mu1,sR1)

# Innovations Likelihood 
Linn = function(para){
  phi0=para[1]; phi1=para[2]; sQ=para[3]; alpha=para[4]
  sR0=para[5]; mu1=para[6]; sR1=para[7]
  sv = SVfilter(num,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1)
  return(sv$like)    
}

# Estimation  
(est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1)))
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimates=est$par, SE)  
rownames(u)=c("phi0","phi1","sQ","alpha","sigv0","mu1","sigv1"); u

# Graphics   (need filters at the estimated parameters)
phi0=est$par[1]; phi1=est$par[2]; sQ=est$par[3]; alpha=est$par[4]
sR0=est$par[5]; mu1=est$par[6]; sR1=est$par[7]
sv = SVfilter(num,y,phi0,phi1,sQ,alpha,sR0,mu1,sR1) 

# densities plot (f is chi-sq, fm is fitted mixture)
x = seq(-15,6,by=.01) 
f = exp(-.5*(exp(x)-x))/(sqrt(2*pi))
f0 = exp(-.5*(x^2)/sR0^2)/(sR0*sqrt(2*pi))
f1 = exp(-.5*(x-mu1)^2/sR1^2)/(sR1*sqrt(2*pi))
fm = (f0+f1)/2
plot(x, f, type="l"); lines(x, fm, lty=2,lwd=2)

dev.new()
par(mfrow=c(2,1))
Time = 801:1000
plot(Time, y[Time], type="l", main="log(Squared NYSE Returns)")
plot(Time, sv$xp[Time],type="l", main="Predicted log-Volatility", ylim=c(-1.5,1.8), ylab="", xlab="")
 lines(Time, sv$xp[Time]+2*sqrt(sv$Pp[Time]), lty="dashed")
 lines(Time, sv$xp[Time]-2*sqrt(sv$Pp[Time]), lty="dashed")


Example 6.19

n.boot = 500        # number of bootstrap replicates
tol = sqrt(.Machine$double.eps)  # convergence limit 

gnpgr = diff(log(gnp))
fit = arima(gnpgr, order=c(1,0,0))
y = as.matrix(log(resid(fit)^2))
num = length(y) 
plot.ts(y, ylab="")

# Initial Parameters
phi1 = .9; sQ = .5; alpha = mean(y); sR0 = 1; mu1 = -3; sR1 = 2.5
init.par = c(phi1, sQ, alpha, sR0, mu1, sR1)

# Innovations Likelihood 
Linn=function(para){
  phi1 = para[1]; sQ = para[2]; alpha = para[3]
  sR0 = para[4]; mu1 = para[5]; sR1 = para[6]
  sv = SVfilter(num, y, 0, phi1, sQ, alpha, sR0, mu1, sR1)
  return(sv$like)    
  }

# Estimation  
(est = optim(init.par, Linn, NULL, method="BFGS", hessian=TRUE, control=list(trace=1,REPORT=1)))
SE = sqrt(diag(solve(est$hessian)))
u = cbind(estimates=est$par, SE)  
rownames(u)=c("phi1","sQ","alpha","sig0","mu1","sig1"); u

# Bootstrap 
para.star = matrix(0, n.boot, 6)  # to store parameter estimates
Linn2 = function(para){   
  phi1 = para[1]; sQ = para[2]; alpha = para[3]
  sR0 = para[4]; mu1 = para[5]; sR1 = para[6]
  sv = SVfilter(num, y.star, 0, phi1, sQ, alpha, sR0, mu1, sR1)
  return(sv$like)    
  }

for (jb in 1:n.boot){  cat("iteration:", jb, "\n")
 phi1 = est$par[1]; sQ = est$par[2]; alpha = est$par[3]    
 sR0 = est$par[4]; mu1 = est$par[5]; sR1 = est$par[6] 
 Q = sQ^2; R0 = sR0^2; R1 = sR1^2
 sv = SVfilter(num, y, 0, phi1, sQ, alpha, sR0, mu1, sR1)
 
 sig0 = sv$Pp+R0 
 sig1 = sv$Pp+R1 
 K0   = sv$Pp/sig0 
 K1   = sv$Pp/sig1
 inn0 = y-sv$xp-alpha; inn1 = y-sv$xp-mu1-alpha
 den1 = (1/sqrt(sig1))*exp(-.5*inn1^2/sig1)
 den0 = (1/sqrt(sig0))*exp(-.5*inn0^2/sig0)
 fpi1 = den1/(den0+den1)
 
 # (start resampling at t=4)
 e0   = inn0/sqrt(sig0)
 e1   = inn1/sqrt(sig1)  
 indx = sample(4:num, replace=TRUE)
 sinn = cbind(c(e0[1:3], e0[indx]), c(e1[1:3], e1[indx]))
 eF   = matrix(c(phi1, 1, 0, 0), 2, 2)
 xi   = cbind(sv$xp,y)    # initialize
   
   for (i in 4:num){    # generate boot sample
   G   = matrix(c(0, alpha+fpi1[i]*mu1), 2, 1)
   h21 = (1-fpi1[i])*sqrt(sig0[i]); h11 = h21*K0[i]
   h22 = fpi1[i]*sqrt(sig1[i]); h12 = h22*K1[i]
   H   = matrix(c(h11,h21,h12,h22),2,2)
   xi[i,] = t(eF%*%as.matrix(xi[i-1,],2) + G + H%*%as.matrix(sinn[i,],2))
   }
   
 # Estimates from boot data 
 y.star = xi[,2]
 phi1 = .9; sQ = .5; alpha = mean(y.star); sR0 = 1; mu1 = -3; sR1 = 2.5
 init.par = c(phi1, sQ, alpha, sR0, mu1, sR1)   #  same as for data
 est.star = optim(init.par, Linn2, NULL, method="BFGS", control=list(reltol=tol)) 
 para.star[jb,] = cbind(est.star$par[1], abs(est.star$par[2]), est.star$par[3], abs(est.star$par[4]), 
                        est.star$par[5], abs(est.star$par[6])) 
}

# Some summary statistics and graphics  
rmse = rep(NA,6)  # SEs from the bootstrap
for(i in 1:6){rmse[i] = sqrt(sum((para.star[,i]-est$par[i])^2)/n.boot)
                cat(i, rmse[i],"\n") 
             }  
dev.new()
phi = para.star[,1]
hist(phi, 15, prob=TRUE, main="", xlim=c(.4,1.2), xlab="")
u = seq(.4, 1.2, by=.01)
lines(u,dnorm(u, mean=est$par[1], sd=SE[1]), lty="dashed", lwd=2)


Example 6.23

#########################################################################################
# This example adapted from code by: Hedibert Freitas Lopes
#########################################################################################

if (!'plyr' %in% installed.packages()[,1]) install.packages(plyr)   # check if plyr installed
library(plyr)
 
# setup
mc.seed    = 90210
burn.in    = 1000   
step.size  = 3
n.mcmc     = 2000*step.size + burn.in  

g = function(arg){c(arg[1],arg[1]/(1+arg[1]^2),cos(1.2*arg[2]))}

fullxt = function(t,yt,xtm1,xt,xtp1,theta,sig,tau){
  zm1 = g(c(xtm1,t-1))
  z   = g(c(xt,t))
  full = dnorm(xt,sum(zm1*theta),tau,log=TRUE)+
         dnorm(xtp1,sum(z*theta),tau,log=TRUE)+
         dnorm(yt,xt^2/20,sig,log=TRUE)
  return(full)
}

fullxn = function(t,yt,xtm1,xt,theta,sig,tau){
  zm1  = g(c(xtm1,t-1))
  z    = g(c(xt,t))
  full = dnorm(xt,sum(zm1*theta),tau,log=TRUE)+
         dnorm(yt,xt^2/20,sig,log=TRUE)
  return(full)
}

pr <- progress_text()           # displays progress
pr$init(n.mcmc)

draw.x = function(y,x,x0,theta,sig,tau,v){
  x1  = rnorm(1,x[1],v)
  num = fullxt(1,y[1],x0,x1  ,x[2],theta,sig,tau)
  den = fullxt(1,y[1],x0,x[1],x[2],theta,sig,tau)
  if (log(runif(1))<(num-den)){x[1]=x1}
  xn  = rnorm(1,x[n],v)
  num = fullxn(n,y[n],x[n-1],xn  ,theta,sig,tau)
  den = fullxn(n,y[n],x[n-1],x[n],theta,sig,tau)
  if (log(runif(1))<(num-den)){x[n]=xn} 
  for (t in 2:(n-1)){
    xt  = rnorm(1,x[t],v)
    num = fullxt(t,y[t],x[t-1],xt  ,x[t+1],theta,sig,tau)
    den = fullxt(t,y[t],x[t-1],x[t],x[t+1],theta,sig,tau)
    if (log(runif(1))<(num-den)){x[t]=xt}
  }
  pr$step()
  return(x)
}

fixedpar = function(y,X,b,A,v,lam){
  n     = length(y)
  k     = ncol(X)
  par1  = (v+n)/2
  var   = solve(crossprod(X,X)+A)
  mean  = matrix(var%*%(crossprod(X,y)+crossprod(A,b)),k,1)
  par2  = v*lam + sum((y-crossprod(t(X),mean))^2)
  par2  = (par2 + crossprod(t(crossprod(mean-b,A)),mean-b))/2
  sig2  = 1/rgamma(1,par1,par2)
  var   = var*sig2
  mean  = mean + crossprod(t(chol(var)),rnorm(k))
  return(c(mean,sig2))
}

quant005   = function(x){quantile(x,0.005)}
quant995   = function(x){quantile(x,0.995)}

# Simulating the data
set.seed(mc.seed)
n         = 100
sig2      = 1
tau2      = 10
sig       = sqrt(sig2)
tau       = sqrt(tau2)
theta     = c(0.5,25,8)
ptr       = c(theta,tau2,sig2)
y         = rep(0,n)
x         = rep(0,n)
x0        = 0
Z         = g(c(x0,0))
x[1]      = rnorm(1,sum(Z*theta),tau)
y[1]      = rnorm(1,x[1]^2/20,sig)
for (t in 2:n){
  Z    = g(c(x[t-1],t-1))
  x[t] = rnorm(1,sum(Z*theta),tau)
  y[t] = rnorm(1,x[t]^2/20,sig)
}
xtr = x

# process graphics
par(mfrow=c(2,2), mar=c(2.5,2.5,1,1), mgp=c(1.6,.6,0))
ts.plot(y, xlab="time", ylab=expression(y[~t]))
ts.plot(x, xlab="time", ylab=expression(x[~t]))
plot(x, y, xlab=expression(x[~t]), ylab=expression(y[~t]))
plot.ts(x[-n], x[-1], ylab=expression(x[~t]), xlab=expression(x[~t-1]), cex=.8)


# Prior hyperparameters
m0        = 0.0
C0        = 10.0
n0        = 6
sig20     = 2/3
theta0    = c(0.5,25,8)
V0        = diag(c(0.25,10,4))/tau2
nu0       = 6
tau20     = 20/3
iV0       = diag(1/diag(V0))
n0sig202  = n0*sig20/2
            c((n0*sig20/2)/(n0/2-1),((n0*sig20/2)/(n0/2-1))^2/(n0/2-2),
            (nu0*tau20/2)/(nu0/2-1),((nu0*tau20/2)/(nu0/2-1))^2/(nu0/2-2))

# MCMC setup
v     = 0.5
niter = n.mcmc   
theta = ptr[1:3]
tau2  = ptr[4]
sig2  = ptr[5]
x     = xtr
xs    = NULL
ps    = NULL
n0n2  = (n0+n)/2
for (iter in 1:niter){	
  x     = draw.x(y,x,x0,theta,sig,tau,v)
  sig2  = 1/rgamma(1,n0n2,n0sig202+sum((y-x^2/20)^2)/2)
  sig   = sqrt(sig2)
  W     = cbind(c(x0,x[1:(n-1)]),0:(n-1))
  Z     = t(apply(W,1,g))
  par   = fixedpar(x,Z,theta0,iV0,nu0,tau20)
  theta = par[1:3]
  tau2  = par[4]
  tau   = sqrt(tau2)
  xs    = rbind(xs,x)
  ps    = rbind(ps,c(par,sig2))
}

M0 = burn.in     
M  = n.mcmc-M0  
L  = step.size   

index = seq((M0+1),niter,by=L) 
length(index)

mx = apply(xs[index,],2,median)
lx = apply(xs[index,],2,quant005)
ux = apply(xs[index,],2,quant995)

# display state estimation
dev.new()
par(mar=c(3,3,1.5,1), mgp=c(1.6,.6,0))  
plot(xtr,type="p",pch=20,ylim=range(lx,ux+10,xtr),xlab="time",ylab="",main="State Estimation")
lines(mx,col=1,type="o")
legend("topleft",legend=c("True","Posterior Median","99%CI"),col=c(1,1,gray(.5, alpha=.4)),
       lty=c(0,1,1), pch=c(20,1,-1), lwd=c(1,1,5), cex=.8)
 xx=c(time(xtr), rev(time(xtr)))
 yy=c(lx, rev(ux))
polygon(xx, yy, border=NA, col=gray(.5, alpha = .4))


# display parameter estimation
names = c(expression(alpha), expression(beta), expression(gamma),
          expression(sigma[w]^2), expression(sigma[v]^2))

dev.new()
par(mfrow=c(3,5), mar=c(2.75,2.75,2,1), mgp=c(1.6,.6,0))
for (i in 1:5){
  plot(1:length(index),ps[index,i],type="l",xlab="iteration",ylab="")
  title(names[i], cex=1.5)
  abline(h=ptr[i],col=2,lwd=4)
}
for (i in 1:5)
  acf(ps[index,i],main="")
for (i in 1:5){
  hist(ps[index,i],prob=TRUE, xlab="",ylab="",main="")
  lines(density(ps[index,i],adjust=2))
   lp = apply(as.matrix(ps[index,i]),2,quant005)
   up = apply(as.matrix(ps[index,i]),2,quant995)
  abline(v=c(ptr[i],lp,up),col=c(2,4,4),lwd=c(4,1,1), lty=c(1,2,2))
}

[+] Chapter 7

Code in Introduction

x = matrix(0, 128, 6)
for (i in 1:6) x[,i] = rowMeans(fmri[[i]])
colnames(x)=c("Brush", "Heat", "Shock", "Brush", "Heat", "Shock")
plot.ts(x, main="")
mtext("Awake", side=3, line=1.2, adj=.05,cex=1.2)
mtext("Sedated", side=3, line=1.2, adj=.85, cex=1.2) 

attach(eqexp)
P = 1:1024; S = P+1024
x = cbind(EQ5[P], EQ6[P], EX5[P], EX6[P], NZ[P], EQ5[S], EQ6[S], EX5[S], EX6[S], NZ[S])
x.name = c("EQ5","EQ6","EX5","EX6","NZ")
colnames(x) = c(x.name, x.name)
plot.ts(x, main="")
mtext("P waves", side=3, line=1.2, adj=.05, cex=1.2)
mtext("S waves", side=3, line=1.2, adj=.85, cex=1.2)


Example 7.1

attach(climhyd)
plot.ts(climhyd)    # figure 7.3
Y = climhyd         # Y holds the transformed series
Y[,6] = log(Y[,6])  # log inflow
Y[,5] = sqrt(Y[,5]) # sqrt precipitation 
                              
L = 25              # setup 
M = 100 
alpha = .001
fdr = .001
nq = 2              # number of inputs  (Temp and Precip)

# Spectral Matrix 
Yspec = mvspec(Y, spans=L, kernel="daniell", taper=.1, plot=FALSE)
n = Yspec$n.used          # effective sample size
Fr = Yspec$freq           # fundamental freqs 
n.freq = length(Fr)       # number of frequencies
Yspec$bandwidth           # = 0.05 
 
# Coherencies (see section 4.7 also)
Fq = qf(1-alpha, 2, L-2); cn = Fq/(L-1+Fq)
plt.name = c("(a)","(b)","(c)","(d)","(e)","(f)")
dev.new()
par(mfrow=c(2,3), cex.lab=1.2) 

# The coherencies are listed as 1,2,...,15=choose(6,2) 
for (i in 11:15){
 plot(Fr,Yspec$coh[,i], type="l", ylab="Sq Coherence", xlab="Frequency", ylim=c(0,1), 
      main=c("Inflow with", names(climhyd[i-10])))
 abline(h = cn); text(.45,.98, plt.name[i-10], cex=1.2)  
} 

 # Multiple Coherency 
coh.15 = stoch.reg(Y, cols.full = c(1,5), cols.red = NULL, alpha, L, M, plot.which = "coh")  
text(.45,.98, plt.name[6], cex=1.2) 
title(main = c("Inflow with", "Temp and Precip")) 

# Partial F (note F-stat is called eF in the code)
numer.df = 2*nq
denom.df = Yspec$df-2*nq

dev.new()
par(mfrow=c(3,1), mar=c(3,3,2,1)+.5, mgp = c(1.5,0.4,0), cex.lab=1.2)  
out.15 = stoch.reg(Y, cols.full = c(1,5), cols.red = 5, alpha, L, M, plot.which = "F.stat")
 eF = out.15$eF 
 pvals = pf(eF, numer.df, denom.df, lower.tail = FALSE)
 pID = FDR(pvals, fdr)
abline(h=c(eF[pID]), lty=2)
title(main = "Partial F Statistic")

# Regression Coefficients
S = seq(from = -M/2+1, to = M/2 - 1, length = M-1)

plot(S, coh.15$Betahat[,1], type = "h", xlab = "", ylab =names(climhyd[1]), 
        ylim = c(-.025, .055), lwd=2)
abline(h=0)
title(main = "Impulse Response Functions")

plot(S, coh.15$Betahat[,2], type = "h", xlab = "Index", ylab = names(climhyd[5]), 
        ylim = c(-.015, .055), lwd=2)
abline(h=0)


Example 7.2

attach(beamd)
tau = rep(0,3) 
 u = ccf(sensor1, sensor2, plot=FALSE)
tau[1] = u$lag[which.max(u$acf)]    #  17
 u = ccf(sensor3, sensor2, plot=FALSE)
tau[3] = u$lag[which.max(u$acf)]    # -22

Y = ts.union(lag(sensor1,tau[1]), lag(sensor2, tau[2]), lag(sensor3, tau[3]))
beam = rowMeans(Y)

par(mfrow=c(4,1), mar=c(0,5.1,0,5.1), oma=c(6,0,5,0))
plot.ts(sensor1, xaxt="no")
title(main="Infrasonic Signals and Beam", outer=TRUE)
plot.ts(sensor2, xaxt="no");  plot.ts(sensor3, xaxt="no")
plot.ts(beam);  title(xlab="Time", outer=TRUE)


Example 7.4

attach(beamd)
L = 9
fdr = .001 
N = 3
Y = cbind(beamd, beam=rowMeans(beamd))
n = nextn(nrow(Y))

Y.fft = mvfft(as.ts(Y))/sqrt(n)
Df = Y.fft[,1:3]   # fft of the data
Bf = Y.fft[,4]     # beam fft 

ssr = N*Re(Bf*Conj(Bf))               # raw signal spectrum
sse = Re(rowSums(Df*Conj(Df))) - ssr  # raw error spectrum

# Smooth
SSE = filter(sse, sides=2, filter=rep(1/L,L), circular=TRUE)
SSR = filter(ssr, sides=2, filter=rep(1/L,L), circular=TRUE)
SST = SSE + SSR

par(mfrow=c(2,1), mar=c(4,4,2,1)+.1)
Fr  = 0:(n-1)/n
nFr = 1:200   # freqs to plot

plot( Fr[nFr], SST[nFr], type="l", ylab="log Power", xlab="", main="Sum of Squares",log="y")
lines(Fr[nFr], SSE[nFr], type="l", lty=2)

eF = (N-1)*SSR/SSE; df1 = 2*L;  df2 = 2*L*(N-1)
pvals = pf(eF, df1, df2, lower=FALSE)  # p values for FDR
pID = FDR(pvals, fdr); Fq = qf(1-fdr, df1, df2)  

plot(Fr[nFr], eF[nFr], type="l", ylab="F-statistic", xlab="Frequency",  main="F Statistic")
abline(h=c(Fq, eF[pID]), lty=1:2)


Example 7.6

attach(beamd)
L  = 9
M  = 100 
M2 = M/2
N  = 3   
Y  = cbind(beamd, beam <- rowMeans(beamd))
n  = nextn(nrow(Y))
n.freq = n/2
 
Y[,1:3] = Y[,1:3]-Y[,4]  # center each series  

Y.fft = mvfft(as.ts(Y))/sqrt(n)
Ef    = Y.fft[,1:3]              # fft of the error
Bf    = Y.fft[,4]                # beam fft
ssr   = N*Re(Bf*Conj(Bf))        # Raw Signal Spectrum 
sse   = Re(rowSums(Ef*Conj(Ef))) # Raw Error Spectrum

# Smooth
SSE = filter(sse, sides=2, filter=rep(1/L,L), circular=TRUE)
SSR = filter(ssr, sides=2, filter=rep(1/L,L), circular=TRUE)

# Estimate Signal and Noise Spectra
fv = SSE/(L*(N-1))          # Equation (7.77)
fb = (SSR-SSE/(N-1))/(L*N)  # Equation (7.78)
fb[fb<0] = 0 

H0 = N*fb/(fv+N*fb)     
H0[ceiling(.04*n):n] = 0    # zero out H0 beyond frequency .04

# Extend components to make it a valid transform
H0 = c(H0[1:n.freq], rev(H0[2:(n.freq+1)]))  
h0 = Re(fft(H0, inverse = TRUE))            # Impulse Response 
h0 = c(rev(h0[2:(M2+1)]), h0[1:(M2+1)])     # center it
h1 = spec.taper(h0, p = .5)                 # taper it
k1 = h1/sum(h1)                             # normalize it
f.beam = filter(Y$beam, filter=k1, sides=2) # filter it

# Graphics  
nFr = 1:50      # freqs to display
Fr = (nFr-1)/n  # frequencies 

layout(matrix(c(1, 2, 4, 1, 3, 4), nc=2)) 
par(mar=c(4,4,2,1) + .1)
plot(10*Fr, fb[nFr], type="l", ylab="Power", xlab="Frequency (Hz)") 
 lines(10*Fr, fv[nFr], lty=2)   
 text(.24, 5, "(a)", cex=1.2)
plot(10*Fr,H0[nFr],type="l",ylab="Frequency Response",xlab="Frequency (Hz)") 
 text(.23, .84, "(b)", cex=1.2)
plot(-M2:M2,k1,type="l",ylab="Impulse Response",xlab="Index", lwd=1.5)
 text(45, .022, "(c)", cex=1.2)
ts.plot(cbind(f.beam,beam), lty=1:2, ylab="beam")  
 text(2040, 2, "(d)", cex=1.2)


Example 7.7

n          = 128               # length of series
n.freq     = 1 + n/2           # number of frequencies
Fr         = (0:(n.freq-1))/n  # the frequencies 
N          = c(5,4,5,3,5,4)    # number of series for each cell
n.subject  = sum(N)            # number of subjects (26)
n.trt      = 6                 # number of treatments
L          = 3                 # for smoothing
num.df     = 2*L*(n.trt-1)     # dfs for F test
den.df     = 2*L*(n.subject-n.trt)


# Design Matrix (Z): 
Z1 = outer(rep(1,N[1]), c(1,1,0,0,0,0))
Z2 = outer(rep(1,N[2]), c(1,0,1,0,0,0))
Z3 = outer(rep(1,N[3]), c(1,0,0,1,0,0))
Z4 = outer(rep(1,N[4]), c(1,0,0,0,1,0)) 
Z5 = outer(rep(1,N[5]), c(1,0,0,0,0,1)) 
Z6 = outer(rep(1,N[6]), c(1,-1,-1,-1,-1,-1)) 

Z  = rbind(Z1, Z2, Z3, Z4, Z5, Z6)
ZZ = t(Z)%*%Z 

SSEF <- rep(NA, n) -> SSER   

HatF = Z%*%solve(ZZ, t(Z))
HatR = Z[,1]%*%t(Z[,1])/ZZ[1,1]

par(mfrow=c(3,3), mar=c(3.5,4,0,0), oma=c(0,0,2,2), mgp = c(1.6,.6,0))
loc.name = c("Cortex 1","Cortex 2","Cortex 3","Cortex 4","Caudate","Thalamus 1","Thalamus 2",
             "Cerebellum 1","Cerebellum 2")

for(Loc in 1:9) {   
 i = n.trt*(Loc-1)   
 Y = cbind(fmri[[i+1]], fmri[[i+2]], fmri[[i+3]], fmri[[i+4]], fmri[[i+5]], fmri[[i+6]])
 Y = mvfft(spec.taper(Y, p=.5))/sqrt(n)	
 Y = t(Y)      # Y is now 26 x 128 FFTs

 # Calculation of Error Spectra 
 for (k in 1:n) {   
  SSY = Re(Conj(t(Y[,k]))%*%Y[,k])
  SSReg = Re(Conj(t(Y[,k]))%*%HatF%*%Y[,k])
 SSEF[k] = SSY - SSReg
  SSReg = Re(Conj(t(Y[,k]))%*%HatR%*%Y[,k])
 SSER[k] = SSY - SSReg  
 }

# Smooth 
sSSEF = filter(SSEF, rep(1/L, L), circular = TRUE)
sSSER = filter(SSER, rep(1/L, L), circular = TRUE)

eF =(den.df/num.df)*(sSSER-sSSEF)/sSSEF

plot(Fr, eF[1:n.freq], type="l", xlab="Frequency", ylab="F Statistic", ylim=c(0,7))
abline(h=qf(.999, num.df, den.df),lty=2)
text(.25, 6.5, loc.name[Loc], cex=1.2)   
}


Example 7.8

n      = 128
n.freq = 1 + n/2  
Fr     = (0:(n.freq-1))/n  
nFr    = 1:(n.freq/2) 
N      = c(5,4,5,3,5,4)  
n.para = 6          # number of parameters
n.subject = sum(N)  # total number of subjects

L = 3   
df.stm = 2*L*(3-1)              # stimulus (3 levels: Brush,Heat,Shock)
df.con = 2*L*(2-1)              # conscious (2 levels: Awake,Sedated)  
df.int = 2*L*(3-1)*(2-1)        # interaction 
den.df = 2*L*(n.subject-n.para) # df for full model 

# Design Matrix:          mu  a1  a2   b  g1  g2   
 Z1 = outer(rep(1,N[1]), c(1,  1,  0,  1,  1,  0)) 
 Z2 = outer(rep(1,N[2]), c(1,  0,  1,  1,  0,  1)) 
 Z3 = outer(rep(1,N[3]), c(1, -1, -1,  1, -1, -1)) 
 Z4 = outer(rep(1,N[4]), c(1,  1,  0, -1, -1,  0)) 
 Z5 = outer(rep(1,N[5]), c(1,  0,  1, -1,  0, -1)) 
 Z6 = outer(rep(1,N[6]), c(1, -1, -1, -1,  1,  1))
 
 Z = rbind(Z1, Z2, Z3, Z4, Z5, Z6)  
 ZZ = t(Z)%*%Z

rep(NA, n)-> SSEF -> SSE.stm -> SSE.con -> SSE.int              
HatF    = Z%*%solve(ZZ,t(Z))     
Hat.stm = Z[,-(2:3)]%*%solve(ZZ[-(2:3),-(2:3)], t(Z[,-(2:3)]))  
Hat.con = Z[,-4]%*%solve(ZZ[-4,-4], t(Z[,-4]))   
Hat.int = Z[,-(5:6)]%*%solve(ZZ[-(5:6),-(5:6)], t(Z[,-(5:6)]))                                                            

par(mfrow=c(5,3),mar=c(3.5,4,0,0),oma=c(0,0,2,2),mgp = c(1.6,.6,0)) 
loc.name = c("Cortex 1","Cortex 2","Cortex 3","Cortex 4","Caudate", "Thalamus 1","Thalamus 2",
             "Cerebellum 1","Cerebellum 2")
for(Loc in c(1:4,9)) {   # only Loc 1 to 4 and 9 used                
  i = 6*(Loc-1)                                                                 
  Y = cbind(fmri[[i+1]], fmri[[i+2]], fmri[[i+3]], fmri[[i+4]], fmri[[i+5]], fmri[[i+6]])                     
  Y = mvfft(spec.taper(Y, p=.5))/sqrt(n)  
  Y = t(Y)  
  for (k in 1:n) {    
    SSY=Re(Conj(t(Y[,k]))%*%Y[,k])
    SSReg= Re(Conj(t(Y[,k]))%*%HatF%*%Y[,k])
  SSEF[k]=SSY-SSReg 
    SSReg=Re(Conj(t(Y[,k]))%*%Hat.stm%*%Y[,k])
  SSE.stm[k] = SSY-SSReg
    SSReg=Re(Conj(t(Y[,k]))%*%Hat.con%*%Y[,k])
  SSE.con[k]=SSY-SSReg  
    SSReg=Re(Conj(t(Y[,k]))%*%Hat.int%*%Y[,k])
  SSE.int[k]=SSY-SSReg   
  }
 # Smooth    
  sSSEF    = filter(SSEF, rep(1/L, L), circular = TRUE)
  sSSE.stm = filter(SSE.stm, rep(1/L, L), circular = TRUE)
  sSSE.con = filter(SSE.con, rep(1/L, L), circular = TRUE)
  sSSE.int = filter(SSE.int, rep(1/L, L), circular = TRUE)  
  eF.stm   = (den.df/df.stm)*(sSSE.stm-sSSEF)/sSSEF
  eF.con   = (den.df/df.con)*(sSSE.con-sSSEF)/sSSEF
  eF.int   = (den.df/df.int)*(sSSE.int-sSSEF)/sSSEF

 plot(Fr[nFr],eF.stm[nFr], type="l", xlab="Frequency", ylab="F Statistic", ylim=c(0,12))
   abline(h=qf(.999, df.stm, den.df),lty=2)       
  if(Loc==1) mtext("Stimulus", side=3, line=.3, cex=1)
  mtext(loc.name[Loc], side=2, line=3, cex=.9)
 plot(Fr[nFr],eF.con[nFr], type="l", xlab="Frequency", ylab="F Statistic", ylim=c(0,12))
  abline(h=qf(.999, df.con, den.df),lty=2)
  if(Loc==1)  mtext("Consciousness", side=3, line=.3, cex=1)   
 plot(Fr[nFr],eF.int[nFr], type="l", xlab="Frequency",ylab="F Statistic", ylim=c(0,12))
  abline(h=qf(.999, df.int, den.df),lty=2)
  if(Loc==1) mtext("Interaction", side=3, line= .3, cex=1)   
} 


Example 7.9

n      = 128
n.freq = 1 + n/2  
Fr     = (0:(n.freq-1))/n
nFr    = 1:(n.freq/2) 
N      = c(5,4,5,3,5,4)
L      = 3
n.subject = sum(N)
   
# Design Matrix            
Z1 = outer(rep(1,N[1]), c(1,0,0,0,0,0)) 
Z2 = outer(rep(1,N[2]), c(0,1,0,0,0,0)) 
Z3 = outer(rep(1,N[3]), c(0,0,1,0,0,0)) 
Z4 = outer(rep(1,N[4]), c(0,0,0,1,0,0)) 
Z5 = outer(rep(1,N[5]), c(0,0,0,0,1,0)) 
Z6 = outer(rep(1,N[6]), c(0,0,0,0,0,1))
Z  = rbind(Z1, Z2, Z3, Z4, Z5, Z6)
ZZ = t(Z)%*%Z 

A      = rbind(diag(1,3), diag(1,3))   # Contrasts:  6 x 3 
nq     = nrow(A)
num.df = 2*L*nq
den.df = 2*L*(n.subject-nq)                   
HatF   = Z%*%solve(ZZ, t(Z))           # full model hat matrix   

rep(NA, n)-> SSEF -> SSER 
eF = matrix(0,n,3)

par(mfrow=c(5,3), mar=c(3.5,4,0,0), oma=c(0,0,2,2), mgp = c(1.6,.6,0)) 

loc.name = c("Cortex 1","Cortex 2","Cortex 3","Cortex 4","Caudate","Thalamus 1","Thalamus 2",
             "Cerebellum 1","Cerebellum 2")
cond.name = c("Brush", "Heat", "Shock")             

for(Loc in c(1:4,9)) {
 i = 6*(Loc-1)                                                                 
 Y = cbind(fmri[[i+1]],fmri[[i+2]],fmri[[i+3]],fmri[[i+4]], fmri[[i+5]],fmri[[i+6]])                       
 Y = mvfft(spec.taper(Y, p=.5))/sqrt(n); Y = t(Y)  
 for (cond in 1:3){
  Q = t(A[,cond])%*%solve(ZZ, A[,cond])
  HR = A[,cond]%*%solve(ZZ, t(Z))      
  for (k in 1:n){   
    SSY = Re(Conj(t(Y[,k]))%*%Y[,k])
    SSReg= Re(Conj(t(Y[,k]))%*%HatF%*%Y[,k])
   SSEF[k]= (SSY-SSReg)*Q 
    SSReg= HR%*%Y[,k]
   SSER[k] = Re(SSReg*Conj(SSReg))  
  }     

# Smooth      
 sSSEF = filter(SSEF, rep(1/L, L), circular = TRUE)
 sSSER = filter(SSER, rep(1/L, L), circular = TRUE)       
 eF[,cond]= (den.df/num.df)*(sSSER/sSSEF)   }                                    
 plot(Fr[nFr], eF[nFr,1], type="l", xlab="Frequency", ylab="F Statistic", ylim=c(0,5))
  abline(h=qf(.999, num.df, den.df),lty=2)       
  if(Loc==1) mtext("Brush", side=3, line=.3, cex=1)
  mtext(loc.name[Loc], side=2, line=3, cex=.9)
 plot(Fr[nFr], eF[nFr,2], type="l", xlab="Frequency", ylab="F Statistic", ylim=c(0,5))
  abline(h=qf(.999, num.df, den.df),lty=2)
  if(Loc==1)  mtext("Heat", side=3, line=.3, cex=1)   
 plot(Fr[nFr], eF[nFr,3], type="l", xlab="Frequency", ylab="F Statistic", ylim=c(0,5))
  abline(h = qf(.999, num.df, den.df) ,lty=2)
  if(Loc==1) mtext("Shock", side=3, line=.3, cex=1)  
}  


Example 7.10

P = 1:1024 
S = P+1024
N = 8
n = 1024
p.dim = 2
m = 10
L = 2*m+1

eq.P   = as.ts(eqexp[P,1:8])
eq.S   = as.ts(eqexp[S,1:8]) 
eq.m   = cbind(rowMeans(eq.P), rowMeans(eq.S))
ex.P   = as.ts(eqexp[P,9:16])
ex.S   = as.ts(eqexp[S,9:16]) 
ex.m   = cbind(rowMeans(ex.P), rowMeans(ex.S)) 
m.diff = mvfft(eq.m - ex.m)/sqrt(n)

eq.Pf = mvfft(eq.P-eq.m[,1])/sqrt(n)
eq.Sf = mvfft(eq.S-eq.m[,2])/sqrt(n)
ex.Pf = mvfft(ex.P-ex.m[,1])/sqrt(n)
ex.Sf = mvfft(ex.S-ex.m[,2])/sqrt(n) 

fv11 = rowSums(eq.Pf*Conj(eq.Pf)) + rowSums(ex.Pf*Conj(ex.Pf))/(2*(N-1))
fv12 = rowSums(eq.Pf*Conj(eq.Sf)) + rowSums(ex.Pf*Conj(ex.Sf))/(2*(N-1))
fv22 = rowSums(eq.Sf*Conj(eq.Sf)) + rowSums(ex.Sf*Conj(ex.Sf))/(2*(N-1))
fv21 = Conj(fv12)

# Equal Means  
T2 = rep(NA, 512)
for (k  in 1:512){
 fvk = matrix(c(fv11[k], fv21[k], fv12[k], fv22[k]), 2, 2)
 dk = as.matrix(m.diff[k,])
 T2[k] = Re((N/2)*Conj(t(dk))%*%solve(fvk,dk))  }
eF = T2*(2*p.dim*(N-1))/(2*N-p.dim-1)
par(mfrow=c(2,2), mar=c(3,3,2,1), mgp = c(1.6,.6,0), cex.main=1.1)
freq = 40*(0:511)/n  # in Hz (cycles per second)  
plot(freq, eF, type="l", xlab="Frequency (Hz)", ylab="F Statistic", main="Equal Means")
abline(h=qf(.999, 2*p.dim, 2*(2*N-p.dim-1)))

# Equal P   
kd    = kernel("daniell",m);  
u     = Re(rowSums(eq.Pf*Conj(eq.Pf))/(N-1))
feq.P = kernapply(u, kd, circular=TRUE)
u     = Re(rowSums(ex.Pf*Conj(ex.Pf))/(N-1))
fex.P =	kernapply(u, kd, circular=TRUE)

plot(freq, feq.P[1:512]/fex.P[1:512], type="l", xlab="Frequency (Hz)", ylab="F Statistic", 
      main="Equal P-Spectra")
abline(h = qf(.999, 2*L*(N-1),  2*L*(N-1))) 

# Equal S   
u     = Re(rowSums(eq.Sf*Conj(eq.Sf))/(N-1))
feq.S = kernapply(u, kd, circular=TRUE)
u     = Re(rowSums(ex.Sf*Conj(ex.Sf))/(N-1))
fex.S =	kernapply(u, kd, circular=TRUE)

plot(freq, feq.S[1:512]/fex.S[1:512], type="l", xlab="Frequency (Hz)", ylab="F Statistic", 
      main="Equal S-Spectra")
abline(h=qf(.999, 2*L*(N-1),  2*L*(N-1))) 

# Equal Spectra  
u      = rowSums(eq.Pf*Conj(eq.Sf))/(N-1)
feq.PS = kernapply(u, kd, circular=TRUE)
u      = rowSums(ex.Pf*Conj(ex.Sf)/(N-1))
fex.PS = kernapply(u, kd, circular=TRUE)
fv11   = kernapply(fv11, kd, circular=TRUE)
fv22   = kernapply(fv22, kd, circular=TRUE)
fv12  = kernapply(fv12, kd, circular=TRUE)

Mi = L*(N-1) 
M  = 2*Mi
TS = rep(NA,512)

for (k  in 1:512){
 det.feq.k = Re(feq.P[k]*feq.S[k] - feq.PS[k]*Conj(feq.PS[k]))
 det.fex.k = Re(fex.P[k]*fex.S[k] - fex.PS[k]*Conj(fex.PS[k]))
 det.fv.k  = Re(fv11[k]*fv22[k] - fv12[k]*Conj(fv12[k]))
 
 log.n1 = log(M)*(M*p.dim)
 log.d1 = log(Mi)*(2*Mi*p.dim)
 log.n2 = log(Mi)*2 +log(det.feq.k)*Mi + log(det.fex.k)*Mi  
 log.d2 = (log(M)+log(det.fv.k))*M
 r = 1 - ((p.dim+1)*(p.dim-1)/6*p.dim*(2-1))*(2/Mi - 1/M)
 TS[k] = -2*r*(log.n1+log.n2-log.d1-log.d2)   
}

plot(freq, TS, type="l", xlab="Frequency (Hz)", ylab="Chi-Sq Statistic", main="Equal Spectral Matrices")
abline(h = qchisq(.9999, p.dim^2)) 


Example 7.11

P     = 1:1024
S     = P+1024
mag.P = log10(apply(eqexp[P,],2,max) - apply(eqexp[P,],2,min))
mag.S = log10(apply(eqexp[S,],2,max) - apply(eqexp[S,],2,min))
eq.P  = mag.P[1:8]
eq.S  = mag.S[1:8]
ex.P  = mag.P[9:16]
ex.S =  mag.S[9:16]
NZ.P  = mag.P[17]
NZ.S  = mag.S[17] 

# Compute linear discriminant function 
cov.eq     = var(cbind(eq.P, eq.S))
cov.ex     = var(cbind(ex.P, ex.S))
cov.pooled = (cov.ex + cov.eq)/2

means.eq   =  colMeans(cbind(eq.P,eq.S)); 
means.ex   =  colMeans(cbind(ex.P,ex.S))
slopes.eq  = solve(cov.pooled, means.eq)
inter.eq   = -sum(slopes.eq*means.eq)/2
slopes.ex  = solve(cov.pooled, means.ex)
inter.ex   = -sum(slopes.ex*means.ex)/2
d.slopes   = slopes.eq - slopes.ex
d.inter    = inter.eq - inter.ex

# Classify new observation 
new.data = cbind(NZ.P, NZ.S)

d = sum(d.slopes*new.data) + d.inter
post.eq = exp(d)/(1+exp(d))

# Print (disc function, posteriors) and plot results  
cat(d.slopes[1], "mag.P +" , d.slopes[2], "mag.S +" , d.inter,"\n")  
cat("P(EQ|data) =", post.eq,  "  P(EX|data) =", 1-post.eq, "\n" )    

plot(eq.P, eq.S, xlim=c(0,1.5), ylim=c(.75,1.25), xlab="log mag(P)", ylab ="log mag(S)", pch = 8, 
      cex=1.1, lwd=2, main="Classification Based on Magnitude Features")
 points(ex.P, ex.S, pch = 6, cex=1.1, lwd=2)
 points(new.data, pch = 3, cex=1.1, lwd=2)
 abline(a = -d.inter/d.slopes[2], b = -d.slopes[1]/d.slopes[2])
 text(eq.P-.07,eq.S+.005, label=names(eqexp[1:8]), cex=.8)
 text(ex.P+.07,ex.S+.003, label=names(eqexp[9:16]), cex=.8)
 text(NZ.P+.05,NZ.S+.003, label=names(eqexp[17]), cex=.8)
 legend("topright",c("EQ","EX","NZ"),pch=c(8,6,3),pt.lwd=2,cex=1.1)

# Cross-validation 
all.data = rbind(cbind(eq.P,eq.S), cbind(ex.P,ex.S))
post.eq <- rep(NA, 8) -> post.ex 

for(j in 1:16) {
 if (j <= 8) {samp.eq = all.data[-c(j, 9:16),];  samp.ex = all.data[9:16,]}
 if (j > 8)  {samp.eq = all.data[1:8,];    samp.ex = all.data[-c(j, 1:8),]}
 
 df.eq      = nrow(samp.eq)-1;  df.ex = nrow(samp.ex)-1
 mean.eq    = colMeans(samp.eq);  mean.ex = colMeans(samp.ex)
 cov.eq     = var(samp.eq);  cov.ex = var(samp.ex)
 cov.pooled = (df.eq*cov.eq + df.ex*cov.ex)/(df.eq + df.ex)
 slopes.eq  = solve(cov.pooled, mean.eq)
 inter.eq   = -sum(slopes.eq*mean.eq)/2
 slopes.ex  = solve(cov.pooled, mean.ex)
 inter.ex   = -sum(slopes.ex*mean.ex)/2
 d.slopes   = slopes.eq - slopes.ex
 d.inter    = inter.eq - inter.ex
 
 d = sum(d.slopes*all.data[j,]) + d.inter
 if (j <= 8) post.eq[j] = exp(d)/(1+exp(d))
 if (j > 8) post.ex[j-8] = 1/(1+exp(d))  
}

Posterior = cbind(1:8, post.eq, 1:8, post.ex)
colnames(Posterior) = c("EQ","P(EQ|data)","EX","P(EX|data)")
# results from cross-validation 
round(Posterior, 3)   


Example 7.12

P = 1:1024
S = P+1024
p.dim = 2
n =1024

eq = as.ts(eqexp[,1:8])
ex = as.ts(eqexp[,9:16])
nz = as.ts(eqexp[,17])
f.eq <- array(dim=c(8,2,2,512)) -> f.ex 
f.NZ = array(dim=c(2,2,512))

# below calculates determinant for 2x2 Hermitian matrix
det.c = function(mat){return(Re(mat[1,1]*mat[2,2]-mat[1,2]*mat[2,1]))}
L = c(15,13,5)  # for smoothing 
for (i in 1:8){     # compute spectral matrices
 f.eq[i,,,] = mvspec(cbind(eq[P,i],eq[S,i]), spans=L, taper=.5, plot=FALSE)$fxx
 f.ex[i,,,] = mvspec(cbind(ex[P,i],ex[S,i]), spans=L, taper=.5, plot=FALSE)$fxx
}
u = mvspec(cbind(nz[P],nz[S]), spans=L, taper=.5) 
f.NZ = u$fxx	 
bndwidth = u$bandwidth*40  # about .75 Hz
fhat.eq = apply(f.eq, 2:4, mean)  # average spectra
fhat.ex = apply(f.ex, 2:4, mean)

# plot the average spectra
par(mfrow=c(2,2), mar=c(3,3,2,1),  mgp = c(1.6,.6,0))
Fr = 40*(1:512)/n
plot(Fr, Re(fhat.eq[1,1,]), type="l", xlab="Frequency (Hz)", ylab="")
plot(Fr, Re(fhat.eq[2,2,]), type="l", xlab="Frequency (Hz)", ylab="")
plot(Fr, Re(fhat.ex[1,1,]), type="l", xlab="Frequency (Hz)", ylab="")
plot(Fr, Re(fhat.ex[2,2,]), type="l", xlab="Frequency (Hz)", ylab="")
mtext("Average P-spectra", side=3, line=-1.5, adj=.2, outer=TRUE)
mtext("Earthquakes", side=2, line=-1, adj=.8,  outer=TRUE)
mtext("Average S-spectra", side=3, line=-1.5, adj=.82, outer=TRUE)
mtext("Explosions", side=2, line=-1, adj=.2, outer=TRUE)

dev.new()
par(fig = c(.75, 1, .75, 1), new = TRUE)
ker = kernel("modified.daniell", L)$coef; ker = c(rev(ker),ker[-1])
plot((-33:33)/40,ker,type="l",ylab="",xlab="",cex.axis=.7,yaxp=c(0,.04,2))

# choose alpha
Balpha = rep(0,19) 
for (i in 1:19){  alf=i/20
for (k in 1:256) {  	
Balpha[i]= Balpha[i] + Re(log(det.c(alf*fhat.ex[,,k] + (1-alf)*fhat.eq[,,k])/ det.c(fhat.eq[,,k]))- 
                           alf*log(det.c(fhat.ex[,,k])/det.c(fhat.eq[,,k])))} }
alf = which.max(Balpha)/20   # = .4  

# calculate information criteria          
rep(0,17) -> KLDiff -> BDiff -> KLeq -> KLex -> Beq -> Bex
for (i in 1:17){
 if (i <= 8) f0 = f.eq[i,,,]
 if (i > 8 & i <= 16) f0 = f.ex[i-8,,,]
 if (i == 17) f0 = f.NZ
 for (k in 1:256) {    # only use freqs out to .25
  tr = Re(sum(diag(solve(fhat.eq[,,k],f0[,,k]))))
  KLeq[i] = KLeq[i] + tr + log(det.c(fhat.eq[,,k])) - log(det.c(f0[,,k]))
  Beq[i] =  Beq[i] + Re(log(det.c(alf*f0[,,k]+(1-alf)*fhat.eq[,,k])/det.c(fhat.eq[,,k])) - 
                         alf*log(det.c(f0[,,k])/det.c(fhat.eq[,,k])))
  tr = Re(sum(diag(solve(fhat.ex[,,k],f0[,,k]))))
  KLex[i] = KLex[i] + tr +  log(det.c(fhat.ex[,,k])) - log(det.c(f0[,,k]))
  Bex[i] = Bex[i] + Re(log(det.c(alf*f0[,,k]+(1-alf)*fhat.ex[,,k])/det.c(fhat.ex[,,k])) - 
                       alf*log(det.c(f0[,,k])/det.c(fhat.ex[,,k]))) 
  }
KLDiff[i] = (KLeq[i] - KLex[i])/n 
BDiff[i] =  (Beq[i] - Bex[i])/(2*n) 
}

x.b = max(KLDiff)+.1; x.a = min(KLDiff)-.1
y.b = max(BDiff)+.01; y.a = min(BDiff)-.01

dev.new()                                 
plot(KLDiff[9:16], BDiff[9:16], type="p", xlim=c(x.a,x.b), ylim=c(y.a,y.b), cex=1.1, lwd=2, 
      xlab="Kullback-Leibler Difference",ylab="Chernoff Difference", main="Classification 
      Based on Chernoff and K-L Distances", pch=6)
points(KLDiff[1:8], BDiff[1:8], pch=8, cex=1.1, lwd=2)
points(KLDiff[17], BDiff[17],  pch=3, cex=1.1, lwd=2)
legend("topleft", legend=c("EQ", "EX", "NZ"), pch=c(8,6,3), pt.lwd=2)
abline(h=0, v=0, lty=2, col="gray")
text(KLDiff[-c(1,2,3,7,14)]-.075, BDiff[-c(1,2,3,7,14)], label=names(eqexp[-c(1,2,3,7,14)]), cex=.7)
text(KLDiff[c(1,2,3,7,14)]+.075, BDiff[c(1,2,3,7,14)], label=names(eqexp[c(1,2,3,7,14)]), cex=.7)


Example 7.13

library(cluster)
P = 1:1024
S = P+1024
p.dim = 2
n =1024

eq = as.ts(eqexp[,1:8])
ex = as.ts(eqexp[,9:16])
nz = as.ts(eqexp[,17])

f = array(dim=c(17,2,2,512))  
L = c(15,15)   # for smoothing 
for (i in 1:8){     # compute spectral matrices
 f[i,,,] = mvspec(cbind(eq[P,i],eq[S,i]), spans=L, taper=.5, plot=FALSE)$fxx
 f[i+8,,,] = mvspec(cbind(ex[P,i],ex[S,i]), spans=L, taper=.5, plot=FALSE)$fxx 
}
f[17,,,] = mvspec(cbind(nz[P],nz[S]), spans=L, taper=.5, plot=FALSE)$fxx	

# calculate symmetric information criteria 
JD = matrix(0,17,17) 
for (i in 1:16){
 for (j in (i+1):17){	
  for (k in 1:256) {    # only use freqs out to .25
    tr1 = Re(sum(diag(solve(f[i,,,k],f[j,,,k]))))
    tr2 = Re(sum(diag(solve(f[j,,,k], f[i,,,k]))))
    JD[i,j] = JD[i,j] + (tr1 + tr2 - 2*p.dim)
  }
 }
}
 
JD = (JD + t(JD))/n
colnames(JD) = c(colnames(eq), colnames(ex), "NZ") 
rownames(JD) = colnames(JD)
cluster.2 = pam(JD, k = 2, diss = TRUE)  

summary(cluster.2)  # print results
par(mgp = c(1.6,.6,0), cex=3/4, cex.lab=4/3, cex.main=4/3)
clusplot(JD, cluster.2$cluster, col.clus=1, labels=3, lines=0, col.p=1, 
               main="Clustering Results for Explosions and Earthquakes")
text(-7,-.5, "Group I", cex=1.1, font=2)
text(1, 5, "Group II", cex=1.1, font=2)


Example 7.14

n = 128
Per = abs(mvfft(fmri1[,-1]))^2/n

par(mfrow=c(2,4), mar=c(3,2,2,1),  mgp = c(1.6,.6,0), oma=c(0,1,0,0))
for (i in 1:8) plot(0:20, Per[1:21,i], type="l", ylim=c(0,8), main=colnames(fmri1)[i+1], 
        xlab="Cycles", ylab="", xaxp=c(0,20,5))
mtext("Periodogram", side=2, line=-.3, outer=TRUE, adj=c(.2,.8))

fxx = mvspec(fmri1[,-1], kernel("daniell", c(1,1)), taper=.5, plot=FALSE)$fxx
l.val = rep(NA,64)
for (k in 1:64) {
 u = eigen(fxx[,,k], symmetric=TRUE, only.values=TRUE)
 l.val[k] = u$values[1]
}

dev.new()
plot(l.val, type="l", xaxp=c(0,64,8), xlab="Cycles (Frequency x 128)", ylab="First Principal Component")
axis(1, seq(4,60,by=8), labels=FALSE)

# at freq k=4
u = eigen(fxx[,,4], symmetric=TRUE)
lam = u$values 
evec = u$vectors
lam[1]/sum(lam) # % of variance explained
sig.e1 = matrix(0,8,8)
for (l in 2:5){  # last 3 evs are 0
 sig.e1 = sig.e1 + lam[l]*evec[,l]%*%Conj(t(evec[,l]))/(lam[1]-lam[l])^2
}
sig.e1 = Re(sig.e1)*lam[1]*sum(kernel("daniell", c(1,1))$coef^2)
p.val = round(pchisq(2*abs(evec[,1])^2/diag(sig.e1), 2, lower.tail=FALSE), 3)
cbind(colnames(fmri1)[-1], abs(evec[,1]), p.val) # print table values


Example 7.15

bhat = sqrt(lam[1])*evec[,1]
Dhat = Re(diag(fxx[,,4] - bhat%*%Conj(t(bhat))))
res = Mod(fxx[,,4] - Dhat - bhat%*%Conj(t(bhat)))


Example 7.16

gr = diff(log(ts(econ5, start=1948, frequency=4))) # growth rate
plot(100*gr, main="Growth Rates (%)")

# scale each series to have variance 1
gr = ts(apply(gr,2,scale), freq=4)   # scaling strips ts attributes
L = c(7,7)   # degree of smoothing 
gr.spec = mvspec(gr, spans=L, detrend=FALSE, taper=.25, plot=FALSE)

dev.new()
plot(kernel("modified.daniell", L))  # view the kernel - not shown 

dev.new()
plot(gr.spec, log="no",col=1, main="Individual Spectra", lty=1:5, lwd=2)
legend("topright", colnames(econ5), lty=1:5, lwd=2)  

dev.new()
plot.spec.coherency(gr.spec, ci=NA,  main="Squared Coherencies")

# PCs
n.freq = length(gr.spec$freq)
lam = matrix(0,n.freq,5)
for (k in 1:n.freq) lam[k,] = eigen(gr.spec$fxx[,,k], symmetric=TRUE, only.values=TRUE)$values 

dev.new()
par(mfrow=c(2,1), mar=c(4,2,2,1), mgp=c(1.6,.6,0))
plot(gr.spec$freq, lam[,1], type="l", ylab="", xlab="Frequency", main="First Eigenvalue")
abline(v=.25, lty=2)
plot(gr.spec$freq, lam[,2], type="l", ylab="", xlab="Frequency", main="Second Eigenvalue")
abline(v=.125, lty=2)
e.vec1 = eigen(gr.spec$fxx[,,10], symmetric=TRUE)$vectors[,1] 
e.vec2 =  eigen(gr.spec$fxx[,,5], symmetric=TRUE)$vectors[,2]
round(Mod(e.vec1), 2);  round(Mod(e.vec2), 3) 


Example 7.18

u = factor(bnrf1ebv)  # first, input the data as factors and then 
x = model.matrix(~u-1)[,1:3]  # make an indicator matrix
# x = x[1:1000,]  # select subsequence if desired

Var = var(x)  # var-cov matrix 
xspec = mvspec(x, spans=c(7,7), detrend=FALSE, plot=FALSE)  
fxxr = Re(xspec$fxx)  # fxxr is real(fxx) 

# compute Q = Var^-1/2 
ev = eigen(Var)  
Q = ev$vectors%*%diag(1/sqrt(ev$values))%*%t(ev$vectors)             

# compute spec env and scale vectors     
num     = xspec$n.used  # sample size used for FFT
nfreq   = length(xspec$freq)  # number of freqs used   
specenv = matrix(0,nfreq,1)   # initialize the spec envelope
beta    = matrix(0,nfreq,3)   # initialize the scale vectors
 
for (k in 1:nfreq){ 
  ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num, symmetric=TRUE)   
  specenv[k] = ev$values[1]   # spec env at freq k/n is max evalue
  b = Q%*%ev$vectors[,1]      # beta at freq k/n 
  beta[k,] = b/sqrt(sum(b^2)) # helps to normalize beta
}  
 
# output and graphics  
frequency = xspec$freq
plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)")

# add significance threshold to plot 
m = xspec$kernel$m
etainv = sqrt(sum(xspec$kernel[-m:m]^2))
thresh = 100*(2/num)*exp(qnorm(.9999)*etainv)*rep(1,nfreq)
lines(frequency, thresh, lty="dashed", col="blue")

# details 
output = cbind(frequency, specenv, beta)
colnames(output) = c("freq","specenv", "A", "C", "G")
round(output,3)


Example 7.19

u = arima(diff(log(gnp)), order=c(0,0,2))$resid  # residuals
x = cbind(u, abs(u), u^2)  # transformation set 
Var = var(x)                    
xspec = mvspec(x, spans=c(5,3), taper=.1, plot=FALSE) 
fxxr = Re(xspec$fxx);  ev = eigen(Var)  
Q = ev$vectors%*%diag(1/sqrt(ev$values))%*%t(ev$vectors)               
num = xspec$n.used;  nfreq = length(xspec$freq)        
specenv = matrix(0, nfreq, 1); beta = matrix(0,nfreq,3)           

for (k in 1:nfreq){ 
  ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num)  
  specenv[k] = ev$values[1]            
  b = Q%*%ev$vectors[,1]               
  beta[k,] = b/b[1]  
}                 
# output and graphics 
frequency = xspec$freq
plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)")
output = cbind(frequency, specenv, beta)
colnames(output) = c("freq","specenv","x", "|x|", "x^2")
round(output,4)  # results (not shown)

# plot transformation
b = output[1, 3:5]; g = function(x) b[1]*x+b[2]*abs(x)+b[3]*x^2
curve(g, -.04, .04)


Example 7.20

set.seed(90210)
u = exp(3*sin(2*pi*1:500*.1) + rnorm(500,0,4))   # the data
spec.pgram(u, spans=c(5,3), taper=.5, log="dB")

dev.new()
x = cbind(u, sqrt(u), u^(1/3))           # transformation set 
Var = var(x)                     
xspec = mvspec(x, spans=c(5,3), taper=.5, plot=FALSE)  
fxxr = Re(xspec$fxx)            
ev = eigen(Var)  
Q = ev$vectors%*%diag(1/sqrt(ev$values))%*%t(ev$vectors)               
num = xspec$n.used               
nfreq = length(xspec$freq)         
specenv = matrix(0,nfreq,1)      
beta=matrix(0,nfreq,3)          
 for (k in 1:nfreq){ 
  ev = eigen(2*Q%*%fxxr[,,k]%*%Q/num) 
  specenv[k] = ev$values[1]            
  b = Q%*%ev$vectors[,1]               
  beta[k,] = b/sign(b[1]) 
}

# output and graphics  
frequency = xspec$freq
plot(frequency, 100*specenv, type="l", ylab="Spectral Envelope (%)")
m = xspec$kernel$m
etainv = sqrt(sum(xspec$kernel[-m:m]^2))
thresh = 100*(2/num)*exp(qnorm(.999)*etainv)*rep(1,nfreq)
lines(frequency,thresh, lty="dashed", col="blue")
output = cbind(frequency, specenv, beta)
colnames(output) = c("freq","specenv","x", "sqrt(x)", "x^(1/3)")
round(output,4) 

# plot transform
dev.new()
b = output[50,3:5]
g = function(x) 4.5 + b[1]*x+b[2]*sqrt(x)+b[3]*x^(1/3)
curve(g, 1, 4000)
lines(log(1:4000), lty=2)