# R/foreca.EM-aux.R In ForeCA: Forecastable Component Analysis

#### Documented in foreca.EM.E_and_M_stepforeca.EM.E_stepforeca.EM.hforeca.EM.M_step

#' @title ForeCA EM auxiliary functions
#' @name foreca.EM-aux
#' @description
#' \code{foreca.EM.one_weightvector} relies on several auxiliary functions:
#'
NULL

#' @rdname foreca.EM-aux
#' @description
#' \code{foreca.EM.E_step} computes the spectral density of
#' \eqn{y_t=\mathbf{U}_t \mathbf{w}} given the weightvector \eqn{\mathbf{w}}
#' and the normalized spectrum estimate \eqn{f_{\mathbf{U}}}.
#' A wrapper around \code{\link{spectrum_of_linear_combination}}.
#' @keywords manip math
#'
#' @inheritParams common-arguments
#' @param weightvector numeric; weights \eqn{\mathbf{w}} for
#' \eqn{y_t = \mathbf{U}_t \mathbf{w}}. Must have unit norm in \eqn{\ell^2}.
#' @export
#' @return
#' \code{foreca.EM.E_step} returns the normalized univariate spectral
#' density (normalized such that its \code{sum} equals \eqn{0.5}).
#' @examples
#' \dontrun{
#' XX <- diff(log(EuStockMarkets)) * 100
#' UU <- whiten(XX)$U #' ff <- mvspectrum(UU, 'mvspec', normalize = TRUE) #' #' ww0 <- initialize_weightvector(num.series = ncol(XX), method = 'rnorm') #' #' f.ww0 <- foreca.EM.E_step(ff, ww0) #' plot(f.ww0, type = "l") #' } foreca.EM.E_step <- function(f.U, weightvector) { stopifnot(!any(is.na(weightvector)), isTRUE(all.equal(target = 1, current = base::norm(weightvector, "2"))), length(dim(f.U)) == 3, dim(f.U)[2] == dim(f.U)[3], # square dim(f.U)[2] == length(weightvector)) # weightvector has the right dimension check_mvspectrum_normalized(f.U) spec.dens.est <- spectrum_of_linear_combination(f.U, weightvector) check_mvspectrum_normalized(spec.dens.est) # is also a test already return(spec.dens.est) } #' @rdname foreca.EM-aux #' @description #' \code{foreca.EM.M_step} computes the minimizing eigenvector #' (\eqn{\rightarrow \widehat{\mathbf{w}}_{i+1}}) of the weighted #' covariance matrix, where the weights equal the negative logarithm of the #' spectral density at the current \eqn{\widehat{\mathbf{w}}_i}. #' @keywords manip #' @inheritParams common-arguments #' @param minimize logical; if \code{TRUE} (default) it returns the eigenvector #' corresponding to the smallest eigenvalue; otherwise to the largest eigenvalue. #' @inheritParams common-arguments #' @export #' @return #' \code{foreca.EM.M_step} returns a list with three elements: #' \itemize{ #' \item \code{matrix}: weighted covariance matrix, where the weights are #' the negative log of the spectral density. If density is estimated #' by discrete probabilities, #' then this \code{matrix} is positive semi-definite, since #' \eqn{-\log(p) \geq 0} for \eqn{p \in [0, 1]}. #' See \code{\link{weightvector2entropy_wcov}}. #' \item \code{vector}: minimizing (or maximizing if #' \code{minimize = FALSE}) eigenvector of \code{matrix}, #' \item \code{value}: corresponding eigenvalue. #' } #' @examples #' \dontrun{ #' one.step <- foreca.EM.M_step(ff, f.ww0, #' entropy.control = list(prior.weight = 0.1)) #' image(one.step$matrix)
#'
#' requireNamespace(LICORS)
#' # if you have the 'LICORS' package use
#' LICORS::image2(one.step$matrix) #' #' ww1 <- one.step$vector
#' f.ww1 <- foreca.EM.E_step(ff, ww1)
#'
#' layout(matrix(1:2, ncol = 2))
#' matplot(seq(0, pi, length = length(f.ww0)), cbind(f.ww0, f.ww1),
#'         type = "l", lwd =2, xlab = "omega_j", ylab = "f(omega_j)")
#' plot(f.ww0, f.ww1, pch = ".", cex = 3, xlab = "iteration 0",
#'      ylab = "iteration 1", main = "Spectral density")
#' abline(0, 1, col = 'blue', lty = 2, lwd = 2)
#'
#' Omega(mvspectrum.output = f.ww0) # start
#' Omega(mvspectrum.output = f.ww1) # improved after one iteration
#' }

foreca.EM.M_step <- function(f.U, f.current, minimize = TRUE,
entropy.control = list()) {
stopifnot(all(f.current >= 0),
is.logical(minimize))

num.freqs <- length(f.current)

check_mvspectrum_normalized(f.U)
check_mvspectrum_normalized(f.current)

entropy.control <- complete_entropy_control(entropy.control,
num.outcomes = 2 * num.freqs)
integrated.spectrum.entropy <-
weightvector2entropy_wcov(NULL,
f.U = f.U,
f.current = f.current,
entropy.control = entropy.control)
EE <- eigen(integrated.spectrum.entropy, symmetric = TRUE)
sel <- ifelse(minimize, which.min(EE$values), which.max(EE$values))

weightvector <- EE$vector[, sel] # make first entry always positive (scale vector always this way) # for consistent results among re-runs weightvector <- weightvector / sign(weightvector[1]) out <- list(matrix = integrated.spectrum.entropy, vector = weightvector, value = EE$values[sel])
return(out)
}

#' @rdname foreca.EM-aux
#' @description
#' \code{foreca.EM.E_and_M_step} is a wrapper around \code{foreca.EM.E_step}
#' followed by \code{foreca.EM.M_step}.
#' @keywords manip
#' @inheritParams common-arguments
#' @export
#' @return
#' Contrary to \code{foreca.EM.M_step}, \code{foreca.EM.E_and_M_step} only returns the optimal
#' weightvector as a numeric.
#' @examples
#' \dontrun{
#' ww0 <- initialize_weightvector(NULL, ff, method = "rnorm")
#' ww1 <- foreca.EM.E_and_M_step(ww0, ff)
#' ww0
#' ww1
#' barplot(rbind(ww0, ww1), beside = TRUE)
#' abline(h = 0, col = "blue", lty = 2)
#' }

foreca.EM.E_and_M_step <- function(weightvector, f.U, minimize = TRUE,
entropy.control = list()) {

f.current <- spectrum_of_linear_combination(f.U, weightvector)
opt.vec <- foreca.EM.M_step(f.U, f.current, minimize = minimize,
entropy.control = entropy.control)$vector return(opt.vec) } #' @rdname foreca.EM-aux #' @description #' \code{foreca.EM.h} evaluates (an upper bound of) the entropy of the spectral density as a function #' of \eqn{\mathbf{w}_i} (or \eqn{\mathbf{w}_{i+1}}). This is the objective funcion that should be #' \code{minimize}d. #' #' @keywords manip #' @param weightvector.new weightvector \eqn{\widehat{\mathbf{w}}_{i+1}} of the new #' iteration (i+1). #' @param f.current numeric; spectral density estimate of #' \eqn{y_t=\mathbf{U}_t \mathbf{w}} for the current estimate #' \eqn{\widehat{\mathbf{w}}_i} (required for #' \code{foreca.EM.M_step}; optional for \code{foreca.EM.h}). #' @param weightvector.current weightvector \eqn{\widehat{\mathbf{w}}_{i}} of the #' current iteration (i). #' @param return.negative logical; if \code{TRUE} it returns the negative #' spectral entropy. This is useful when maximizing forecastibility which is #' equivalent (up to an additive constant) to maximizing negative entropy. #' Default: \code{FALSE}. #' @return #' \code{foreca.EM.h} returns non-negative real value (see References for details): #' \itemize{ #' \item entropy, if \code{weightvector.new = weightvector.current}, #' \item an upper bound of that entropy for \code{weightvector.new}, #' otherwise. #' } #' @seealso #' \code{\link{weightvector2entropy_wcov}} #' @export #' @examples #' \dontrun{ #' foreca.EM.h(ww0, ff) # iteration 0 #' foreca.EM.h(ww1, ff, ww0) # min eigenvalue inequality #' foreca.EM.h(ww1, ff) # KL divergence inequality #' one.step$value
#'
#' # by definition of Omega, they should equal 1 (modulo rounding errors)
#' Omega(mvspectrum.output = f.ww0) / 100 + foreca.EM.h(ww0, ff)
#' Omega(mvspectrum.output = f.ww1) / 100 + foreca.EM.h(ww1, ff)
#' }

foreca.EM.h <- function(weightvector.new, f.U,
weightvector.current = weightvector.new,
f.current = NULL,
entropy.control = list(),
return.negative = FALSE) {

# short as quadratic_form(apply(f.U*-log(weightvector.current), 2:3, sum),
# weightvector.new)
stopifnot(length(weightvector.current) == length(weightvector.new),
all(f.current >= 0) || is.null(f.current))

check_mvspectrum_normalized(f.U)
if (is.null(f.current)) {
f.current <- foreca.EM.E_step(f.U, weightvector.current)
} else {
check_mvspectrum_normalized(f.current)
}

num.series <- length(weightvector.new)
num.freqs <- length(f.current)

entropy.control <- complete_entropy_control(entropy.control,
num.outcomes = 2 * num.freqs)
integrated.spectrum.entropy <-
weightvector2entropy_wcov(NULL,
f.U = f.U,
f.current = f.current,
entropy.control = entropy.control)