ForeCA (version 0.2.4)

foreca.EM-aux: ForeCA EM auxiliary functions

Description

foreca.EM.one_weightvector relies on several auxiliary functions:

foreca.EM.E_step computes the spectral density of \(y_t=\mathbf{U}_t \mathbf{w}\) given the weightvector \(\mathbf{w}\) and the normalized spectrum estimate \(f_{\mathbf{U}}\). A wrapper around spectrum_of_linear_combination.

foreca.EM.M_step computes the minimizing eigenvector (\(\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 \(\widehat{\mathbf{w}}_i\).

foreca.EM.E_and_M_step is a wrapper around foreca.EM.E_step followed by foreca.EM.M_step.

foreca.EM.h evaluates (an upper bound of) the entropy of the spectral density as a function of \(\mathbf{w}_i\) (or \(\mathbf{w}_{i+1}\)). This is the objective funcion that should be minimized.

Usage

foreca.EM.E_step(f.U, weightvector)
foreca.EM.M_step(f.U, f.current, minimize = TRUE, entropy.control = list())
foreca.EM.E_and_M_step(weightvector, f.U, minimize = TRUE,
  entropy.control = list())
foreca.EM.h(weightvector.new, f.U, weightvector.current = weightvector.new,
  f.current = NULL, entropy.control = list(), return.negative = FALSE)

Arguments

f.U

multivariate spectrum of class 'mvspectrum' with normalize = TRUE.

weightvector

numeric; weights \(\mathbf{w}\) for \(y_t = \mathbf{U}_t \mathbf{w}\). Must have unit norm in \(\ell^2\).

f.current

numeric; spectral density estimate of \(y_t=\mathbf{U}_t \mathbf{w}\) for the current estimate \(\widehat{\mathbf{w}}_i\) (required for foreca.EM.M_step; optional for foreca.EM.h).

minimize

logical; if TRUE (default) it returns the eigenvector corresponding to the smallest eigenvalue; otherwise to the largest eigenvalue.

entropy.control

list; control settings for entropy estimation. See complete_entropy_control for details.

weightvector.new

weightvector \(\widehat{\mathbf{w}}_{i+1}\) of the new iteration (i+1).

weightvector.current

weightvector \(\widehat{\mathbf{w}}_{i}\) of the current iteration (i).

return.negative

logical; if 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: FALSE.

Value

foreca.EM.E_step returns the normalized univariate spectral density (normalized such that its sum equals \(0.5\)).

foreca.EM.M_step returns a list with three elements:

• matrix: weighted covariance matrix, where the weights are the negative log of the spectral density. If density is estimated by discrete probabilities, then this matrix is positive semi-definite, since \(-\log(p) \geq 0\) for \(p \in [0, 1]\). See weightvector2entropy_wcov.

• vector: minimizing (or maximizing if minimize = FALSE) eigenvector of matrix,

• value: corresponding eigenvalue.

Contrary to foreca.EM.M_step, foreca.EM.E_and_M_step only returns the optimal weightvector as a numeric.

foreca.EM.h returns non-negative real value (see References for details):

• entropy, if weightvector.new = weightvector.current,

• an upper bound of that entropy for weightvector.new, otherwise.

See Also

weightvector2entropy_wcov

Examples

# NOT RUN {
XX <- diff(log(EuStockMarkets)) * 100
UU <- whiten(XX)$U
ff <- mvspectrum(UU, 'wosa', normalize = TRUE)

ww0 <- initialize_weightvector(num.series = ncol(XX), method = 'rnorm')

f.ww0 <- foreca.EM.E_step(ff, ww0)
plot(f.ww0, type = "l")

one.step <- foreca.EM.M_step(ff, f.ww0, 
                             entropy.control = list(prior.weight = 0.1))
image(one.step$matrix)
# }
# NOT RUN {
requireNamespace(LICORS)
# if you have the 'LICORS' package use
LICORS::image2(one.step$matrix)
# }
# NOT RUN {
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

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.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)

# }