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

#### Documented in initialize_weightvector

#' @title Initialize weightvector for iterative ForeCA algorithms
#' @description
#' \code{initialize_weightvector} returns a unit norm (in \eqn{\ell^2})
#' vector \eqn{\mathbf{w}_0 \in R^K} that can be used as the starting
#' point for any iterative ForeCA algorithm, e.g.,
#' quickly computable heuristics are available via the \code{method} argument.
#'
#' @keywords manip
#' @inheritParams common-arguments
#' @inheritParams sfa
#' @param num.series positive integer; number of time series \eqn{K} (determines the length
#' of the weightvector). If \code{num.series = 1} it simply returns
#' a 1 \eqn{\times} 1 array equal to \code{1}.
#' @param method string; which heuristics should be used to generate a good starting \eqn{\mathbf{w}_0}?
#' Default: \code{"rnorm"}; see Details.
#' @param seed non-negative integer; seed for random initialization which will be
#' returned for reproducibility. By default it sets a random seed.
#' @details
#' The \code{method} argument specifies the heuristics that is used to get a good
#' starting vector \eqn{\mathbf{w}_0}:
#'
#' \itemize{
#'  \item{\code{"max"}}{ vector with all \eqn{0}s, but a \eqn{1} at the position
#'  of the maximum forecastable series in \code{U}.}
#'  \item{\code{"rcauchy"}}{ random start using \code{rcauchy(k)}.}
#'  \item{\code{"rnorm"}}{ random start using \code{rnorm(k, 0, 1)}.}
#'  \item{\code{"runif"}}{ random start using \code{runif(k, -1, 1)}.}
#'  \item{\code{"PCA.large"}}{ first eigenvector of PCA (largest variance signal).}
#'  \item{\code{"PCA.small"}}{ last eigenvector of PCA (smallest variance signal).}
#'  \item{\code{"PCA"}}{ checks both small and large, and chooses the one with higher
#'             forecastability as computed by \code{\link{Omega}}..}
#'  \item{\code{"SFA.fast"}}{ last eigenvector of SFA (fastest signal).}
#'  \item{\code{"SFA.slow"}}{ first eigenvector of SFA (slowest signal).}
#'  \item{\code{"SFA"}}{ checks both slow and fast, and chooses the one with higher
#'             forecastability as computed by \code{\link{Omega}}.}
#' }
#'
#' Each vector has length K and is automatically normalized to have unit norm
#' in \eqn{\ell^2}.
#'
#' For the \code{'SFA*'} methods see \code{\link{sfa}}.
#' Note that maximizing (or minimizing) the lag \eqn{1} auto-correlation does
#' not necessarily yield the most forecastable signal, but it's a good start.
#' @return
#' numeric; a vector of length \eqn{K} with unit norm in \eqn{\ell^2}.
#' @examples
#' XX <- diff(log(EuStockMarkets))
#' \dontrun{
#' initialize_weightvector(U = XX, method = "SFA")
#' }
#' initialize_weightvector(num.series = ncol(XX), method = "rnorm")
#' @export
#'
initialize_weightvector <- function(U = NULL, f.U = NULL,
num.series = ncol(U),
method = c("rnorm", "max", "SFA", "PCA",
"rcauchy", "runif", "SFA.slow", "SFA.fast",
"PCA.large", "PCA.small"),
seed = sample(1e6, 1), ...) {

stopifnot(is.null(U) || is.array(U) || is.ts(U) || is.matrix(U),
is.null(f.U) || is.array(f.U),
num.series >= 1,
is.character(method))

if (!is.null(U)) {
# if U is provided it must be whitened
U <- check_whitened(U)
}

method <- match.arg(method)

if (is.null(num.series) && is.null(f.U)) {
stop("You must provide either 'num.series', 'series', or 'f.U'.")
}

if (method %in% c("SFA", "SFA.slow", "SFA.fast",
"PCA", "PCA.large", "PCA.small") && is.null(U)) {
stop("For SFA- or PCA-type methods you must provide data via the 'U' argument.")
}

if (is.null(num.series)) {
if (is.null(f.U)) {
num.series <- ncol(U)
} else if (is.null(U)) {
num.series <- dim(f.U)[2]
} else {
stop("Something went wrong in initialize_weightvector().")
}
}

if (num.series == 1) {
return(cbind(1))
}

set.seed(seed)
if (method == "rcauchy") {
ww0 <- rcauchy(num.series)
} else if (method == "runif") {
ww0 <- runif(num.series, -1, 1)
} else if (method == "rnorm") {
ww0 <- rnorm(num.series)
} else if (method == "max") {
ww0 <- rep(0, num.series)
Omega.tmp <- apply(get_spectrum_from_mvspectrum(f.U), 2,
function(x) Omega(mvspectrum.output = x))
ww0[which.max(Omega.tmp)] <- 1
} else if (any(method == c("SFA", "SFA.slow", "SFA.fast"))) {
sfa.est <- ForeCA::sfa(U, ...)
ww.slow <- sfa.est$loadings[, 1] ww.fast <- sfa.est$loadings[, num.series]

if (method == "SFA.slow") {
ww0 <- ww.slow
} else if (method == "SFA.fast") {
ww0 <- ww.fast
} else if (method == "SFA") {
stopifnot(!is.null(f.U))

omega.slow <- Omega(mvspectrum.output = foreca.EM.E_step(f.U, ww.slow))
omega.fast <- Omega(mvspectrum.output = foreca.EM.E_step(f.U, ww.fast))
if (omega.fast > omega.slow) {
ww0 <- ww.fast
} else {
ww0 <- ww.slow
}
}
} else if (method %in% c("PCA", "PCA.large", "PCA.small")) {
pca.est <- princomp(U, ...)
ww.large <- pca.est$loadings[, 1] ww.small <- pca.est$loadings[, num.series]

if (method == "PCA.large") {
ww0 <- ww.large
} else if (method == "PCA.small") {
ww0 <- ww.small
} else if (method == "PCA") {
stopifnot(!is.null(f.U))

omega.large <- Omega(mvspectrum.output = foreca.EM.E_step(f.U, ww.large))
omega.small <- Omega(mvspectrum.output = foreca.EM.E_step(f.U, ww.small))
if (omega.small > omega.large) {
ww0 <- ww.small
} else {
ww0 <- ww.large
}
}
} else {
stop("Method '", method, "' is not available for initialize_weightvector().")
}

if (ww0[1] != 0) {
# make the first entry always positive (if it's not zero)
ww0 <- ww0 * sign(ww0[1])
}
# normalize
ww0 <- rbind(ww0 / base::norm(ww0, "2"))
rownames(ww0) <- NULL
return(ww0)
}


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