foreca.EM-aux: ForeCA EM auxiliary functions

Description Usage Arguments Value See Also Examples

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

1
2
3
4
5
6
7
8
9
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:

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

See Also

weightvector2entropy_wcov

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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)

## End(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)

ForeCA documentation built on May 29, 2017, 9:09 a.m.