## Example 7.14.  


bold = ts(read.table("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/fmri.dat"), frequency = 1)
par(mfrow=c(2,1))
ts.plot(bold[,2:5], lty = 1:4, xlab = "time (x 2 seconds)", ylab = "BOLD")
ts.plot(bold[,6:9], lty = 1:4, xlab = "time (x 2 seconds)", ylab = "BOLD")
title(main = "\nBlood oxygenation-level dependent signal intensity", outer = T)

data = matrix(nrow = nrow(bold), ncol = 8)  # Standardized values of BOLD; common in PCA
for (j in 2:9) data[,j-1] = bold[,j]/sqrt(var(bold[,j]))


L = 5
kd=kernel(coef = c(3,2,1)/9, m=2)
kd


# SPEC[i,j,k] is the spectrum between the i-th and j-th series at frequency 2k/n'=k/128:
SPEC = array(dim = c(ncol(data),ncol(data),nextn(nrow(data))/2)) 
for(i in 1:ncol(data)) { 
for (j in i:ncol(data)) {
#	power = spec.pgram(data[,c(i,j)], kernel("daniell",(L-1)/2), plot = F) # Uniform smoothing
	power = spec.pgram(data[,c(i,j)], kernel(coef = c(3,2,1)/9, m=2), plot = F) # Nonuniform smoothing
	SPEC[i,i, ] = power$spec[,1]
	SPEC[j,j, ] = power$spec[,2]
	coh.ij = power$coh 
	phase.ij = power$phase 
	SPEC[i,j, ] = sqrt(coh.ij*power$spec[,1]*power$spec[,2])*exp(1i*phase.ij)
	SPEC[j,i, ] = Conj(SPEC[i,j, ])
	}}


## Plot the 8 individual spectra:
win.graph()
par(mfrow = c(2,4))
for (j in 1:8) {
plot(c(1:64)/128, Re(SPEC[j,j,]), xlab = "cycles per unit time", ylab = paste("Series", j), type = "l")
abline(v=4/128, col = 2)
}
title(main = "\nSample spectra; vertical lines at frequency omega = 4/128", outer = T)


# Compute and plot the first pc's at each frequency 1/128, ... , 64/128:
lambda = vector(length = 64)
for (j in 1:64) {
f = SPEC[,,j]
eig = eigen(f, only.values = T)
lambda[j] = eig$values[1]
}
win.graph()
plot(c(1:64)/128, lambda, xlab = "cycles per unit time", ylab = "First PC", type = "l")
abline(v=4/128, col = 2)
title(main = "\nSpectral density of first PC at each frequency 1/128, ..., 64/128", outer = T)


f.32 = SPEC[,,4]  # Spectral density  matrix at omega = 4/128 = 1/32
eig = eigen(f.32)
cat("Eigenvalues at frequency 4/128 are", round(eig$values, 2), "\n")
cat("Contribution of the first PC is", 100*round(eig$values[1]/sum(eig$values), 4), "%", "\n")
cat("Coefficients of first PC are", "\n", round(eig$vectors[,1],2), "\n", 
    "with moduli", round(abs(eig$vectors[,1]),2), "\n")

x = t(data)
e1 = eig$vectors[,1]
pc1 = t(Conj(e1))%*%x
pc1 = as.vector(pc1)

win.graph()
par(mfrow=c(2,2))
ts.plot(Re(pc1), main = "Real part")
ts.plot(Im(pc1), main = "Imaginary part")
ts.plot(Mod(pc1), main = "Modulus")
phase = atan(Im(pc1)/Re(pc1))
ts.plot(phase, main = "Phase")
title(main = "First pc at signal frequency 1/32", outer = TRUE, line=-1)


e2 = eig$vectors[,2]
pc2 = t(Conj(e2))%*%x
pc2 = as.vector(pc2)

win.graph()
par(mfrow=c(2,2))
ts.plot(Re(pc2), main = "Real part")
ts.plot(Im(pc2), main = "Imaginary part")
ts.plot(Mod(pc2), main = "Modulus")
phase = atan(Im(pc2)/Re(pc2))
ts.plot(phase, main = "Phase")
title(main = "Second pc at signal frequency 1/32", outer = TRUE, line=-1)