Permalink
Browse files

Add recurrent forecast routine.

Thanks goes to Alex Shlemov for the implementation
  • Loading branch information...
1 parent ae864c7 commit af9c90e13320196e99f6ad6441c252bccc5897fa @asl committed Feb 19, 2012
Showing with 156 additions and 2 deletions.
  1. +4 −1 NAMESPACE
  2. +2 −0 R/common.R
  3. +85 −0 R/forecast.R
  4. +1 −1 man/reconstruct.Rd
  5. +64 −0 man/rforecast.Rd
View
@@ -35,7 +35,9 @@ export(clone,
tmatmul,
tcols,
trows,
- is.tmat)
+ is.tmat,
+# Forecast stuff
+ rforecast)
S3method("clone", ssa)
S3method("decompose", ssa)
@@ -72,3 +74,4 @@ S3method("wcor", "toeplitz-ssa")
S3method("plot", wcor.matrix)
S3method("clusterify", wcor.matrix)
S3method("clusterify", ssa)
+S3method("rforecast", "1d-ssa")
View
@@ -112,6 +112,8 @@ clusterify <- function(this, ...)
UseMethod("clusterify");
calc.v <- function(this, ...)
UseMethod("calc.v");
+rforecast <- function(this, ...)
+ UseMethod("rforecast")
.object.size <- function(this, ...)
UseMethod(".object.size")
View
@@ -0,0 +1,85 @@
+# R package for Singular Spectrum Analysis
+# Copyright (c) 2012 Alexandr Shlemov <shlemovalex@gmail.com>
+# Copyright (c) 2012 Anton Korobeynikov <asl@math.spbu.ru>
+#
+# 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.
+
+basis2lrf <- function(U) {
+ N <- nrow(U);
+ lpf <- U %*% t(U[N, , drop = FALSE]);
+
+ (lpf[-N]) / (1 - lpf[N])
+}
+
+apply.lrf <- function(F, lrf, len = 1) {
+ N <- length(F)
+ r <- length(lrf)
+
+ # Sanity check of inputs
+ if (r > N)
+ stop("Wrong length of LRF")
+
+ # Run the actual LRF
+ F <- c(F, rep(NA, len))
+ for (i in 1:len)
+ F[N+i] <- sum(F[(N+i-r) : (N+i-1)]*lrf)
+
+ F
+}
+
+"rforecast.1d-ssa" <- function(this, groups, len = 1,
+ base = c("reconstructed", "original"),
+ ..., cache = TRUE) {
+ base <- match.arg(base)
+ if (missing(groups))
+ groups <- as.list(1:min(nlambda(this), nu(this)))
+
+ # Determine the upper bound of desired eigentriples
+ desired <- max(unlist(groups));
+
+ # Continue decomposition, if necessary
+ if (desired > min(nlambda(this), nu(this)))
+ decompose(this, ..., neig = desired);
+
+ # Grab the reconstructed series if we're basing on them
+ if (identical(base, "reconstructed"))
+ r <- reconstruct(this, groups = groups, ..., cache = cache)
+
+ out <- list()
+ for (i in seq_along(groups)) {
+ group <- groups[[i]]
+
+ # Calculate the LRF corresponding to group
+ U <- .get(this, "U")[, group, drop = FALSE]
+ lrf <- basis2lrf(U)
+
+ # Calculate the forecasted values
+ out[[i]] <- apply.lrf(if (identical(base, "reconstructed")) r[[i]] else .get(this, "F"),
+ lrf, len)
+ # FIXME: try to fixup the attributes
+ }
+
+ names(out) <- paste(sep = "", "F", 1:length(groups))
+
+ out
+}
+
+rforecast.ssa <- function(x, groups, len = 1,
+ base = c("reconstructed", "original"),
+ ..., cache = TRUE) {
+ stop("generic recurrent forecast not implemented yet!")
+}
View
@@ -24,7 +24,7 @@
\value{
List of vectors of reconstructed series. Elements of the list are
- named 'F1'', 'F2', and so on.
+ named 'F1', 'F2', and so on.
}
\examples{
View
@@ -0,0 +1,64 @@
+\name{rforecast}
+\alias{rforecast}
+\alias{rforecast.default}
+\alias{rforecast.ssa}
+\alias{rforecast.1d-ssa}
+\title{Perform recurrent SSA forecasting of the series}
+
+\description{
+ Perform recurrent SSA forecasting of the series.
+}
+
+\usage{
+\method{rforecast}{ssa}(x, groups, len = 1, base = c("reconstructed", "original"), \dots, cache = TRUE)
+}
+
+\arguments{
+ \item{x}{SSA object holding the decomposition}
+ \item{groups}{the grouping of eigentriples to be used in the forecast}
+ \item{len}{the desired length of the forecasted series}
+ \item{base}{series used as a 'seed' of forecast: original or
+ reconstructed according to the value of \code{groups} argument}
+ \item{\dots}{additional arguments passed to \code{\link{reconstruct}} routines}
+ \item{cache}{logical, if 'TRUE' then intermediate results will be
+ cached in the SSA object.}
+}
+
+\details{
+ The routines applies the recurrent SSA forecasting algorithm to
+ produce the new series which are expected to 'continue' the current
+ series on the basis of the decomposition given. The algorithm
+ sequentialy projects the incomplete embedding vectors (either original
+ or from reconstructed series) onto the subspace spanned by the
+ selected eigentriples of the decomposition to derive the missed
+ (ending) values of the such vectors.
+
+ In such a way the forecasted elements of the series are produced on
+ one-by-one basis.
+}
+
+\note{
+ In fact the current implementation of the algorithm instead of using
+ the direct projections of the series calculates the so-called Linear
+ Recurrent Formula which governs the series and uses it for the
+ forecast.
+}
+
+\value{
+ List of vectors of forecasted series. Elements of the list are
+ named 'F1', 'F2', and so on.
+}
+
+\seealso{
+ \code{\link[svd:svd]{svd}},
+ \code{\link[Rssa:new.ssa]{new.ssa}},
+ \code{\link[Rssa:reconstruct]{reconstruct}}
+}
+
+\examples{
+# Decompose 'co2' series with default parameters
+s <- new.ssa(co2)
+# Produce 5 forecasted values of the series using the first 3
+# eigentriples as a base space for the forecast.
+rforecast(s, groups = list(1:3), len = 5)
+}

0 comments on commit af9c90e

Please sign in to comment.