Lagged and Residual Coherence
================================

LagCoh <- function(XY,u){ 
#Input XY Two Dim Time Series.
#Input Integer Lag u.
#Output: Lin. Coh.^2, Lagged Coh.
#Modif. Daniell (3,5,7,9) Spectra.
######################
X <- ts(XY[,1])
Y <- ts(XY[,2])
X <- X-mean(X)
Y <- Y-mean(Y)
#######################
#Two Dim Series XY Of Length: length(X)-u
XY <- cbind(X[1:(length(X)-u)], Y[1:(length(Y)-u)])
XY <- ts(XY)
########################
#Define Lag Process
Xu <- X[1:(length(X)-u)]*X[(1+u):length(X)]
Xu <- Xu-mean(Xu)
########################
#Two Dim Series XuY Of Length length(X)-u
XuY <- cbind(Xu,Y[1:(length(Y)-u)])
XuY <- ts(XuY)
########################
#Two Dim Series XuX Of Length length(X)-u
XuX <- cbind(Xu,X[1:(length(X)-u)])
XuX <- ts(XuX)
########################
#XY Spectra:
XYsp <- spec.pgram(XY, spans=c(3,5,7,9))
FX <- XYsp$spec[,1]
FY <- XYsp$spec[,2]
FX <- 10^(FX/10.) #Unlog!!!
FY <- 10^(FY/10.) #Unlog!!!
phaseXY <- XYsp$phase
cohXY <- sqrt(XYsp$coh) 
FXY <- sqrt(FX*FY)*cohXY*complex(real=cos(phaseXY), imag=sin(phaseXY))
########################
#XuY Spectra:
XuYsp <- spec.pgram(XuY, spans=c(3,5,7,9))
FXu <- XuYsp$spec[,1]
FXu <- 10^(FXu/10.) #Unlog!!!
phaseXuY <- XuYsp$phase
cohXuY <- sqrt(XuYsp$coh) 
FXuY <- sqrt(FXu*FY)*cohXuY*complex(real=cos(phaseXuY), imag= sin(phaseXuY))
########################
#XuX Spectra:
XuXsp <- spec.pgram(XuX, spans=c(3,5,7,9))
phaseXuX <- XuXsp$phase
cohXuX <- sqrt(XuXsp$coh) 
FXuX <- sqrt(FXu*FX)*cohXuX*complex(real=cos(phaseXuX), imag= sin(phaseXuX))
########################
#Lagged Coherence:
B <- (FX*FXuY-FXuX*FXY)/(FX*FXu-abs(FXuX)^2)
A <- (FXu -(abs(FXuX)^2)/FX)*(abs(B)^2)/FY
LaggedCoh <- cohXY^2 + A 
#AllCoh <- cbind(cohXY^2,LaggedCoh)
#tsplot(AllCoh, type="l",main="Sq. Lin and Lag Coh")
plot(2*pi*XYsp$freq, cohXY^2,type="l",xlab="w",ylab="",ylim=c(0,1.3))
lines(2*pi*XYsp$freq, LaggedCoh,lty=2)
legend(0,1.2,c("Sq. Coh","Lag. Coh"), lty=1:2, bty="n")
title("Squared and Lagged Coherence",cex=1.4)
MaxA <- max(A)
list(Max_Resid_Coh = MaxA)
}


x <- rnorm(1000) ###<---- Better!!!!


lag1 <- c(0,x[-1000])
lag3 <- c(0,0,0,x[-c(998,999,1000)])
lag5 <- c(0,0,0,0,0,x[-c(996,997,998,999,1000)])

###Use lag processes at lags: 0,1,3,5 and linear part:
y <- 0.1*x + 0.3*x^2 + x*lag1 + x*lag3 + 5*x*lag5 + 3*runif(1000,-0.5,0.5)
xy <- cbind(x,y)

LagCoh(xy,0)
LagCoh(xy,1)
LagCoh(xy,5)


## Barplot of max residual coherence
A <- 0
for( u in 0:15){
ttt <- LagCoh(xy,u)
A[u+1] <- ttt$Max_Resid_Coh}
names(A) <- c(0:15)
barplot(A, xlab="u",
cex.axis = par("cex.axis"),border = "dark blue",col=rainbow(20))