library("Rssa")
data("AustralianWine", package = "Rssa")
wine <- window(AustralianWine, end = time(AustralianWine)[174])

### Begin fragment: fort_rec
fort <- wine[, "Fortified"]
s.fort <- ssa(fort, L = 84, kind = "1d-ssa")
r.fort <- reconstruct(s.fort, groups = list(Trend = 1,
                                            Seasonality = 2:11))
plot(r.fort, add.residuals = TRUE, add.original = TRUE,
     plot.method = "xyplot",
     superpose = TRUE, auto.key = list(columns = 2))
### End fragment: fort_rec

### Begin fragment: fort_identific
plot(s.fort, type = "vectors", idx = 1:8)
plot(s.fort, type = "paired", idx = 2:11, plot.contrib = FALSE)
parestimate(s.fort, groups = list(2:3, 4:5), method = "esprit-ls")
plot(wcor(s.fort, groups = 1:30),
          scales = list(at = c(10, 20, 30)))
plot(reconstruct(s.fort, add.residuals = FALSE, add.original = FALSE,
                 groups = list(G12 = 2:3, G4 = 4:5, G6 = 6:7, G2.4 = 8:9)))
### End fragment: fort_identific

### Begin fragment: fort_forecast
f.fort <- vforecast(s.fort,
                    groups = list(Trend = 1, Signal = 1:11),
                    len = 60, only.new = TRUE)
plot(cbind(fort, f.fort$Signal, f.fort$Trend),
     plot.type = "single", col = c("black", "red", "blue"), ylab = NULL)
### End fragment: fort_forecast

### Begin fragment: wineFortDry_rec
wineFortDry <- wine[, c("Fortified", "Drywhite")]
L <- 84
s.wineFortDry <- ssa(wineFortDry, L = L, kind = "mssa")
r.wineFortDry <- reconstruct(s.wineFortDry,
                             groups = list(Trend = c(1, 6),
                                           Seasonality = c(2:5, 7:12)))
plot(r.wineFortDry, add.residuals = FALSE,
     plot.method = "xyplot",
     superpose = TRUE, auto.key = list(columns = 3))
### End fragment: wineFortDry_rec

### Begin fragment: wineFortDry_summary
summary(s.wineFortDry)
### End fragment: wineFortDry_summary

### Begin fragment: wineFortDry_identific
plot(s.wineFortDry, type = "vectors", idx = 1:8)
plot(s.wineFortDry, type = "paired", idx = 2:11, plot.contrib = FALSE)
parestimate(s.wineFortDry, groups = list(2:3, 4:5), method = "esprit-ls")
plot(wcor(s.wineFortDry, groups = 1:30), scales = list(at = c(10, 20, 30)))
### End fragment: wineFortDry_identific

### Begin fragment: wineFortDry_forecast
f.wineFortDry <- rforecast(s.wineFortDry, groups = list(1, 1:12),
                           len = 60, only.new = TRUE)
par(mfrow = c(2, 1))
plot(cbind(wineFortDry[, "Fortified"], f.wineFortDry$F2[, "Fortified"]),
     plot.type = "single",
     col = c("black", "red"), ylab = "Fort")
plot(cbind(wineFortDry[, "Drywhite"], f.wineFortDry$F2[, "Drywhite"]),
     plot.type = "single",
     col = c("black", "red"), ylab = "Dry")
par(mfrow = c(1, 1))
### End fragment: wineFortDry_forecast

### Begin fragment: wineFortDry_factor
library("lattice")
L <- 24
s.wineFortDrya <- ssa(wineFortDry, L = L, kind = "mssa")
r.wineFortDrya <- reconstruct(s.wineFortDrya, groups = list(Trend = 1))
tp1 <- plot(r.wineFortDrya,
            add.residuals = FALSE, add.original = TRUE,
            plot.method = "xyplot", aspect = 0.3, superpose = TRUE,
            scales = list(y = list(draw = FALSE)),
            auto.key = "", xlab = "",
            col = c("blue", "violet", "blue", "violet"))
tp2 <- plot(s.wineFortDrya, type = "vectors", vectors = "factor", idx = 1,
            aspect = 0.5, superpose = TRUE,
            scales = list(x = list(draw = TRUE), y = list(draw = FALSE)),
            auto.key = list(columns = 2))
tp3 <- plot(s.wineFortDrya, type = "vectors", vectors = "eigen", idx = 1,
            aspect = 1,
            scales = list(x = list(draw = TRUE), y = list (draw = FALSE)))
plot(tp3, split = c(1, 1, 1, 3), more = TRUE)
plot(tp2, split = c(1, 2, 1, 3), more = TRUE)
plot(tp1, split = c(1, 3, 1, 3), more = FALSE)
### End fragment: wineFortDry_factor

### Begin fragment: wineFortRose_scales
wineFortRose <- wine[, c("Fortified", "Rose")]
summary(wineFortRose)
norm.wineFortRosen <- sqrt(colMeans(wineFortRose^2))
wineFortRosen <- sweep(wineFortRose, 2, norm.wineFortRosen, "/")
L <- 84
s.wineFortRosen <- ssa(wineFortRosen, L = L, kind = "mssa")
r.wineFortRosen <- reconstruct(s.wineFortRosen,
                               groups = list(Trend = c(1, 12, 14),
                                             Seasonality = c(2:11, 13)))
s.wineFortRose <- ssa(wineFortRose, L = L, kind = "mssa")
r.wineFortRose <- reconstruct(s.wineFortRose,
                              groups = list(Trend = 1,
                                            Seasonality = 2:11))
wrap.plot <- function(rec, component = 1, series, xlab = "", ylab, ...)
  plot(rec, add.residuals = FALSE, add.original = TRUE,
       plot.method = "xyplot", superpose = TRUE,
       scales = list(y = list(tick.number = 3)),
       slice = list(component = component, series = series),
       xlab = xlab, ylab = ylab, auto.key = "", ...)
trel1 <- wrap.plot(r.wineFortRosen, series = 2, ylab = "Rose, norm")
trel2 <- wrap.plot(r.wineFortRosen, series = 1, ylab = "Fort, norm")
trel3 <- wrap.plot(r.wineFortRose, series = 2, ylab = "Rose")
trel4 <- wrap.plot(r.wineFortRose, series = 1, ylab = "Fort")
plot(trel1, split = c(1, 1, 2, 2), more = TRUE)
plot(trel2, split = c(1, 2, 2, 2), more = TRUE)
plot(trel3, split = c(2, 1, 2, 2), more = TRUE)
plot(trel4, split = c(2, 2, 2, 2))
### End fragment: wineFortRose_scales

### Begin fragment: wineFortRose_forecast
f.wineFortRosen <- rforecast(s.wineFortRosen, groups = list(1:14),
                             len = 13, only.new = TRUE)[, "Rose"]
f.wineFortRosen_long <- c(rep(NA, 174),
                          norm.wineFortRosen["Rose"] * f.wineFortRosen)
xyplot(AustralianWine[100:187, "Rose"] + f.wineFortRosen_long[100:187] ~
         time(AustralianWine)[100:187],
       type = "l", xlab = "Time", ylab = "Rose", lty = c(1, 2))
### End fragment: wineFortRose_forecast

### Begin fragment: winetotal_forecast
FilledRoseAustralianWine <- AustralianWine
FilledRoseAustralianWine[175:176, "Rose"] <- f.wineFortRosen_long[175:176]
mainsales <- ts(rowSums(FilledRoseAustralianWine[, -1]))
total <- FilledRoseAustralianWine[, "Total"]
L = 84
s.totalmain1 <- ssa(list(mainsales[12:187], total[1:176]),
                    L = L, kind = "mssa")
f.totalmain1 <- rforecast(s.totalmain1, groups = list(1:14),
                          len = 11, only.new = TRUE)
s.totalmain2 <- ssa(list(100 * mainsales[12:187], total[1:176]),
                    L = L, kind = "mssa")
f.totalmain2 <- rforecast(s.totalmain2, groups = list(1:14),
                          len = 11, only.new = TRUE)
s.total <- ssa(total[1:176], L = L, kind = "1d-ssa")
f.total <- rforecast(s.total, groups = list(1:14),
                     len = 11, only.new = TRUE)
xtime <- time(AustralianWine)[177:187]
xtime.labels <- paste(month.abb[round(xtime * 12) %% 12 + 1],
                      floor(xtime), sep = ", ")
xyplot(f.total + f.totalmain1[[2]] + f.totalmain2[[2]] +
         mainsales[177:187] ~ xtime,
       type = "l", xlab = "Time", ylab = "Total",
       scales = list(x = list(labels = xtime.labels)),
       auto.key = list(text = c("Separate forecast of 'Total'",
                  "Forecast of 'Total' using 'Main'",
                  "Forecast of 'Total' using 'Main' with weight 100",
                  "Known `Main' sales"),
                       lines = TRUE, points = FALSE))
### End fragment: winetotal_forecast

### Begin fragment: wine_full
L <- 163
norm.wine <- sqrt(colMeans(wine[, -1]^2))
winen <- sweep(wine[, -1], 2, norm.wine, "/")
s.winen <- ssa(winen, L = L, kind = "mssa")
r.winen <- reconstruct(s.winen,
                       groups = list(Trend = c(1, 2, 5),
                                     Seasonality = c(3:4, 6:12)))
plot(r.winen, add.residuals = FALSE,
     plot.method = "xyplot",
     slice = list(component = 1), screens = list(colnames(winen)),
     col = c("blue", "green", "red", "violet", "black", "green4"),
     lty = rep(c(1, 2), each = 6),
     scales = list(y = list(draw = FALSE)),
     layout = c(1, 6))
plot(r.winen, plot.method = "xyplot", add.original = FALSE,
     add.residuals = FALSE, slice = list(component = 2),
     col = c("blue", "green", "red", "violet", "black", "green4"),
     scales = list(y = list(draw = FALSE)),
     layout = c(1, 6))
### End fragment: wine_full

set.seed(1)
### Begin fragment: mssa_comparison
N <- 71
sigma <- 5
Ls <- c(12, 24, 36, 48, 60)
len <- 24
signal1 <- 30 * cos(2*pi * (1:(N + len)) / 12)
signal2 <- 30 * cos(2*pi * (1:(N + len)) / 12 + pi / 4)
signal <- cbind(signal1, signal2)
R <- 10
mssa.errors <- function(Ls) {
  f1 <- signal1[1:N] + rnorm(N, sd = sigma)
  f2 <- signal2[1:N] + rnorm(N, sd = sigma)
  f <- cbind(f1, f2)
  err.rec <- numeric(length(Ls)); names(err.rec) <- Ls
  err.for <- numeric(length(Ls)); names(err.for) <- Ls
  for (l in seq_along(Ls)) {
    L <- Ls[l]
    s <- ssa(f, L = L, kind = "mssa")
    rec <- reconstruct(s, groups = list(1:2))[[1]]
    err.rec[l] <- mean((rec - signal[1:N, ])^2)
    pred <- vforecast(s, groups = list(1:2), direction = "row",
                      len = len, drop = TRUE)
    err.for[l] <- mean((pred - signal[-(1:N), ])^2)
  }
  list(Reconstruction = err.rec, Forecast = err.for)
}
mres <- replicate(R, mssa.errors(Ls))
err.rec <- rowMeans(simplify2array(mres["Reconstruction", ]))
err.for <- rowMeans(simplify2array(mres["Forecast", ]))
err.rec
err.for
signal <- signal1 + 1i*signal2
cssa.errors <- function(Ls) {
  f1 <- signal1[1:N] + rnorm(N, sd = sigma)
  f2 <- signal2[1:N] + rnorm(N, sd = sigma)
  f <- f1 + 1i*f2
  err.rec <- numeric(length(Ls)); names(err.rec) <- Ls
  err.for <- numeric(length(Ls)); names(err.for) <- Ls
  for (l in seq_along(Ls)) {
    L <- Ls[l]
    s <- ssa(f, L = L, kind = "cssa", svd.method = "svd")
    rec <- reconstruct(s, groups = list(1:2))[[1]]
    err.rec[l] <- mean(abs(rec - signal[1:N])^2)
    pred <- vforecast(s, groups = list(1:2), len = len,
                      drop = TRUE)
    err.for[l] <- mean(abs(pred - signal[-(1:N)])^2)
  }
  list(Reconstruction = err.rec, Forecast = err.for)
}
cres <- replicate(R, cssa.errors(Ls))
err.rec <- rowMeans(simplify2array(cres["Reconstruction", ]))
err.for <- rowMeans(simplify2array(cres["Forecast", ]))
err.rec
err.for
### End fragment: mssa_comparison

### Begin fragment: Mars_input
library("lattice")
data("Mars", package = "Rssa")
### End fragment: Mars_input

### Begin fragment: Mars_25_dec_svd
# ssa(Mars, kind = "2d-ssa", L = c(25, 25), svd.method = "svd")
### End fragment: Mars_25_dec_svd

### Begin fragment: Mars_25_dec_eigen
system.time(ssa(Mars, kind = "2d-ssa", L = c(25, 25),
                svd.method = "eigen"))
### End fragment: Mars_25_dec_eigen

### Begin fragment: Mars_25_dec
system.time(s.Mars.25 <- ssa(Mars, kind = "2d-ssa", L = c(25, 25)))
### End fragment: Mars_25_dec

### Begin fragment: Mars_25_rec
r.Mars.25 <- reconstruct(s.Mars.25,
                         groups = list(Noise = c(12, 13, 15, 16)))
plot(r.Mars.25, cuts = 255, layout = c(3, 1))
### End fragment: Mars_25_rec

### Begin fragment: Mars_25_ident
plot(s.Mars.25, type = "vectors", idx = 1:20,
     cuts = 255, layout = c(10, 2),
     plot.contrib = FALSE)
plot(wcor(s.Mars.25, groups = 1:30), scales = list(at = c(10, 20, 30)))
### End fragment: Mars_25_ident

### Begin fragment: Mars_160_80_rec
system.time(s.Mars.160.80 <-
  ssa(Mars, kind = "2d-ssa", L = c(160, 80)))
r.Mars.160.80.groups <- list(Noise = c(36, 37, 42, 43))
r.Mars.160.80 <- reconstruct(s.Mars.160.80, groups = r.Mars.160.80.groups)
plot(r.Mars.160.80, cuts = 255, layout = c(3, 1))
### End fragment: Mars_160_80_rec

### Begin fragment: Mars_160_80_esprit
pe.Mars.160.80 <- parestimate(s.Mars.160.80, groups = r.Mars.160.80.groups)
pe.Mars.160.80
pe.Mars.160.80[[1]]
pe.Mars.160.80[[2]]
plot(pe.Mars.160.80, col = c(11, 12, 13, 14))
plot(s.Mars.160.80, type = "vectors", idx = r.Mars.160.80.groups$Noise,
     cuts = 255, layout = c(4, 1), plot.contrib = FALSE)
### End fragment: Mars_160_80_esprit

### Begin fragment: 2dssa_plots
plot2d <- function(x) {
  regions <- list(col = colorRampPalette(grey(c(0, 1))));
  levelplot(t(x[seq(nrow(x), 1, -1), ]), aspect = "iso",
            par.settings = list(regions = regions), colorkey = FALSE, 
            scales = list(draw = FALSE, relation = "same"), 
            xlab = "", ylab = "")
}
centered.mod.fft <- function(x) {
  N <- dim(x)
  shift.exp <- exp(2i*pi * floor(N/2) / N)
  shift1 <- shift.exp[1]^(0:(N[1] - 1))
  shift2 <- shift.exp[2]^(0:(N[2] - 1))
  Mod(t(mvfft(t(mvfft(outer(shift1, shift2) * x)))))
}
### End fragment: 2dssa_plots

### Begin fragment: 2dssa_brecon_dec
library("raster")
UK <- getData("alt", country = "GB", mask = TRUE)
brecon <- crop(UK, extent(UK, 1040, 1119, 590, 689))
m.brecon <- as.matrix(brecon)
s.brecon <- ssa(m.brecon, kind = "2d-ssa", L = c(8, 8),
                svd.method = "eigen")
plot(s.brecon, type = "vectors", idx = 1:32,
     cuts = 255, layout = c(8, 4), plot.contrib = FALSE)
plot(wcor(s.brecon, groups = 1:32), scales = list(at = c(10, 20, 30)))
### End fragment: 2dssa_brecon_dec

### Begin fragment: 2dssa_brecon_rec
r.brecon <- reconstruct(s.brecon, groups = list(1:3, 4:8, 9:17))
plot(r.brecon, cuts = 255, layout = c(5, 1), 
     par.strip.text = list(cex = 0.75))
plot(r.brecon, cuts = 255, layout = c(5, 1), 
     par.strip.text = list(cex = 0.75), type = "cumsum", at = "free")
brecon.F1 <- raster(r.brecon$F1, template = brecon)
brecon.F2 <- raster(r.brecon$F2, template = brecon)
### End fragment: 2dssa_brecon_rec

### Begin fragment: 2dssa_brecon_rec_dft
library("lattice")
plot2d(centered.mod.fft(m.brecon-r.brecon$F1))
plot2d(centered.mod.fft(m.brecon-r.brecon$F1-r.brecon$F2))
plot2d(centered.mod.fft(m.brecon-r.brecon$F1-r.brecon$F2-r.brecon$F3))
### End fragment: 2dssa_brecon_rec_dft

### Begin fragment: Mars_shaped_dec
mask.Mars.0 <- (Mars != 0)
mask.Mars.1 <- (Mars != 255)
is.na(Mars[!mask.Mars.0]) <- TRUE
system.time(s.Mars.shaped <-
  ssa(Mars, kind = "2d-ssa",
      mask = mask.Mars.1, wmask = circle(15)))
mask.Mars.res <- (s.Mars.shaped$weights > 0)
plot2d(mask.Mars.0)
plot2d(mask.Mars.1)
plot2d(mask.Mars.res)
### End fragment: Mars_shaped_dec

### Begin fragment: Mars_shaped_rec
r.Mars.shaped <- reconstruct(s.Mars.shaped,
                             groups = list(Noise = c(7, 8, 9, 10)))
plot(r.Mars.shaped, cuts = 255, layout = c(3, 1), fill.color = "green")
### End fragment: Mars_shaped_rec

### Begin fragment: Mars_shaped_ident
plot(s.Mars.shaped, type = "vectors", idx = 1:20, fill.color = "green",
     cuts = 255, layout = c(10, 2), plot.contrib = FALSE)
plot(wcor(s.Mars.shaped, groups = 1:30),
     scales = list(at = c(10, 20, 30)))
### End fragment: Mars_shaped_ident

### Begin fragment: Mars_rect_vs_shaped
Mars.sh <- r.Mars.shaped$Noise
Mars.rect.sh <- Mars.rect <- r.Mars.25$Noise
is.na(Mars.rect.sh[is.na(Mars.sh)]) <- TRUE
library("latticeExtra")
p.part.rect <- plot2d(Mars.rect[60:110, 200:250]) +
                 layer(panel.fill(col = "green", alpha = 0.2), under = FALSE) +
                 plot2d(Mars.rect.sh[60:110, 200:250])
p.part.shaped <- plot2d(r.Mars.shaped[[1]][60:110, 200:250]) +
                   layer(panel.fill(col = "green"), under = TRUE)
plot(c(Rectangular = p.part.rect, Shaped = p.part.shaped))
### End fragment: Mars_rect_vs_shaped