#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)] }