R/liftHurstC.R

In nunesmatt/CliftLRD: Complex-Valued Wavelet Lifting Estimators of the Hurst Exponent for Irregularly Sampled Time Series

liftHurstC <-
function (x, grid = 1:length(x), model = "FGN", ntraj = 50, 
    cutoffs = 0, cut.fine = TRUE, efun = meanmoC, afun = idj, 
    altype = 1, tail = TRUE, normalise = TRUE, level = 0.05, 
    bc = TRUE, vc = TRUE, jsc = TRUE, BHonly = TRUE, ...) 
{
    levvar <- function(nj, max = 5000) {
        z2njh <- sum(1/((0:max) + nj/2)^2)
        v <- z2njh/(log(2))^2
        v
    }
    gj <- function(nj) {
        	bias <- digamma(nj/2)/log(2) - log2(nj/2)
        bias
    }
    if (!is.matrix(x)) {
        x <- matrix(x)
    }
    if (ncol(x) == 2) {
        grid <- x[, 1]
        x <- x[, 2]
    }
    n <- length(x)
    trajmat <- matrix(0, ntraj, n - 2)
    beta <- NULL
    betamat <- NULL
    emat <- scmat <- NULL
    if (bc) {
        cat("doing bias correction...\n")
    }
    if (vc) {
        cat("doing weighted regression...\n")
    }
    if (jsc) {
        cat("doing j* coefficient computation...\n")
    }


        for (i in 1:ntraj) {
            trajmat[i, ] <- sample(1:n, n - 2, replace = FALSE)
        }
        cat("generated trajectories.\n")
        for (k in 1:ntraj) {
            if ((k%%5) == 0) {
                cat(k, "...")
            }
            xlift <- fwtnppermC(grid, x, mod = trajmat[k, ],...)
		
		W <- xlift$W
		Gpre = W %*% t(Conj(W))	

            rem <- xlift$removelist	# same as trajmat[k,]
            coeff <- xlift$coeffv
	if (normalise) {
            coeff <- coeff/sqrt(diag(Gpre))
        }

            scales <- xlift$lengthsremove
            scalesa <- rep(NA, n)
            scalesa[rem] <- scales
            al <- artificial.levels(scales, rem, grid, tail = tail, 
                type = altype)
            levno <- length(al)
            energies <- energiesu <- NULL
            scalesx <- NULL
            l2vec <- NULL
            levvec <- NULL
            nj <- sapply(al, length)

	# modify here for complex coeffs
            for (j in 1:levno) {
              	energies[j] <- efun(coeff[al[[j]]])
                scalesx[j] <- afun(scalesa[al[[j]]], j = j)
                if (jsc) {
                  l2vec <- c(l2vec, log2(scalesa[al[[j]]]))
                  levvec <- c(levvec, rep(scalesx[j], times = nj[j]))
                }
            }

            if (jsc) {
                jstarc <- lm(l2vec ~ log2(levvec))$coef[2]
            }
            if (bc) {
                gg <- sapply(nj, gj)
            }
            else {
                gg <- 0
            }
            if (vc) {
               	levw <- 1/sapply(nj, levvar)
            }
            else {
                levw <- rep(1, times = levno)
            }
            for (i in 1:length(cutoffs)) {
                if (cut.fine) {
                  index <- (1 + cutoffs[i]):levno
                }
                else {
                  index <- 1:(levno - cutoffs[i])
                }
                beta[i] <- lm((log2(energies) - gg)[index] ~ 
                  log2(scalesx[index]), weights = levw[index])$coef[2]
                if (jsc) {
                  beta[i] <- beta[i] * jstarc
                }
            }
            betamat <- rbind(betamat, beta)
        }
        cat("\n")
        beta <- apply(betamat, 2, mean)
        Hs <- Hfrombeta(betamat, model = model)
        sds <- apply(Hs, 2, sd)
        H <- apply(Hs, 2, mean)
        ci <- t(apply(Hs, 2, bootci, level = level))
        bandH <- cbind(beta, H, sds, ci)

    if (BHonly) {
        return(bandH)
    }
    else {
        return(list(bandH = bandH, energies = (log2(energies) - 
            gg)[index], scales = log2(scalesx[index])))
    }
}