wave <- function(data, length, N=2, j1=2, j2=8, plotflag = T)
#Function to do the Wavelet estimator of H.
#Based on work by Ardry and Flandrin.
#Written by Vadim Teverovsky 1997.
#Calls function wd, available in the statlib package on wavelets.
#also fucntion accessD, from same package.
#data	--	Data
#length	--	Length of data to be used (must be power of 2)	
#j1,j2	--	Beginning and ending octave for estimation.
#noctave --	Number of octaves available.
#bound.effect -- Do not use first and last few coefficients, due to bound. 
#			effects.
#
{
        transform <- wd(data[1:(length)], filter = N)
        noctave <- log(length, base = 2) - 1
        statistic <- rep(0, noctave)
        bound.effect <- ceiling(log(2 * N, base = 2))
        if(j2 > noctave - bound.effect) {
               cat("WARNING Upper bound too high, resetting j2 to ", noctave - 
                        bound.effect, "\n")
                j2 <- noctave - bound.effect
        }
        for(j in 1:(noctave - bound.effect)) {
                statistic[j] <- log(mean((accessD(transform, lev = (noctave + 
		    1 - j))[N:(2^(noctave + 1 - j) - N)])^2), base = 2) 
        }
        slope <- lsfit(c(j1:j2), statistic[j1:j2])$coef
        if(plotflag == T) {
                plot(c(1:(noctave - bound.effect)), statistic[1:(noctave - 
                        bound.effect)], xlab = "Octave", ylab = " ", pch = 3)
                abline(slope[1], slope[2])
        }
        result <- 0.5 * (slope[2] + 1)
        names(result) <- NULL
        return(result)
}
