cfband<-function(x,y,lagmax=10, limmin=-0.5,limpl=0.5) {
# The Direct Method as presented by Brillinger [4]
# required arguments; x,y. lagmax is the number of cross-cov estimated
# limmin and limpl are the axis in the produced plots.
# example of use: cfband(VI,EC,lagmax=25, limmin=-1, limpl=1)
c3pos<-matrix(0,nrow=lagmax, ncol=1)
f3pos<-c3pos
c3neg<-matrix(0,nrow=lagmax, ncol=1)
f3neg<-c3neg
# Standardize series
x1<-(as.vector(x)-mean(as.vector(x)))/sqrt(var(as.vector(x)))
x2<-(as.vector(y)-mean(as.vector(y)))/sqrt(var(as.vector(y)))
u1<-seq(0,(lagmax-1))
u2<--u1
# for each lag calculate the product series and estimate its mean
for (k in (1:lagmax)) {
  x3<-x1[(1+u1[k]):length(x1)]*x2[1:(length(x2)-u1[k])]
  c3pos[k]<-mean(x3)
  #pad
  x3<-ts(c(x3,rep(0,(length(x1)-length(x3)))))
  x3p<-x3
  #taper
  x3p[1:floor(0.1*length(x3))]<-(1-cos(pi*(time(x3p[1:floor(0.1*length(x3))])-0.5)/(0.1*length(x3))))*x3p[1:floor(0.1*length(x3))]
  x3p[(floor(0.9*length(x3))+1):length(x3)]<-(1-cos(pi*(length(x3)-time(x3p[(floor(0.9*length(x3))+1):length(x3)])+0.5)/(0.1*length(x3))))*x3p[(floor(0.9*length(x3))+1):length(x3)]
  #spectrum estimate
  f33<-fft(x3p)
  f33<-c(f33[(length(f33)/2+2):length(f33)],f33[1:(length(f33)/2+1)])
  f33p<-(Mod(f33)^2/(2*pi*length(x3)))
  f3pos[k]<-mean(f33p[(length(f33)/2-5):(length(f33)/2+5)])
  # neg lag
  x3<-x1[1:(length(x1)-u1[k])]*x2[(1+u1[k]):length(x2)]
  c3neg[k]<-mean(x3)
  x3<-ts(c(x3,rep(0,(length(x1)-length(x3)))))
  x3p<-x3
  x3p[1:floor(0.1*length(x3))]<-(1-cos(pi*(time(x3p[1:floor(0.1*length(x3))])-0.5)/(0.1*length(x3))))*x3p[1:floor(0.1*length(x3))]
  x3p[(floor(0.9*length(x3))+1):length(x3)]<-(1-cos(pi*(length(x3)-time(x3p[(floor(0.9*length(x3))+1):length(x3)])+0.5)/(0.1*length(x3))))*x3p[(floor(0.9*length(x3))+1):length(x3)]
  f33<-fft(x3p)
  f33<-c(f33[(length(f33)/2+2):length(f33)],f33[1:(length(f33)/2+1)])
  f33p<-(Mod(f33)^2/(2*pi*length(x3)))
  f3neg[k]<-mean(f33p[(length(f33)/2-5):(length(f33)/2+5)])
}
#plots
par(mfrow=c(2,2))
#upper and lower bounds
uppb<-qt(0.975,22)*sqrt(2*pi*f3pos/length(x1))
lowb<-(-qt(0.975,22)*sqrt(2*pi*f3pos/length(x1)))
plot(u1,c3pos, type="h",ylim=c(limmin,limpl), ylab="cross-cov", xlab="tao")
lines(u1,uppb,lty=2)
lines(u1,lowb, lty=2)
abline(h=0)
title(main="r_{Y,X}(tao), tao > 0")
e1<-ts(x1)
e2<-ts(x2)
a<-acf(ts.union(e1,e2),plot=F, lag.max=lagmax)
plot(u1,a$acf[1:lagmax,1,2], type="h",ylim=c(limmin,limpl), ylab="cross-cov", xlab="tao")
abline(h=2/sqrt(length(e1)), lty=2)
abline(h=(-2/sqrt(length(e1))), lty=2)
abline(h=0)
title(main="standard est. method")
#
uppb<-qt(0.975,22)*sqrt(2*pi*f3neg/length(x1))
lowb<-(-qt(0.975,22)*sqrt(2*pi*f3neg/length(x1)))
plot(u2,c3neg, type="h",ylim=c(limmin,limpl), ylab="cross-cov", xlab="tao")
lines(u2,uppb,lty=2)
lines(u2,lowb, lty=2)
abline(h=0)
title(main="r_{Y,X}(tao), tao < 0")
plot(u2,a$acf[1:lagmax,2,1], type="h",ylim=c(limmin,limpl), ylab="cross-cov", xlab="tao")
abline(h=2/sqrt(length(e1)), lty=2)
abline(h=(-2/sqrt(length(e1))), lty=2)
abline(h=0)
title(main="standard est. method")
}