FDJ2<-function(x,y, lagmax=10, limmin=-0.5, limpl=0.5) {
# FD jackknife, supposedly small sample stable
# Ref [8]
# linear interpolation of crossperiodograms
# example: FDJ2(x=VI,y=EC)
# required arguments; x and y
#
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
# unbiased estimates of the cross-covariances
x3pos<-matrix(0,nrow=lagmax,ncol=1)
x3neg<-matrix(0,nrow=lagmax,ncol=1)
for (k in (1:lagmax)) {
  x3pos[k]<-mean(x1[(1+u1[k]):length(x1)]*x2[1:(length(x2)-u1[k])])
  x3neg[k]<-mean(x1[1:(length(x1)-u1[k])]*x2[(1+u1[k]):length(x2)])
}
c33pos<-matrix(0, nrow=lagmax, ncol=floor((length(x)-1)/2))
c33neg<-c33pos
ntild<-floor((length(x)-1)/2)
nprim<-2*length(x)
#taper and pad
x1p<-c(rep(0,floor(length(x)/2)),x1,rep(0,ceiling(length(x)/2)))
x1p[1:floor(0.1*length(x1p))]<-(1-cos(pi*(time(x1p[1:floor(0.1*length(x1p))])-0.5)/(0.1*length(x1p))))*x1p[1:floor(0.1*length(x1p))]
x1p[(floor(0.9*length(x1p))+1):length(x1p)]<-(1-cos(pi*(length(x1p)-time(x1p[(floor(0.9*length(x1p))+1):length(x1p)])+0.5)/(0.1*length(x1p))))*x1p[(floor(0.9*length(x1p))+1):length(x1p)]
d11<-fft(x1p)
x2p<-c(rep(0,floor(length(x)/2)),x2,rep(0,ceiling((length(x))/2)))
x2p[1:floor(0.1*length(x2p))]<-(1-cos(pi*(time(x2p[1:floor(0.1*length(x2p))])-0.5)/(0.1*length(x2p))))*x2p[1:floor(0.1*length(x2p))]
x2p[(floor(0.9*length(x2p))+1):length(x2p)]<-(1-cos(pi*(length(x2p)-time(x2p[(floor(0.9*length(x2p))+1):length(x2p)])+0.5)/(0.1*length(x2p))))*x2p[(floor(0.9*length(x2p))+1):length(x2p)]
d22<-fft(x2p)
#
I11<-abs(d11)^2/(2*pi*length(x1))
I22<-abs(d22)^2/(2*pi*length(x1))
I12<-d11*Conj(d22)/(2*pi*length(x1))
I21<-Conj(d11)*d22/(2*pi*length(x1))
c3pos<-Re(fft(2*pi*I12/length(x1p), inverse=T)[1:lagmax])
c3neg<-Re(fft(2*pi*I21/length(x1p), inverse=T)[1:lagmax])
fourfreq<-seq(0,(length(x1p)-1))*2*pi/length(x1p)
for (k in (1:ntild)) {
  I12j<-I12
  I21j<-I21
  # ntild jacknife terms
  kk<-c(2*k,2*k+1,2*k+2)
  omega1<-(rep(fourfreq[2*k+3],3)-fourfreq[kk])/(fourfreq[2*k+3]-fourfreq[2*k-1])  
  omega2<-c(1,1,1)-omega1
  # linear interpolation 
  I12j[kk]<-omega1*I12[2*k-1]+omega2*I12[2*k+3]
  I21j[kk]<-omega1*I21[2*k-1]+omega2*I21[2*k+3]
  kkn<-sort(rep(length(x1p),3)-(kk-rep(1,3)))
  omega1<-(rep(fourfreq[kkn[3]+1],3)-fourfreq[kkn])/(fourfreq[kkn[3]+1]-fourfreq[kkn[1]-1])  
  omega2<-c(1,1,1)-omega1
  I12j[kkn]<-omega1*I12[kkn[1]-1]+omega2*I12[kkn[3]+1]
  I21j[kkn]<-omega1*I21[kkn[1]-1]+omega2*I21[kkn[3]+1]
  c33pos[,k]<-Re(fft(2*pi*I12j/length(x1p), inverse=T)[1:lagmax])
  c33neg[,k]<-Re(fft(2*pi*I21j/length(x1p), inverse=T)[1:lagmax])
}
c33pos<-apply(c33pos,2,function(c33pos,vec) {c33pos-vec}, vec=c3pos)
c33neg<-apply(c33neg,2,function(c33neg,vec2) {c33neg-vec2}, vec2=c3neg)
c33pos<-c33pos*c33pos
c33neg<-c33neg*c33neg
c3p<-(length(x)/ntild)*0.623*apply(c33pos,1,sum)
c3n<-(length(x)/ntild)*0.623*apply(c33neg,1,sum)
# plot
par(mfrow=c(2,2))
uppb<-2*sqrt(c3p)
lowb<-(-uppb)
plot(u1,x3pos, 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<-2*sqrt(c3n)
lowb<-(-uppb)
plot(u2,x3neg, 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")
}