########################################################

spec.trap<-function(x){ # 2011 paper trapezoid estimator \hat f 
n<-length(x)
m<- emp.bw(x)
L <- lambda(m)
xcov <- acf(x, lag.max = n, type = "covariance", plot = F)$acf
xcov <- as.vector(xcov) # xcov[1] is lag 0 !
taper<-c(1,L,c(1:(2*n))*0)

taper<- taper[1:n]
f<- Re(fft(taper*xcov)) /(2*pi)
# length(f)= n; it represents Fourier freq 0,2pi/n,2pi*2/n up to 2pi
f[f<0]<- f[f<0]*0 # take positive part
f }
 
#####################################################
spec.NL<-function(x, zeta=2){ # 2011 nonlinear \tilde f 
n<-length(x)
n2<- floor(n/2)
xcov <- acf(x, lag.max=n, type = "covariance", plot = F)$acf
xcov <- as.vector(xcov) # xcov[1] is lag 0 !
Tn<- zeta*2*xcov[1]*sqrt(log(n, 10)/n)

xcov[n2:n]<- xcov[n2:n]*0  # just use first n/2 acf's

xcov[abs(xcov) < Tn] <-xcov[abs(xcov) < Tn]*0 
f<- Re(fft(xcov)) /(2*pi)

# length(f)= n; it represents Fourier freq 0,2pi/n,2pi*2/n up to 2pi
f[f<0]<- f[f<0]*0 # take positive part
f }
 
#####################################################
spec.NL2<-function(x, delta=0){ # 2011 hybrid trapezoid \check f 
n<-length(x)
m<- emp.bw(x)
L <- lambda(m)
xcov <- acf(x, lag.max = n, type = "covariance", plot = F)$acf
xcov <- as.vector(xcov) # xcov[1] is lag 0 !
Tn<- 2*xcov[1]*sqrt(log(n, 10)/n)

xcov[abs(xcov) < delta*Tn] <-xcov[abs(xcov) < delta*Tn]*0 
taper<-c(1,L,c(1:(2*n))*0)

taper<- taper[1:n]
f<- Re(fft(taper*xcov)) /(2*pi)
# length(f)= n; it represents Fourier freq 0,2pi/n,2pi*2/n up to 2pi
f[f<0]<- f[f<0]*0 # take positive part
f }
 
########################################################

emp.bw <- function(x, thresh)
{  
# Empirical bandwidth choice for flat-top spectral density
# from Politis (2003); R function courtesy of Tim McMurry.  
# function returns a value AFTER which the acf is zero
# i.e. for white noise it returns 0
# for MA(1) it returns 1
        n <- length(x)
        Kn <- 5 # as in paper
        # alternatively: Kn<- max(5, 3*sqrt(log(n, 10)-1))
        c <- 2  # as in paper
        
        if(missing(thresh))
                thresh <- c*sqrt(log(n, 10)/n)
                
        ac <- as.vector(acf(x, type="correlation", plot=FALSE, 
                lag.max = floor(n/2))$acf)
        l <- length(ac)
        
        pos <- 1
        
        while(pos < n/2)
        {
                npos <- match(TRUE, abs(ac[pos:l]) < thresh)
                pos <- pos+npos-1

                if(all(abs(ac[pos:(pos+Kn-1)]) < thresh)){
                        return(pos-2)
                } else {
                        npos <- match(FALSE, abs(ac[pos:(pos+Kn-1)]) < thresh)
                        pos <- pos + npos - 1
                }
        }
        stop("No bandwidth found")      
}

########################################################