From acb845ab8d2df2a499ffde3f35fd096ced707e6f Mon Sep 17 00:00:00 2001 From: Khromov Nikita Date: Thu, 30 Oct 2025 15:37:48 +0300 Subject: [PATCH 1/2] Implement complex multichannel ssa --- NAMESPACE | 4 + R/cmhankel.R | 332 +++++++++++++++++++++++++++++++++++++++++++++++++++ R/ssa.R | 16 ++- 3 files changed, 347 insertions(+), 5 deletions(-) create mode 100644 R/cmhankel.R diff --git a/NAMESPACE b/NAMESPACE index fbe7b8ea..e909ff92 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -27,6 +27,7 @@ export(clone, new.hmat, hmatmul, hankel, + mhankel, hcols, hrows, is.hmat, @@ -73,6 +74,7 @@ S3method("clone", ssa) S3method("decompose", "ssa") S3method("decompose", "toeplitz.ssa") S3method("decompose", "cssa") +S3method("decompose", "cmssa") S3method("decompose", "pssa") S3method("decompose", "ossa") S3method("reconstruct", ssa) @@ -80,6 +82,7 @@ S3method("residuals", ssa) S3method("residuals", "ssa.reconstruction") S3method("calc.v", "ssa") S3method("calc.v", "cssa") +S3method("calc.v", "cmssa") S3method("$", ssa) S3method("print", ssa) S3method("print", ossa) @@ -92,6 +95,7 @@ S3method("plot", "2d.ssa.reconstruction") S3method("plot", "nd.ssa.reconstruction") S3method("plot", "mssa.reconstruction") S3method("plot", "cssa.reconstruction") +S3method("plot", "cmssa.reconstruction") S3method("plot", "grouping.auto.wcor") S3method("plot", "grouping.auto.pgram") S3method("wcor", "default") diff --git a/R/cmhankel.R b/R/cmhankel.R new file mode 100644 index 00000000..3068bff0 --- /dev/null +++ b/R/cmhankel.R @@ -0,0 +1,332 @@ +# R package for Singular Spectrum Analysis +# Copyright (c) 2013 Anton Korobeynikov +# +# This program is free software; you can redistribute it +# and/or modify it under the terms of the GNU General Public +# License as published by the Free Software Foundation; +# either version 2 of the License, or (at your option) +# any later version. +# +# This program is distributed in the hope that it will be +# useful, but WITHOUT ANY WARRANTY; without even the implied +# warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +# PURPOSE. See the GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public +# License along with this program; if not, write to the +# Free Software Foundation, Inc., 675 Mass Ave, Cambridge, +# MA 02139, USA. + +.traj.dim.cmssa <- function(x) { + Ldim <- sum(x$wmask) + if (Ldim == 0) + Ldim <- x$window + + Kdim <- sum(x$fmask) + if (Kdim == 0) + Kdim <- sum(x$length - x$window + 1) + + c(Ldim, Kdim) +} + +.chmat.striped <- function(x, fft.plan, re = TRUE) { + N <- x$length; L <- x$window + + if (re) + F <- lapply(.F(x), Re) + else + F <- lapply(.F(x), Im) + + field <- matrix(0., max(N), length(N)) + + weights <- .get(x, "weights") + if (!is.null(weights)) + mask <- weights > 0 + else + mask <- matrix(TRUE, max(N), length(N)) + + for (idx in seq_along(N)) { + imask <- which(mask[seq_len(N[idx]), idx]) + field[imask, idx] <- F[[idx]][imask] + } + + new.hbhmat(field, L = c(L, 1), + wmask = NULL, + fmask = .get(x, "fmask"), + weights = weights) +} + +.get.or.create.cmhmat <- function(x, re = TRUE) { + .get.or.create(x, paste0("hmat_", if (re) "r" else "i"), + .chmat.striped(x, fft.plan = .get.or.create.cfft.plan(x), re = re)) +} + +.get.or.create.trajmat.cmssa <- .get.or.create.cmhmat + +# TODO: add diagonal averaging similarly to the hankel function. +# TODO: efficient implementation similarly to new.hbhmat. +mhankel <- function(X, L) { + if (!inherits(X, "series.list")) + X <- .to.series.list(X) + + Reduce(cbind, lapply(X, hankel, L = L)) +} + +decompose.cmssa <- function(x, + neig = NULL, + ..., + force.continue = FALSE) { + # Check, whether continuation of decomposition is requested + if (!force.continue && nsigma(x) > 0) + stop("Continuation of decomposition is not supported for this method.") + + if (is.null(neig)) + neig <- .default.neig(x, ...) + + if (identical(x$svd.method, "svd")) { + S <- svd(mhankel(.F(x), L = x$window), nu = neig, nv = neig) + .set.decomposition(x, sigma = S$d, U = S$u, V = S$v) + } else if (identical(x$svd.method, "eigen")) { + h <- mhankel(.F(x), L = x$window) + + ## FIXME: Build the complex L-covariance matrix properly + S <- eigen(tcrossprod(h, Conj(h)), symmetric = TRUE) + + ## Fix small negative values + S$values[S$values < 0] <- 0 + + ## Save results + .set.decomposition(x, + sigma = sqrt(S$values[seq_len(neig)]), + U = S$vectors[, seq_len(neig), drop = FALSE]) + } else if (identical(x$svd.method, "primme")) { + if (!requireNamespace("PRIMME", quietly = TRUE)) + stop("PRIMME package is required for SVD method `primme'") + + R <- .get.or.create.cmhmat(x, re = TRUE) + I <- .get.or.create.cmhmat(x, re = FALSE) + + matmul <- function(x, y) { + if (is.matrix(y)) + apply(y, 2, ematmul, emat = x, transposed = FALSE) + else + ematmul(x, y, transposed = FALSE) + } + + tmatmul <- function(x, y) { + if (is.matrix(y)) + apply(y, 2, ematmul, emat = x, transposed = TRUE) + else + ematmul(x, y, transposed = TRUE) + } + + A <- function(x, trans) { + rX <- Re(x) + iX <- Im(x) + if (identical(trans, "c")) { + (tmatmul(R, rX) + tmatmul(I, iX)) + + 1i * (tmatmul(R, iX) - tmatmul(I, rX)) + } else { + (matmul(R, rX) - matmul(I, iX)) + + 1i * (matmul(R, iX) + matmul(I, rX)) + } + } + + S <- PRIMME::svds(A, NSvals = neig, m = nrow(R), n = ncol(R), isreal = FALSE, ...) + .set.decomposition(x, sigma = S$d, U = S$u, V = S$v) + } else if (identical(x$svd.method, "irlba")) { + if (!requireNamespace("irlba", quietly = TRUE)) + stop("irlba package is required for SVD method `irlba'") + + h <- mhankel(.F(x), L = x$window) + S <- irlba::irlba(h, nv = neig, ...) + .set.decomposition(x, sigma = S$d, U = S$u, V = S$v) + } else if (identical(x$svd.method, "rsvd")) { + if (!requireNamespace("irlba", quietly = TRUE)) + stop("irlba package is required for SVD method `rsvd'") + + h <- mhankel(.F(x), L = x$window) + S <- irlba::svdr(h, k = neig, ...) + .set.decomposition(x, sigma = S$d, U = S$u, V = S$v) + } else + stop("unsupported SVD method") + + x +} + +calc.v.cmssa <- function(x, idx, ...) { + nV <- nv(x) + + V <- matrix(NA_complex_, .traj.dim(x)[2], length(idx)) + idx.old <- idx[idx <= nV] + idx.new <- idx[idx > nV] + + if (length(idx.old) > 0) { + V[, idx <= nV] <- .V(x)[, idx.old] + } + + if (length(idx.new) > 0) { + sigma <- .sigma(x)[idx.new] + + if (any(sigma <= .Machine$double.eps)) { + sigma[sigma <= .Machine$double.eps] <- Inf + warning("some sigmas are equal to zero. The corresponding vectors will be zeroed") + } + + U <- .U(x)[, idx.new, drop = FALSE] + + + tmatmul <- function(x, y) { + if (is.matrix(y)) + apply(y, 2, ematmul, emat = x, transposed = TRUE) + else + ematmul(x, y, transposed = TRUE) + } + + Rh <- .get.or.create.cmhmat(x, re = TRUE) + Ih <- .get.or.create.cmhmat(x, re = FALSE) + + hcrossprod <- function(x) { + rx <- Re(x) + ix <- Im(x) + (tmatmul(Rh, rx) + tmatmul(Ih, ix)) + + 1i * (tmatmul(Rh, ix) - tmatmul(Ih, rx)) + } + + V[, idx > nV] <- apply(U, 2, hcrossprod) / sigma + + } + + invisible(V) +} + +.hankelize.one.cmssa <- function(x, U, V, fft.plan = .get.or.create.cfft.plan(x)) { + hankelize.part <- function(x, U, V) { + h <- .get.or.create.cmhmat(x) + storage.mode(U) <- storage.mode(V) <- "double" + F <- .Call("hbhankelize_one_fft", U, V, h@.xData) + + ## FIXME: This is ugly + N <- x$length + mN <- max(N) + cidx <- unlist(lapply(seq_along(N), function(idx) + seq_len(N[idx]) + mN * (idx - 1))) + + F[cidx] + } + + R1 <- hankelize.part(x, Re(U), Re(V)) + R2 <- hankelize.part(x, Im(U), Im(V)) + I1 <- hankelize.part(x, Re(U), Im(V)) + I2 <- hankelize.part(x, Im(U), Re(V)) + + (R1 + R2) + 1i * (-I1 + I2) +} + +.elseries.cmssa <- function(x, idx, ...) { + if (max(idx) > nsigma(x)) + stop("Too few eigentriples computed for this decomposition") + + N <- x$length + sigma <- .sigma(x) + U <- .U(x) + F <- .F(x) + + res <- complex(sum(N)) + for (i in idx) { + if (nv(x) >= i) { + # FIXME: Check, whether we have factor vectors for reconstruction + # FIXME: Get rid of .get call + V <- .V(x)[, i] + } else { + # No factor vectors available. Calculate them on-fly. + V <- calc.v(x, i) + } + + res <- res + sigma[i] * .hankelize.one(x, U = U[, i], V = V) + } + + cN <- c(0, cumsum(N)) + sres <- list() + for (i in seq_along(N)) { + sres[[i]] <- res[(cN[i]+1):cN[i+1]] + attr(sres[[i]], "na.action") <- attr(F[[i]], "na.action") + } + class(sres) <- "series.list" + + sres +} + +plot.cmssa.reconstruction <- function(x, ...) { + rx = lapply(x, Re) + ix = lapply(x, Im) + + a <- attributes(x) + ar <- a + ar$series <- Re(a$series) + ar$residuals <- Re(a$residuals) + attributes(rx) <- ar + ai <- a + ai$series <- Im(a$series) + ai$residuals <- Im(a$residuals) + attributes(ix) <- ai + + dots <- list(...) + if (is.null(dots$main)) + title <- "Reconstructed Series" + else + title <- dots$main + dots$main <- NULL + + oldpar <- par(mfrow = c(1, 2)) + do.call(plot.mssa.reconstruction, c(list(x = rx, main = paste(title, "(Real part)")), dots)) + do.call(plot.mssa.reconstruction, c(list(x = ix, main = paste(title, "(Imaginary part)")), dots)) + par(oldpar) +} + +.init.fragment.cmssa <- function(this) + expression({ + if (any(circular)) + stop("Circular variant of complex multichannel SSA isn't implemented yet") + + # We assume that we have mts-like object. With series in the columns. + # Coerce input to series.list object + # Note that this will correctly remove leading and trailing NA's + x <- .to.series.list(x, na.rm = TRUE) + + # Sanity check - the input series should be complex + if (!all(sapply(x, is.complex))) + stop("complex SSA should be performed on complex time series") + + # Grab the inner attributes, if any + iattr <- lapply(x, attributes) + + N <- sapply(x, length) + + # If L is provided it should be length 1 + if (missing(L)) { + L <- (min(N) + 1) %/% 2 + } else { + if (length(L) > 1) + warning("length of L is > 1, only the first element will be used") + L <- L[1] + } + + wmask <- NULL + if (!all(N == max(N)) || any(sapply(x, anyNA))) { + K <- N - L + 1 + + weights <- matrix(0, max(N), length(N)) + fmask <- matrix(FALSE, max(K), length(N)) + wmask <- rep(TRUE, L) + for (idx in seq_along(N)) { + mask <- !is.na(x[[idx]]) + fmask[seq_len(K[idx]), idx] <- .factor.mask.1d(mask, wmask) + weights[seq_len(N[idx]), idx] <- .field.weights.1d(wmask, fmask[seq_len(K[idx]), idx]) + } + } else { + fmask <- weights <- NULL + } + + column.projector <- row.projector <- NULL + }) diff --git a/R/ssa.R b/R/ssa.R index 759959bc..8f2aab74 100644 --- a/R/ssa.R +++ b/R/ssa.R @@ -28,7 +28,7 @@ .determine.svd.method <- function(x, kind, neig = NULL, ..., - svd.method = (if (identical(kind, "cssa")) "eigen" else "nutrlan")) { + svd.method = (if ("cssa" %in% kind) "eigen" else "nutrlan")) { tjdim <- .traj.dim(x) L <- tjdim[1]; K <- tjdim[2] @@ -42,7 +42,7 @@ if (L < 500) { truncated <- FALSE svd.method <- "eigen" - } else if (neig > L /2) { + } else if (neig > L / 2) { # Check, whether desired eigentriples amount is too huge if (L < 500) { svd.method <- "eigen" @@ -63,7 +63,7 @@ ssa <- function(x, column.projector = "none", row.projector = "none", column.oblique = "identity", row.oblique = "identity", ..., - kind = c("1d-ssa", "2d-ssa", "nd-ssa", "toeplitz-ssa", "mssa", "cssa"), + kind = c("1d-ssa", "2d-ssa", "nd-ssa", "toeplitz-ssa", "mssa", "cssa", "cmssa"), circular = FALSE, svd.method = c("auto", "nutrlan", "propack", "svd", "eigen", "rspectra", "primme", "irlba", "rsvd"), force.decompose = TRUE) { @@ -82,9 +82,13 @@ ssa <- function(x, ## Provide some sane defaults, e.g. complex inputs should default to cssa if (missing(kind)) { - if (is.complex(x)) + ismssa <- (inherits(x, "mts") || inherits(x, "data.frame") || inherits(x, "list") || inherits(x, "series.list")) + iscssa <- is.complex(x) || ismssa && any(sapply(x, is.complex)) + if (iscssa && ismssa) + kind <- "cmssa" + else if (iscssa) kind <- "cssa" - else if (inherits(x, "mts") || inherits(x, "data.frame") || inherits(x, "list") || inherits(x, "series.list")) + else if (ismssa) kind <- "mssa" else if (is.matrix(x)) kind <- "2d-ssa" @@ -105,6 +109,8 @@ ssa <- function(x, kind <- c("2d-ssa", "nd-ssa") else kind <- "nd-ssa" + } else if (identical(kind, "cmssa")) { + kind <- c("cmssa", "mssa", "cssa") } else if (identical(kind, "mssa")) { ## Nothing special here (yet!) } else if (identical(kind, "cssa")) { From 328e1579a427a22497a6f0ace335d832ee86445a Mon Sep 17 00:00:00 2001 From: Khromov Nikita Date: Sat, 29 Nov 2025 22:21:28 +0300 Subject: [PATCH 2/2] Add documentation and simple tests for cmssa --- NAMESPACE | 2 ++ R/cmhankel.R | 2 -- R/forecast.R | 33 ++++++++++++++++++++++ man/calcv.Rd | 2 ++ man/decompose.Rd | 2 ++ man/forecast.Rd | 8 ++++++ man/hankel.Rd | 1 + man/plot.reconstruct.Rd | 5 +++- man/reconstruct.Rd | 4 +-- man/rforecast.Rd | 17 +++++++---- man/ssa.Rd | 48 +++++++++++++++++-------------- man/vforecast.Rd | 8 +++--- tests/testthat/test-cmssa.R | 56 +++++++++++++++++++++++++++++++++++++ 13 files changed, 152 insertions(+), 36 deletions(-) create mode 100644 tests/testthat/test-cmssa.R diff --git a/NAMESPACE b/NAMESPACE index e909ff92..e2b95298 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -119,10 +119,12 @@ S3method("forecast", "toeplitz.ssa") S3method("predict", "1d.ssa") S3method("predict", "toeplitz.ssa") S3method("predict", "mssa") +S3method("predict", "cmssa") S3method("rforecast", "1d.ssa") S3method("rforecast", "toeplitz.ssa") S3method("rforecast", "mssa") S3method("rforecast", "cssa") +S3method("rforecast", "cmssa") S3method("rforecast", "pssa.1d.ssa") S3method("vforecast", "1d.ssa") S3method("vforecast", "toeplitz.ssa") diff --git a/R/cmhankel.R b/R/cmhankel.R index 3068bff0..1226eb6c 100644 --- a/R/cmhankel.R +++ b/R/cmhankel.R @@ -278,10 +278,8 @@ plot.cmssa.reconstruction <- function(x, ...) { title <- dots$main dots$main <- NULL - oldpar <- par(mfrow = c(1, 2)) do.call(plot.mssa.reconstruction, c(list(x = rx, main = paste(title, "(Real part)")), dots)) do.call(plot.mssa.reconstruction, c(list(x = ix, main = paste(title, "(Imaginary part)")), dots)) - par(oldpar) } .init.fragment.cmssa <- function(this) diff --git a/R/forecast.R b/R/forecast.R index 7b8b17c3..e384f6cd 100644 --- a/R/forecast.R +++ b/R/forecast.R @@ -271,6 +271,25 @@ rforecast.mssa <- function(x, groups, len = 1, invisible(out) } +rforecast.cmssa <- function(x, + groups, + len = 1, + base = c("reconstructed", "original"), + direction = c("row", "column"), + only.new = TRUE, + ..., + drop = TRUE, + drop.attributes = FALSE, + cache = TRUE) { + direction <- match.arg(direction) + if (identical(direction, "row")) + stop("Recurrent Forecast for CMSSA with direction row is not implemented yet.") + + rforecast.mssa(x = x, groups = groups, len = len, base = base, direction = direction, + only.new = only.new, ..., drop = drop, drop.attributes = drop.attributes, + cache = cache) +} + .shift.matrix.1d <- function(U, ...) { wmask <- rep(TRUE, nrow(U)) .shift.matrix(U, wmask, ndim = 1, ...) @@ -433,6 +452,19 @@ vforecast.mssa <- function(x, groups, len = 1, invisible(out) } +vforecast.cmssa <- function(x, + groups, + len = 1, + base = c("reconstructed", "original"), + direction = c("row", "column"), + only.new = TRUE, + ..., + drop = TRUE, + drop.attributes = FALSE, + cache = TRUE) { + stop("Vector Forecast for CMSSA is not implemented yet.") +} + bforecast.1d.ssa <- function(x, groups, len = 1, R = 100, level = 0.95, type = c("recurrent", "vector"), @@ -610,6 +642,7 @@ forecast.1d.ssa <- function(object, "bforecast.toeplitz.ssa" <- `bforecast.1d.ssa`; "forecast.toeplitz.ssa" <- `forecast.1d.ssa`; "predict.toeplitz.ssa" <- `predict.1d.ssa`; +"predict.cmssa" <- `predict.mssa`; "lrr.mssa" <- `lrr.1d.ssa` diff --git a/man/calcv.Rd b/man/calcv.Rd index 2f933c92..80c6c19a 100644 --- a/man/calcv.Rd +++ b/man/calcv.Rd @@ -2,6 +2,7 @@ \alias{calc.v} \alias{calc.v.ssa} \alias{calc.v.cssa} +\alias{calc.v.cmssa} \title{Calculate Factor Vector(s)} @@ -13,6 +14,7 @@ \usage{ \method{calc.v}{ssa}(x, idx, \dots) \method{calc.v}{cssa}(x, idx, \dots) +\method{calc.v}{cmssa}(x, idx, \dots) } \arguments{ diff --git a/man/decompose.Rd b/man/decompose.Rd index 10102fa2..d15d5468 100644 --- a/man/decompose.Rd +++ b/man/decompose.Rd @@ -3,6 +3,7 @@ \alias{decompose.default} \alias{decompose.ssa} \alias{decompose.cssa} +\alias{decompose.cmssa} \alias{decompose.toeplitz.ssa} \title{Perform SSA Decomposition} @@ -14,6 +15,7 @@ \method{decompose}{ssa}(x, neig = NULL, \dots, force.continue = FALSE) \method{decompose}{toeplitz.ssa}(x, neig = NULL, \dots, force.continue = FALSE) \method{decompose}{cssa}(x, neig = NULL, \dots, force.continue = FALSE) +\method{decompose}{cmssa}(x, neig = NULL, \dots, force.continue = FALSE) } \arguments{ diff --git a/man/forecast.Rd b/man/forecast.Rd index 9f042dcd..ad96c6ac 100644 --- a/man/forecast.Rd +++ b/man/forecast.Rd @@ -5,6 +5,7 @@ \alias{predict.ssa} \alias{predict.1d.ssa} \alias{predict.mssa} +\alias{predict.cmssa} \alias{predict.toeplitz.ssa} \title{Perform SSA forecasting of series} @@ -47,8 +48,15 @@ direction = c("column", "row"), \dots, drop = TRUE, drop.attributes = FALSE, cache = TRUE) +\method{predict}{cmssa}(object, + groups, len = 1, + method = c("recurrent", "vector"), + direction = c("column", "row"), + \dots, + drop = TRUE, drop.attributes = FALSE, cache = TRUE) } + \arguments{ \item{object}{SSA object holding the decomposition} \item{groups}{list, the grouping of eigentriples to be used in the forecast} diff --git a/man/hankel.Rd b/man/hankel.Rd index e9f8d53b..377002b7 100644 --- a/man/hankel.Rd +++ b/man/hankel.Rd @@ -4,6 +4,7 @@ \alias{hcols} \alias{hrows} \alias{hankel} +\alias{mhankel} \alias{hmatmul} \title{Hankel matrices operations.} diff --git a/man/plot.reconstruct.Rd b/man/plot.reconstruct.Rd index 8d33a553..f9d37610 100644 --- a/man/plot.reconstruct.Rd +++ b/man/plot.reconstruct.Rd @@ -4,6 +4,7 @@ \alias{plot.1d.ssa.reconstruction} \alias{plot.toeplitz.ssa.reconstruction} \alias{plot.mssa.reconstruction} +\alias{plot.cmssa.reconstruction} \alias{plot.2d.ssa.reconstruction} \alias{plot.nd.ssa.reconstruction} @@ -33,6 +34,7 @@ base.series = NULL, add.original = TRUE, add.residuals = TRUE) +\method{plot}{cmssa.reconstruction}(x, \dots) \method{plot}{2d.ssa.reconstruction}(x, \dots, type = c("raw", "cumsum"), base.series = NULL, @@ -55,7 +57,8 @@ multidimentional array to draw.} \item{type}{Type of the plot (see 'Details' for more information)} \item{\dots}{Arguments to be passed to methods, such as graphical - parameters} + parameters. In `cmssa` case all arguments will be propagated to the + reconstruction plot method of the `mssa` class} \item{plot.method}{Plotting method to use: either ordinary all-in-one via matplot or xyplot, or native plotting method of the input time series} diff --git a/man/reconstruct.Rd b/man/reconstruct.Rd index 2823577c..8805b556 100644 --- a/man/reconstruct.Rd +++ b/man/reconstruct.Rd @@ -32,7 +32,7 @@ Algorithm 5.2 for 2D-SSA and Algorithm 5.6 for Shaped 2D-SSA. Fast implementation of reconstruction with the help of FFT is described in Korobeynikov (2010) - for the 1D case and in Section 6.2 (Rank-one quasi-hankelization) of Golyandina et al (2015) + for the 1D case and in Section 6.2 (Rank-one quasi-hankelization) of Golyandina et al (2015) for the general case. } @@ -67,7 +67,7 @@ \seealso{ \code{\link{Rssa}} for an overview of the package, as well as, - \code{\link[Rssa:ssa-input]{ssa-input}}, + \code{\link[Rssa:ssa-input]{ssa-input}}, \code{\link[Rssa:ssa]{ssa}}, \code{\link[Rssa:plot.reconstruction]{plot.reconstruction}}, } diff --git a/man/rforecast.Rd b/man/rforecast.Rd index 62ec4221..f0f865a0 100644 --- a/man/rforecast.Rd +++ b/man/rforecast.Rd @@ -6,6 +6,7 @@ \alias{rforecast.toeplitz.ssa} \alias{rforecast.mssa} \alias{rforecast.cssa} +\alias{rforecast.cmssa} \alias{rforecast.pssa.1d.ssa} \title{Perform recurrent SSA forecasting of the series} @@ -26,6 +27,9 @@ \method{rforecast}{cssa}(x, groups, len = 1, base = c("reconstructed", "original"), only.new = TRUE, reverse = FALSE, \dots, drop = TRUE, drop.attributes = FALSE, cache = TRUE) +\method{rforecast}{cmssa}(x, groups, len = 1, base = c("reconstructed", "original"), + direction = c("row", "column"), only.new = TRUE, \dots, drop = TRUE, + drop.attributes = FALSE, cache = TRUE) \method{rforecast}{pssa.1d.ssa}(x, groups, len = 1, base = c("reconstructed", "original"), only.new = TRUE, reverse = FALSE, \dots, drop = TRUE, drop.attributes = FALSE, cache = TRUE) @@ -38,7 +42,8 @@ \item{base}{series used as a 'seed' of forecast: original or reconstructed according to the value of \code{groups} argument} \item{direction}{direction of forecast in multichannel SSA case, "column" - stands for so-called L-forecast and "row" stands for K-forecast} + stands for so-called L-forecast and "row" stands for K-forecast + (for complex case only "column" variant is impemented)} \item{only.new}{logical, if 'TRUE' then only forecasted values are returned, whole series otherwise} \item{reverse}{logical, direction of forecast in 1D SSA case, 'FALSE' @@ -61,12 +66,12 @@ sequentialy projects the incomplete embedding vectors (from either the original or the reconstructed series) onto the subspace spanned by the selected eigentriples of the decomposition to derive the missed - (last) values of the such vectors. Then the filled value + (last) values of the such vectors. Then the filled value In such a way the forecasting elements are produced on one-by-one basis. It is shown in Golyandina et al (2001) that this approach corresponds - to application of a linear recurrence formula (the same formula as + to application of a linear recurrence formula (the same formula as described in \code{\link{lrr}}) to initial data taken from either the original or the reconstructed series. @@ -82,9 +87,9 @@ and is based on a multivariate LRR. Forecast uses the formulae from Golyandina and Stepanov (2005) and Golyandina et.al (2015). - For details of 1D-SSA recurrent forecasting, see Section 3.2.1.2 and + For details of 1D-SSA recurrent forecasting, see Section 3.2.1.2 and Algorithm 3.5 in Golyandina et al (2018). - For details of MSSA recurrent forecasting, see Section 4.3.1.2 and + For details of MSSA recurrent forecasting, see Section 4.3.1.2 and Algorithm 4.4 (column forecasting). } @@ -101,7 +106,7 @@ Golyandina, N., Nekrutkin, V. and Zhigljavsky, A. (2001): \emph{Analysis of Time Series Structure: SSA and related techniques.} Chapman and Hall/CRC. ISBN 1584881941 - + Golyandina, N., Korobeynikov, A., Shlemov, A. and Usevich, K. (2015): \emph{Multivariate and 2D Extensions of Singular Spectrum Analysis with the Rssa Package}. Journal of Statistical Software, Vol. 67, Issue 2. diff --git a/man/ssa.Rd b/man/ssa.Rd index abd6a4ab..62286f34 100644 --- a/man/ssa.Rd +++ b/man/ssa.Rd @@ -18,9 +18,9 @@ ssa(x, column.oblique = "identity", row.oblique = "identity", \dots, - kind = c("1d-ssa", "2d-ssa", "nd-ssa", "toeplitz-ssa", "mssa", "cssa"), + kind = c("1d-ssa", "2d-ssa", "nd-ssa", "toeplitz-ssa", "mssa", "cssa", "cmssa"), circular = FALSE, - svd.method = c("auto", "nutrlan", "propack", "svd", "eigen", + svd.method = c("auto", "nutrlan", "propack", "svd", "eigen", "rspectra", "primme", "irlba", "rsvd"), force.decompose = TRUE) } @@ -40,7 +40,7 @@ ssa(x, \item{\dots}{further arguments passed to \code{\link[Rssa:decompose.ssa]{decompose}} } \item{kind}{SSA method. This includes ordinary 1d SSA, 2d SSA, Toeplitz variant of 1d SSA, multichannel variant of SSA and complex - SSA} + SSA (1d and multichannel variants)} \item{circular}{logical vector of one or two elements, describes series topology for 1d SSA and Toeplitz SSA or field topology for 2d SSA. 'TRUE' means series circularity for 1d case or circularity by a corresponding @@ -50,7 +50,7 @@ ssa(x, \item{column.projector, row.projector}{column and row signal subspaces projectors for SSA with projection. See `Details' for information about methods of projectors specification} - \item{column.oblique, row.oblique}{column and row matrix weights for + \item{column.oblique, row.oblique}{column and row matrix weights for Weighted Oblique SSA. See `Details' for information about how to use this feature} \item{force.decompose}{logical, if 'TRUE' then the decomposition is performed before return.} @@ -67,7 +67,7 @@ ssa(x, constructs the SSA object filling all necessary internal structures and performing the decomposition if necessary. For the comprehensive description of SSA modifications and their algorithms - see Golyandina et al (2018). + see Golyandina et al (2018). \subsection{Variants of SSA}{ The following implementations of the SSA method are supported @@ -84,6 +84,7 @@ ssa(x, several time series (possible of unequal length). See Golyandina and Stepanov (2005).} \item{cssa}{Complex variant of 1d SSA.} + \item{cmssa}{Complex variant of multichannel SSA.} \item{2d-ssa}{2d SSA for decomposition of images and arrays. See Golyandina and Usevich (2009) and Golyandina et.al (2015) for more information.} \item{nd-ssa}{Multidimensional SSA decomposition for arrays (tensors).} @@ -103,12 +104,12 @@ ssa(x, } These functions are not exported, they defined only for \code{wmask} expression. If one has objects with the same names and wants to use them rather than these functions, - one should use special wrapper function \code{I()} (see 'Examples'). + one should use special wrapper function \code{I()} (see 'Examples'). } \subsection{Projectors specification for SSA with projection}{ Projectors are specified by means of \code{column.projector} and - \code{row.projector} arguments (see Golyandina and Shlemov (2017)). + \code{row.projector} arguments (see Golyandina and Shlemov (2017)). Each may be a matrix of orthonormal (otherwise QR orthonormalization process will be perfomed) basis of projection subspace, or single integer, which will be interpreted as @@ -133,15 +134,15 @@ ssa(x, } \subsection{Weighted Oblique SSA}{ - Corresponding matrix norm weights may be specified for ordinary 1D SSA case - by means of \code{column.oblique} and \code{row.oblique} arguments. These - arguments should be either 'identical' or positive numeric vectors of length - \code{L} and \code{N - L + 1} for \code{column.oblique} and + Corresponding matrix norm weights may be specified for ordinary 1D SSA case + by means of \code{column.oblique} and \code{row.oblique} arguments. These + arguments should be either 'identical' or positive numeric vectors of length + \code{L} and \code{N - L + 1} for \code{column.oblique} and \code{row.oblique} respectively. - Weighted Oblique SSA inside \code{\link[Rssa:cadzow]{Cadzow}} iterations - may improve finite-rank estimation of signal - (see e.g. Cadzow(alpha) iterations in Zvonarev and Golyandina (2017) + Weighted Oblique SSA inside \code{\link[Rssa:cadzow]{Cadzow}} iterations + may improve finite-rank estimation of signal + (see e.g. Cadzow(alpha) iterations in Zvonarev and Golyandina (2017) for more information). } @@ -216,7 +217,7 @@ ssa(x, Proceedings of the 5th St.Petersburg Workshop on Simulation, June 26-July 2, 2005, St. Petersburg State University, St. Petersburg, 293--298. \url{https://www.gistatgroup.com/gus/mssa2.pdf} - + Golyandina, N. and Usevich, K. (2009): \emph{2D-extensions of singular spectrum analysis: algorithm and elements of theory.} In Matrix Methods: Theory, Algorithms, Applications. World Scientific @@ -234,21 +235,21 @@ ssa(x, Analysis for time series}. Springer Briefs in Statistics. Springer. Shlemov, A. and Golyandina, N. (2014): \emph{Shaped extensions of singular - spectrum analysis}. 21st International Symposium on Mathematical - Theory of Networks and Systems, July 7-11, 2014. Groningen, + spectrum analysis}. 21st International Symposium on Mathematical + Theory of Networks and Systems, July 7-11, 2014. Groningen, The Netherlands. p.1813-1820. \url{https://arxiv.org/abs/1507.05286} - + Golyandina, N., Korobeynikov, A., Shlemov, A. and Usevich, K. (2015): \emph{Multivariate and 2D Extensions of Singular Spectrum Analysis with the Rssa Package}. Journal of Statistical Software, Vol. 67, Issue 2. \doi{10.18637/jss.v067.i02} - Golyandina, N. and Shlemov, A. (2017): \emph{Semi-nonparametric singular - spectrum analysis with projection}. Statistics and its Interface, Vol. 10, + Golyandina, N. and Shlemov, A. (2017): \emph{Semi-nonparametric singular + spectrum analysis with projection}. Statistics and its Interface, Vol. 10, Issue 1, p. 47-57. \url{https://arxiv.org/abs/1401.4980} - Zvonarev, N. and Golyandina, N. (2017): \emph{Iterative algorithms for weighted and + Zvonarev, N. and Golyandina, N. (2017): \emph{Iterative algorithms for weighted and unweighted finite-rank time-series approximations}. Statistics and its Interface, Vol. 10, Issue 1, p. 5-18. \url{https://arxiv.org/abs/1507.02751} @@ -322,6 +323,11 @@ s <- ssa(EuStockMarkets[, 1] + 1.0i*EuStockMarkets[, 2], kind = "cssa") r <- reconstruct(s, groups = list(Trend = 1:2)) plot(r) +# CMSSA-based trend extraction +s <- ssa(EuStockMarkets[, 1:2] + 1.0i*EuStockMarkets[, 3:4], kind = "cmssa") +r <- reconstruct(s, groups = list(Trend = 1:2)) +plot(r) + # `co2' decomposition with double projection to linear functions s <- ssa(co2, column.projector = "centering", row.projector = "centering") plot(reconstruct(s, groups = list(trend = seq_len(nspecial(s))))) diff --git a/man/vforecast.Rd b/man/vforecast.Rd index 66f150d6..06043352 100644 --- a/man/vforecast.Rd +++ b/man/vforecast.Rd @@ -55,13 +55,13 @@ formula as described in \code{\link{lrr}} is used for constructing of the last components of the new vectors) and then derives the series out of this extended set of vectors. - - For multichannel SSA, forecast can be constructed in two versions, + + For multichannel SSA, forecast can be constructed in two versions, row and column ones; it uses the formulae from Golyandina et al (2015). - For details of 1D-SSA recurrent forecasting, see Section 3.2.1.3 and + For details of 1D-SSA recurrent forecasting, see Section 3.2.1.3 and Algorithm 3.6 in Golyandina et al (2018). - For details of MSSA recurrent forecasting, see Section 4.3.1.3 and + For details of MSSA recurrent forecasting, see Section 4.3.1.3 and Algorithm 4.5 (column forecasting). } diff --git a/tests/testthat/test-cmssa.R b/tests/testthat/test-cmssa.R new file mode 100644 index 00000000..944722a8 --- /dev/null +++ b/tests/testthat/test-cmssa.R @@ -0,0 +1,56 @@ +library(testthat) +library(Rssa) +context("CMSSA") + +test.cases <- c("linear", "periodic") + +for (case in test.cases) { +test_that(sprintf("Complex MSSA reconstruction is correct, %s case", case), { + N <- 100; P <- 5; coefs <- seq(from = 1, to = P) + v <- switch(case, + linear = (1 + 1i) * seq(N) %o% coefs, + periodic = exp(1i * seq(N) / 12) %o% coefs) + s <- ssa(v, kind = "cmssa") + eps <- sqrt(.Machine$double.eps) + + rank <- switch(case, + linear = 2, + periodic = 1) + + rec <- reconstruct(s, groups = list(1:rank))$F1 + expect_true(sqrt(mean(abs(rec - v)^2)) < eps) +}) + +test_that(sprintf("Complex MSSA recurrent column forecast is correct, %s case", case), { + N <- 110; P <- 5; coefs <- seq(from = 1, to = P); len <- 10 + v <- switch(case, + linear = (1 + 1i) * seq(N) %o% coefs, + periodic = exp(1i * seq(N) / 12) %o% coefs) + s <- ssa(v[1:(N - len)], kind = "cmssa") + eps <- sqrt(.Machine$double.eps) + + rank <- switch(case, + linear = 2, + periodic = 1) + + pred <- rforecast(s, groups = list(1:rank), direction = "column", len = len) + expect_true(sqrt(mean(abs(pred - v[(N - len + 1):N])^2)) < eps) +}) +} + +test_that("Built-in SVD and Primme SVD yield same reconstructions in Complex MSSA", { + N <- 100; P <- 5; + v <- matrix(rnorm(N * P) + 1i * rnorm(N * P), ncol = P) + s_default <- ssa(v, kind = "cmssa") + s_primme <- ssa(v, kind = "cmssa", svd = "primme") + eps <- sqrt(.Machine$double.eps) + + ranks <- 1:5 + + for (rank in ranks) { + rec_default <- reconstruct(s_default, groups = list(1:rank))$F1 + rec_primme <- reconstruct(s_primme, groups = list(1:rank))$F1 + expect_true(sqrt(mean(abs(rec_default - rec_primme)^2)) < eps) + } + +}) \ No newline at end of file