#These are the functions performing the local Whittle algorithm.
#	rfunc	--	Function which is minimized in the local Whittle
#			estimator. 
#	locwhitt  --	Basic local whittle.  Finds periodogram (uses
#			function   per  ), then minimizes rfunc.
#	locwhitt.cont -- Continuous version of local Whittle.  Calculates
#			local Whittle with various values of m, and returns
#			a vector of results. 
#	locwhitt.cont.theory  -- Takes as argument the theoretical spectral
#			density of a time series model instead of an actual
#			time series, and produces the same vector as in
#			locwhitt.cont.  Useful for predicting the BIAS of
#			this estimator for known time series models. 
#
# These are functions which perform an algorithm very similar to Local
# Whittle, but using the cumulative periodogram instead of the periodogram,
# and using the appropriate minimizing function.
# 
#	rfunc2	--	Function which is minimized.
#	locwhitt2  --	Local whittle.  Finds cumulative periodogram (uses
#			function   per  ), then minimizes rfunc2.
#	locwhitt2.cont -- Continuous version of locwhitt2.	
#	locwhitt2.cont.theory  -- Takes as argument the theoretical spectral
#			density of a time series model instead of an actual
#			time series, and produces the same vector as in
#			locwhitt2.cont.  Useful for predicting the BIAS of
#			this estimator for known time series models. 


rfunc_function(h, len, im, peri)
#	h	-- Starting H value for minimization.
#	len	-- Length of time series.
#	im	-- Use only len/im frequencies.
#	peri	-- Periodogram of data.
{
        m <- len %/% im
        peri <- peri[2:(m + 1)]
        z <- c(1:m)
        freq <- (2 * pi)/len * z
        result <- log(sum(freq^(2 * h - 1) * peri)) - (2 * h)/m * sum(log(freq)
                )       #       cat("H = ", h, "R = ", result, "\n")
        drop(result)
}


locwhitt_function(data, h = 0.5, im = 2)
#	data	--	Time series.
#	h	-- Starting H value for minimization.
#	im	-- Use only N/im frequencies where N is length of series.

{
        peri <- per(data)
        len <- length(data)
        return(nlminb(start = h, obj = rfunc, len = len, im = im, peri = peri)$
                parameters)
}

locwhitt.cont_function(data, h = 0.5, i = 20, minm = 50)
#	data	--	Time series.
#	h	-- Starting H value for minimization.
#	i	-- Number of different m's to use.
#	minm	-- Minimum number of frequencies to use in estimation.

{
        peri <- per(data)
        len <- length(data)
        len1 <- len/minm
        cat("len = ", len, "len1 = ", len1, "\n")
        hvec <- rep(0, i)
        if(i > len1)
                i <- len1
        for(j in 1:i) {
                im <- as.integer(len1/i * j)
                hvec[j] <- nlminb(start = h, obj = rfunc, len = len, im = im, 
                        peri = peri)$parameters
                cat("im = ", im, "\n\n")
        }
        return(hvec)
}



locwhitt.cont.theory_function(specdens, h = 0.5, i = 20, minm = 50)
#	specdens -- Spectral density of a known time series model.
#	h	-- Starting H value for minimization.
#	i	-- Number of different m's to use.
#	minm	-- Minimum number of frequencies to use in estimation.

{
        peri <- c(0, specdens)
        len <- 2 * length(specdens)
        len1 <- len/minm
        cat("len = ", len, "len1 = ", len1, "\n")
        hvec <- rep(0, i)
        if(i > len1)
                i <- len1
        for(j in 1:i) {
                im <- as.integer(len1/i * j)
                hvec[j] <- nlminb(start = h, obj = rfunc, len = len, im = im, 
                        peri = peri)$parameters
                cat("im = ", im, "\n\n")
        }
        return(hvec)
}

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

rfunc2_function(h, len, il, m, peri)
#	h	-- Starting H value for minimization.
#	len	-- Length of time series.
#	il	-- Start after il frequencies.
#	im	-- Use only len/im frequencies.
#	peri	-- Periodogram of data.
{
        z <- c((il + 1):m)
        freq <- (2 * pi)/len * z
        freq1 <- (2 * pi)/len * il
        res1 <- 1/(m - il) * sum(log10(freq^(2 - 2 * h) - freq1^(2 - 2 * h)))
        result <- log10(sum(peri/(freq^(2 - 2 * h) - freq1^(2 - 2 * h)))) + 
                res1
        drop(result)
}


locwhitt2_function(data, h = 0.5, im = 2, il = 1)
#	data	--	Time series.
#	h	-- Starting H value for minimization.
#	len	-- Length of time series.
#	il	-- Start after il frequencies.
#	im	-- Use only len/im frequencies.
{
        peri.1 <- per(data)
        len <- length(data)
        m <- len %/% im
        peri <- cumsum(peri.1[(il + 2):(m + 1)])
        return(nlminb(start = h, obj = rfunc2, lower = 0.0001, upper = 0.9999, 
                len = len, il = il, m = m, peri = peri)$parameters)
}


locwhitt2.cont_function(data, h = 0.5, i = 20, minm = 50)
#	data	--	Time series.
#	h	-- Starting H value for minimization.
#	i	-- Number of different m's to use.
#	minm	-- Minimum number of frequencies to use in estimation.
{
        peri1 <- per(data)
        peri <- cumsum(peri1)
        len <- length(data)
        len1 <- len/minm
        cat("len = ", len, "len1 = ", len1, "\n")
        hvec <- rep(0, i)
        if(i > len1)
                i <- len1
        for(j in 1:i) {
                im <- as.integer(len1/i * j)
                hvec[j] <- nlminb(start = h, obj = rfunc2, len = len, im = im, 
                        peri = peri)$parameters
                cat("im = ", im, "\n\n")
        }
        return(hvec)
}


locwhitt2.cont.theory_function(specdens, h = 0.5, i = 20, minm = 50)
#	specdens -- Spectral density of a known time series model.
#	h	-- Starting H value for minimization.
#	i	-- Number of different m's to use.
#	minm	-- Minimum number of frequencies to use in estimation.
{
        peri <- c(0, cumsum(specdens))
        len <- 2 * length(specdens)
        len1 <- len/minm
        cat("len = ", len, "len1 = ", len1, "\n")
        hvec <- rep(0, i)
        if(i > len1)
                i <- len1
        for(j in 1:i) {
                im <- as.integer(len1/i * j)
                hvec[j] <- nlminb(start = h, obj = rfunc2, len = len, im = im, 
                        peri = peri)$parameters
                cat("im = ", im, "\n\n")
        }
        return(hvec)
}


#Function to compute periodogram.

per_function(z)
{
        n <- length(z)
        (Mod(fft(z))^2/(2 * pi * n))[1:(n %/% 2 + 1)]
}
