# R/spectral_entropy.R In ForeCA: Forecastable Component Analysis

#### Documented in spectral_entropy

#' @title Estimates spectral entropy of a time series
#'
#' @description
#' Estimates \emph{spectral entropy} from a univariate (or multivariate)
#' normalized spectral density.
#'
#' @details
#' The \emph{spectral entropy} equals the Shannon entropy of the spectral density
#' \eqn{f_x(\lambda)} of a stationary process \eqn{x_t}:
#' \deqn{
#' H_s(x_t) = - \int_{-\pi}^{\pi} f_x(\lambda) \log f_x(\lambda) d \lambda,
#' }
#' where the density is normalized such that
#' \eqn{\int_{-\pi}^{\pi} f_x(\lambda) d \lambda = 1}. An estimate of \eqn{f(\lambda)}
#'  can be obtained
#' by the (smoothed) periodogram (see \code{\link{mvspectrum}}); thus using discrete, and
#' not continuous entropy.
#'
#' @inheritParams common-arguments
#' @param series univariate time series of length \eqn{T}.  In the rare case
#' that users want to call this for a multivariate time \code{series}, note
#' that the estimated spectrum is in general \emph{not} normalized for the computation.
#' Only if the original data is whitened, then it is normalized.
#' @param mvspectrum.output optional; one can directly provide an estimate of
#' the spectrum of \code{series}. Usually the output of \code{\link{mvspectrum}}.
#' @return
#' A non-negative real value for the spectral entropy \eqn{H_s(x_t)}.
#' @references
#' Jerry D. Gibson and Jaewoo Jung (2006). \dQuote{The
#' Interpretation of Spectral Entropy Based Upon Rate Distortion Functions}.
#' IEEE International Symposium on Information Theory, pp. 277-281.
#'
#' L. L. Campbell, \dQuote{Minimum coefficient rate for stationary random processes},
#' Information and Control, vol. 3, no. 4, pp. 360 - 371, 1960.
#'
#' @keywords ts univar math
#' @export
#' @examples
#'
#' set.seed(1)
#' eps <- rnorm(100)
#' spectral_entropy(eps)
#'
#' phi.v <- seq(-0.95, 0.95, by = 0.1)
#' kMethods <- c("mvspec", "pgram")
#' SE <- matrix(NA, ncol = length(kMethods), nrow = length(phi.v))
#' for (ii in seq_along(phi.v)) {
#'   xx.tmp <- arima.sim(n = 200, list(ar = phi.v[ii]))
#'   for (mm in seq_along(kMethods)) {
#'     SE[ii, mm] <- spectral_entropy(xx.tmp, spectrum.control =
#'                                      list(method = kMethods[mm]))
#'   }
#' }
#'
#' matplot(phi.v, SE, type = "l", col = seq_along(kMethods))
#' legend("bottom", kMethods, lty = seq_along(kMethods),
#'        col = seq_along(kMethods))
#'
#' # AR vs MA
#' SE.arma <- matrix(NA, ncol = 2, nrow = length(phi.v))
#' SE.arma[, 1] <- SE[, 2]
#'
#' for (ii in seq_along(phi.v)){
#'   yy.temp <- arima.sim(n = 1000, list(ma = phi.v[ii]))
#'   SE.arma[ii, 2] <-
#'     spectral_entropy(yy.temp, spectrum.control = list(method = "mvspec"))
#' }
#'
#' matplot(phi.v, SE.arma, type = "l", col = 1:2, xlab = "parameter (phi or theta)",
#'         ylab = "Spectral entropy")
#' abline(v = 0, col = "blue", lty = 3)
#' legend("bottom", c("AR(1)", "MA(1)"), lty = 1:2, col = 1:2)
#'

spectral_entropy <- function(series = NULL, spectrum.control = list(),
entropy.control = list(),
mvspectrum.output = NULL, ...){

stopifnot(xor(is.null(series), is.null(mvspectrum.output)),
is.null(mvspectrum.output) ||
length(dim(mvspectrum.output)) >= 0)

if (is.null(mvspectrum.output)) {
series <- as.matrix(series)
is.whitened <- (isTRUE(all.equal(cov(series), diag(1, ncol(series)))) &&
isTRUE(all.equal(colMeans(series), 0)))
attr(series, "whitened") <- is.whitened
num.series <- ncol(series)
spectrum.control <- complete_spectrum_control(spectrum.control)
mvspectrum.output <- mvspectrum(series,
method = spectrum.control$method, smoothing = spectrum.control$smoothing,
normalize = attr(series, "whitened"),
...)
}
if (is.null(dim(mvspectrum.output))) {
dim(mvspectrum.output) <- c(length(mvspectrum.output), 1, 1)
}
num.series <- dim(mvspectrum.output)[2]

if (num.series > 1) {
num.freqs <- dim(mvspectrum.output)[1]
} else {
num.freqs <- length(mvspectrum.output)
}

entropy.control$base <- NULL entropy.control <- complete_entropy_control(entropy.control, num.outcomes = num.freqs * 2) if (num.series > 1) { spec.ent <- mvspectrum2wcov(mvspectrum.output * -log(mvspectrum.output, base = entropy.control$base))
} else {
# this has sum() = 1
spec.dens <- c(rev(c(mvspectrum.output)), c(mvspectrum.output))
spec.dens[spec.dens < 0] <- 0
spec.dens <- spec.dens / sum(spec.dens)
entropy.control <- complete_entropy_control(entropy.control,
num.outcomes = length(spec.dens))
entropy.control$prior.probs <- NULL spec.ent <- do.call("discrete_entropy", c(list(probs = spec.dens), entropy.control)) } attr(spec.ent, "nfrequencies") <- num.freqs attr(spec.ent, "base") <- entropy.control$base
return(spec.ent)
}


## Try the ForeCA package in your browser

Any scripts or data that you put into this service are public.

ForeCA documentation built on July 1, 2020, 7:50 p.m.