#Estimates the n x n covariance matrix and its inverse
CovMat <- function(x, pos.def = TRUE, inverse=TRUE) {
        kappaTrap <- function(x, l){
                if(l > 0) {
                        ifelse(abs(x/l) > 2, 0,
                                ifelse(abs(x/l) < 1, 1, 2 - abs(x/l)))
                } else {
                        c(1, rep(0, length(x)-1))
                }
        }

        pos.part <- function(x, epsilon) {
                ifelse(x > epsilon, x, epsilon)
        }

        adap.bw <- function(x, Kn = 5, c = 2)
        {
                n <- length(x)
                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")
        }

        form.matrix <- function(n, diags)
        {
                mat <- matrix(0,nrow=n, ncol=n)
                len <- length(diags)

                for(i in 1L:n)
                {
                        for(j in 1L:n)
                        {
                                if(abs(i-j) + 1 > len)
                                        mat[i,j] <- 0
                                else
                                        mat[i,j] <- diags[abs(i-j)+1]
                        }
                }
                return(mat)
        }

        n <- length(x)
        x.acf <- as.vector(acf(x, type="covariance", plot=FALSE, lag.max=n)$acf)
        l <- adap.bw(x)

        diags <- kappaTrap((0:(length(x.acf)-1)),l)*x.acf
        Cov.mat <- form.matrix(n, diags)

        if(!pos.def) return(list(S=Cov.mat, Sinv = NA))

        Cov.eig <- eigen(Cov.mat, symmetric=TRUE)
        eig.valsmod <- pos.part(Cov.eig$values, diags[1]/n)
        Spd <- Cov.eig$vectors %*% diag(eig.valsmod) %*%
                t(Cov.eig$vectors)

        if(!inverse) return(list(S=Spd, Sinv = NA))

        Spdinv <- Cov.eig$vectors %*% diag(1/eig.valsmod) %*%
                t(Cov.eig$vectors)

        return(list(S = Spd, Sinv = Spdinv))
}


#LPB takes the original time series (x) as input and returns an n x R matrix of resamples,
#with each column a bootstrap sample.


LPB <- function(x, R) {
        xbar <- mean(x)
        x <- x - xbar
        n <- length(x)

        SSinv <- CovMat(x)

        Ssqrt <- t(chol(SSinv$S))
        Ssqrtinv <- forwardsolve(Ssqrt, diag(1,n))

        Z <- Ssqrtinv %*% x
        Z <- (Z- mean(Z))/sqrt(((n-1)/n)*var(as.vector(Z)))

        onebs <- function() {
                Zst <- sample(Z, n, replace=TRUE)
                xbar + Ssqrt %*% Zst
        }

        replicate(R, onebs())
}