Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Forecast #5

Closed
wants to merge 8 commits into from

2 participants

@eodus

Two commits with some simple corections.

Really, I only add "asp=1" in plot.lrf functon and change method of evaluating roots to "eigenvalues of companion-matrix" (polyroot is unstable for L >= 100). And add semicolons and some comments.

UPD:

  • 5 commits
@asl
Owner

Should we close this?

@asl asl closed this
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
Showing with 103 additions and 75 deletions.
  1. +103 −75 R/forecast.R
View
178 R/forecast.R
@@ -1,5 +1,5 @@
# R package for Singular Spectrum Analysis
-# Copyright (c) 2012 Alexandr Shlemov <shlemovalex@gmail.com>
+# Copyright (c) 2012 Alexander Shlemov <shlemovalex@gmail.com>
# Copyright (c) 2012 Anton Korobeynikov <asl@math.spbu.ru>
#
# This program is free software; you can redistribute it
@@ -22,61 +22,83 @@ basis2lrf <- function(U) {
N <- nrow(U);
lpf <- U %*% t(U[N, , drop = FALSE]);
- (lpf[-N]) / (1 - lpf[N])
+ (lpf[-N]) / (1 - lpf[N]);
}
"lrf.1d-ssa" <- function(this, group, ...) {
- U <- .get(this, "U")[, group, drop = FALSE]
+ # Determine the upper bound of desired eigentriples
+ desired <- max(group);
+
+ # Continue decomposition, if necessary
+ if (desired > nu(this))
+ decompose(this, ..., neig = desired);
+
+ U <- .get(this, "U")[, group, drop = FALSE];
- res <- basis2lrf(U)
- class(res) <- "lrf"
+ res <- basis2lrf(U);
+ class(res) <- "lrf";
- res
+ res;
+}
+
+
+companion.matrix.lrf <- function(x) {
+ n <- length(x) + 1;
+ res <- matrix(0, n, n);
+ res[-n, n] <- x;
+ res[seq(from = 2, by = n + 1, length.out = n - 1)] <- 1;
+
+ res;
}
roots.lrf <- function(x) {
- polyroot(c(-x, 1))
+ # polyroot(c(-x, 1));
+ # Much more complicated but much more stable
+ eigen(companion.matrix.lrf(x), only.values = TRUE)$values;
}
plot.lrf <- function(x, ..., raw = FALSE) {
- r <- roots(x)
+ r <- roots(x);
if (raw) {
- plot(r, ...)
+ plot(r, ...);
} else {
- xlim <- range(c(Re(r), +1, -1))
- ylim <- range(c(Im(r), +1, -1))
+ xlim <- range(c(Re(r), +1, -1));
+ ylim <- range(c(Im(r), +1, -1));
plot(r, ...,
xlim = xlim, ylim = ylim,
main = "Roots of Linear Recurrence Formula",
xlab = "Real Part",
- ylab = "Imaginary Part")
- symbols(0, 0, circles = 1, add = TRUE, inches = FALSE)
+ ylab = "Imaginary Part",
+ asp = 1);
+ symbols(0, 0, circles = 1, add = TRUE, inches = FALSE);
}
}
apply.lrf <- function(F, lrf, len = 1) {
- N <- length(F)
- r <- length(lrf)
+ N <- length(F);
+ r <- length(lrf);
# Sanity check of inputs
if (r > N)
- stop("Wrong length of LRF")
+ stop("Wrong length of LRF");
# Run the actual LRF
- F <- c(F, rep(NA, len))
+ F <- c(F, rep(NA, len));
+
+ # TODO Rewrite this function on C
for (i in 1:len)
- F[N+i] <- sum(F[(N+i-r) : (N+i-1)]*lrf)
+ F[N+i] <- sum(F[(N+i-r) : (N+i-1)]*lrf);
- F
+ F;
}
"rforecast.1d-ssa" <- function(this, groups, len = 1,
base = c("reconstructed", "original"),
..., cache = TRUE) {
- base <- match.arg(base)
+ base <- match.arg(base);
if (missing(groups))
- groups <- as.list(1:min(nlambda(this), nu(this)))
+ groups <- as.list(1:min(nlambda(this), nu(this)));
# Determine the upper bound of desired eigentriples
desired <- max(unlist(groups));
@@ -87,52 +109,40 @@ apply.lrf <- function(F, lrf, len = 1) {
# Grab the reconstructed series if we're basing on them
if (identical(base, "reconstructed"))
- r <- reconstruct(this, groups = groups, ..., cache = cache)
+ r <- reconstruct(this, groups = groups, ..., cache = cache);
- out <- list()
+ out <- list();
for (i in seq_along(groups)) {
- group <- groups[[i]]
+ group <- groups[[i]];
# Calculate the LRF corresponding to group
- lrf <- lrf(this, group)
+ lf <- lrf(this, group);
# Calculate the forecasted values
out[[i]] <- apply.lrf(if (identical(base, "reconstructed")) r[[i]] else .get(this, "F"),
- lrf, len)
+ lf, len);
# FIXME: try to fixup the attributes
+ # FIXME: ASH It's impossible --- different lengths =(
}
- names(out) <- paste(sep = "", "F", 1:length(groups))
+ names(out) <- paste(sep = "", "F", 1:length(groups));
- out
+ # Forecasted series can be pretty huge...
+ invisible(out);
}
+
"vforecast.1d-ssa" <- function(this, groups, len = 1,
- ...) {
- L <- this$window
+ ...,
+ cache = TRUE) {
+ L <- this$window;
K <- this$length - L + 1;
- L.s <- min(L, K);
N <- K + L - 1 + len + L - 1;
N.res <- K + L - 1 + len;
- dv <- c(1:(L.s-1), rep(L.s, N-2*L.s+2), (L.s-1):1);
-
- convolve.open <- function(F, G) {
- NF <- length(F)
- NG <- length(G)
-
- NN <- nextn(NF+NG-1)
- ZFZ <- c(rep(0, NG-1), F, rep(0, NN-(NF+NG-1)));
- GZ <- c(G, rep(0, NN-NG));
-
- res <- fft(fft(ZFZ)*Conj(fft(GZ)), inverse = TRUE)/NN;
-
- Re(res)[1:(NF+NG-1)];
- }
-
if (missing(groups))
- groups <- as.list(1:min(nlambda(this), nu(this)))
+ groups <- as.list(1:min(nlambda(this), nu(this)));
# Determine the upper bound of desired eigentriples
desired <- max(unlist(groups));
@@ -144,40 +154,58 @@ apply.lrf <- function(F, lrf, len = 1) {
lambda <- .get(this, "lambda");
U <- .get(this, "U");
- V <- if (nv(this) >= desired) .get(this, "V") else NULL
-
- out <- list()
+ if (nv(this) >= desired) {
+ V <- .get(this, "V");
+ } else {
+ if (cache){
+ V <- cbind(.get(this, "V", allow.null = TRUE),
+ calc.v(this, idx = (nv(this) + 1):desired));
+ .set(this, "V", V);
+ } else {
+ V <- NULL;
+ }
+ }
+
+ # Make hankel matrix for fast hankelization (we use it for plan)
+ h <- new.hmat(double(N), L);
+
+ out <- list();
for (i in seq_along(groups)) {
- ET <- unique(groups[[i]])
-
- Vl <- matrix(NA, N.res, length(ET));
- Uet <- U[ , ET, drop=FALSE];
-
- Vl[1:K, ] <- (if (is.null(V))
- calc.v(this, idx = ET)
- else V[ , ET]) %*% diag(lambda[ET], nrow = length(ET))
-
- U.head <- Uet[-L, , drop=FALSE];
- U.tail <- Uet[-1, , drop=FALSE];
- P <- solve(t(U.head) %*% U.head, t(U.head) %*% U.tail);
-
- for (j in (K+1):(K+len+L-1)) {
- Vl[j, ] <- P %*% Vl[j-1, ];
+ group <- unique(groups[[i]]);
+
+ Uet <- U[, group, drop = FALSE];
+ Vet <- if (is.null(V)) calc.v(this, idx = group) else V[, group, drop = FALSE];
+
+ Z <- rbind(t(lambda[group] * t(Vet)), matrix(NA, len + L - 1, length(group)));
+
+ U.head <- Uet[-L, , drop = FALSE];
+ U.tail <- Uet[-1, , drop = FALSE];
+ Pi <- Uet[L, ];
+ tUhUt <- t(U.head) %*% U.tail;
+ P <- tUhUt + 1 / (1 - sum(Pi^2)) * Pi %*% (t(Pi) %*% tUhUt);
+
+ for(j in (K + 1):(K + len + L - 1)) {
+ Z[j, ] <- P %*% Z[j - 1, ];
}
-
- res <- rep(0, N);
- for (j in 1:length(ET)) {
- res <- res + convolve.open(Vl[ , j], rev(Uet[ , j]));
+
+ res <- double(N);
+ for(j in seq_along(group)) {
+ res <- res + .hankelize.one.hankel(Uet[ , j], Z[ , j], h);
}
-
- out[[i]] <- (res/dv)[1:N.res];
- };
-
+
+ out[[i]] <- res[1:N.res];
+ }
+
names(out) <- paste(sep = "", "F", 1:length(groups));
-
- out
+
+ # Forecasted series can be pretty huge...
+ invisible(out);
}
+"lrf.toeplitz-ssa" <- get("lrf.1d-ssa");
+"vforecast.toeplitz-ssa" <- get("vforecast.1d-ssa");
+"rforecast.toeplitz-ssa" <- get("rforecast.1d-ssa");
+
rforecast.ssa <- function(x, groups, len = 1,
base = c("reconstructed", "original"),
..., cache = TRUE) {
Something went wrong with that request. Please try again.