statboot<-function(p=0.05,B=500,x,y, lagmax=10, limmin=-0.5,limpl=0.5) {
# Romano and Politis 91
# for a chosen p a stationary bootstrap sample is generated Ref [10]
# required arguments; x and y
# example of use statboot(x=VI,y=EC, B=250, p=0.1)
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)])
}
c3pos<-matrix(0, nrow=lagmax, ncol=B)
c3neg<-c3pos
# bootstrap
for (l in (1:B)) {
  xboot<-x1
  yboot<-x2
  N<-seq(1:length(x))
  pp<-rbinom(length(x),1,(1-p))
  pp[1]<-0
  for (ll in (1:length(x))) {
     if (pp[ll]==0) {
        ind<-sample(N,1)
        xboot[ll]<-x1[ind]
        yboot[ll]<-x2[ind] }
     if (pp[ll]==1) {
        ind<-ind+1
        if (ind>length(x)) {
           ind<-sample(N,1)
	 }
        xboot[ll]<-x1[ind]
        yboot[ll]<-x2[ind] }
   }
  for (k in (1:lagmax)) {
   c3pos[k,l]<-mean(xboot[(1+u1[k]):length(xboot)]*yboot[1:(length(yboot)-u1[k])])
   c3neg[k,l]<-mean(xboot[1:(length(xboot)-u1[k])]*yboot[(1+u1[k]):length(yboot)])
   # unbiased bootstrap estimates of cov(k)
   }
}
c3pos<-sqrt(length(x))*abs(apply(c3pos,2,function(c3pos,vec) {c3pos-vec}, vec=x3pos))
c3neg<-sqrt(length(x))*abs(apply(c3neg,2,function(c3neg,vec) {c3neg-vec}, vec=x3neg))
c3p<-apply(c3pos,1,quantile, probs=0.95)
c3n<-apply(c3neg,1,quantile, probs=0.95)
#
par(mfrow=c(2,2))
uppb<-c3p/sqrt(length(x))
lowb<-(-c3p/sqrt(length(x)))
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<-c3n/sqrt(length(x))
lowb<-(-c3n/sqrt(length(x)))
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")
}