# foreca.EM-aux: ForeCA EM auxiliary functions In ForeCA: Forecastable Component Analysis

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

• 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) ≥q 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.

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