Chapter 4: SSA for multivariate time series

Here you can find the code listings in R language from the corresponding chapter of the book.

You need to load Rssa, ssabook, lattice, latticeExtra, plyr, fma to run these examples.

Contents

  1. Chapter 4: SSA for multivariate time series
    1. Fragment 4.1.1 (‘Stocks’: Reconstruction)
      1. Produced output
    2. Fragment 4.2.1 (‘FORT’ and `DRY’: Reconstruction)
      1. Produced output
    3. Fragment 4.2.2 (‘FORT’ and ‘DRY’: Identification)
      1. Produced output
    4. Fragment 4.3.1 (‘FORT’ and ‘DRY’: Forecast)
      1. Produced output
    5. Fragment 4.3.2 (Simulation for accuracy estimation)
      1. Produced output
    6. Fragment 4.4.1 (‘FORT’ and ‘ROSE’: Influence of series scales)
      1. Produced output
    7. Fragment 4.4.2 (‘FORT’ and ‘ROSE’: Filling-in the missing data in ‘ROSE’)
      1. Produced output
    8. Fragment 4.4.3 (‘Total’ and ‘Mainsales’: Forecast to fill-in ‘Total’)
      1. Produced output
    9. Fragment 4.4.4 (‘Australian wines’: Simultaneous decomposition by MSSA)
      1. Produced output

Fragment 4.1.1 (‘Stocks’: Reconstruction)

s <- ssa(EuStockMarkets[, 1] + 1i*EuStockMarkets[, 2], 
         kind = "cssa", svd.method = "svd")
r <- reconstruct(s, groups = list(Trend = 1:2))
plot(r, plot.method = "xyplot", layout = c(2, 3))
plot(s, type = "vectors", idx = 1:8)
len = 2
print(rforecast(s, groups = list(Trend = 1:2), len = len)[1:len])

Produced output

Reconstruction ‘Stocks’: Reconstructed trends. Eigenvectors ‘Stocks’: Eigenvectors, real and imaginary parts.

[1] 6156.061+8492.425i 6169.006+8507.808i

Fragment 4.2.1 (‘FORT’ and `DRY’: Reconstruction)

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))

Produced output

Reconstruction ‘FORT’ and ‘DRY’: Reconstructed trend and seasonality.

Fragment 4.2.2 (‘FORT’ and ‘DRY’: Identification)

plot(s.wineFortDry, type = "vectors", idx = 1:8)
plot(s.wineFortDry, type = "paired", idx = 2:11, 
     plot.contrib = FALSE)
print(parestimate(s.wineFortDry, groups = list(2:3, 4:5), 
                  method = "esprit"))
plot(wcor(s.wineFortDry, groups = 1:30),
     scales = list(at = c(10, 20, 30)))
si.wineFortDry <- iossa(s.wineFortDry, 
                    nested.groups = list(c(1,6), c(2:5, 7:12)))
plot(si.wineFortDry, type = "vectors", idx = 1:8)

Produced output

Eigenvectors ‘FORT’ and ‘DRY’: 1D graphs of eigenvectors.

2D Scatterplot ‘FORT’ and ‘DRY’: 2D scatterplots of eigenvectors.

Eigenvectors ‘FORT’ and ‘DRY’: 1D graphs of eigenvectors after Iterative O-SSA.

$F1
   period     rate   |    Mod     Arg  |     Re        Im
   12.128  -0.004789 |  0.99522   0.52 |  0.86463   0.49283
  -12.128  -0.004789 |  0.99522  -0.52 |  0.86463  -0.49283

$F2
   period     rate   |    Mod     Arg  |     Re        Im
    4.007  -0.001226 |  0.99877   1.57 |  0.00279   0.99877
   -4.007  -0.001226 |  0.99877  -1.57 |  0.00279  -0.99877

Fragment 4.3.1 (‘FORT’ and ‘DRY’: Forecast)

f.wineFortDry <- rforecast(s.wineFortDry, 
                           groups = list(1, 1:12),
                           len = 60, only.new = TRUE)
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))

Produced output

Forecast ‘FORT’ and ‘DRY’: Forecast of the signal.

Fragment 4.3.2 (Simulation for accuracy estimation)

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", ]))
print(err.rec)
print(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", ]))
print(err.rec)
print(err.for)

Produced output

      12       24       36       48       60 
2.869683 1.587789 1.248881 1.153730 1.855115 
      12       24       36       48       60 
2.671251 2.578059 1.501565 2.595378 4.564218 
      12       24       36       48       60 
7.349316 4.298144 4.101666 4.298144 7.349316 
      12       24       36       48       60 
24.67425 13.60116 14.54819 11.72135 15.86380 

Fragment 4.4.1 (‘FORT’ and ‘ROSE’: Influence of series 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", main = NULL)
trel2 <- wrap.plot(r.wineFortRosen, series = 1, 
                   ylab = "Fort, norm", main = NULL)
trel3 <- wrap.plot(r.wineFortRose, series = 2,
                   ylab = "Rose", main = NULL)
trel4 <- wrap.plot(r.wineFortRose, series = 1,
                   ylab = "Fort", main = NULL)
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))

Produced output

Trends ‘FORT’ and ‘ROSE’: Trends with normalization (ET1,12,14) and without (ET1).

Fragment 4.4.2 (‘FORT’ and ‘ROSE’: Filling-in the missing data in ‘ROSE’)

wineFortRose <- AustralianWine[, c("Fortified", "Rose")]
L <- 84
wineFortRose <- AustralianWine[, c("Fortified", "Rose")]
norm.wineFortRosen <- 
  sqrt(colMeans(wineFortRose^2, na.rm = TRUE))
wineFortRosen <- 
  sweep(wineFortRose, 2, norm.wineFortRosen, "/")
s.wineFortRosen <- ssa(wineFortRosen, L = L, kind = "mssa")
g.wineFortRosen <- 
  gapfill(s.wineFortRosen, groups = list(1:14))
ig.wineFortRosen <- 
  igapfill(s.wineFortRosen, groups = list(1:14))
ig.wineFortRose <- 
  norm.wineFortRosen["Rose"] * ig.wineFortRosen
g.wineFortRose <- 
  norm.wineFortRosen["Rose"] * g.wineFortRosen
xyplot(AustralianWine[100:187, "Rose"] + 
         ig.wineFortRose[100:187, "Rose"] +
         g.wineFortRose[100:187, "Rose"] ~
         time(AustralianWine)[100:187],
       type = "l", xlab = "Time", ylab = "Rose", 
       lty = c(1, 2, 1), lwd = c(2, 1, 1),
       auto.key = list(text = c("`Rose'",
                                "Iterative gap filling",
                                "Subspace-based gap-filling")))

Produced output

Forecast ‘FORT’ and ‘ROSE’: Filling-in of `ROSE’ by two methods.

Fragment 4.4.3 (‘Total’ and ‘Mainsales’: Forecast to fill-in ‘Total’)

FilledRose <- AustralianWine
FilledRose[175:176, "Rose"] <- g.wineFortRose[175:176]
mainsales <- ts(rowSums(FilledRose[, -1]))
tsp(mainsales) <- tsp(AustralianWine)
wine.add.mainsales <- cbind(FilledRose, mainsales)
colnames(wine.add.mainsales) <- 
  c(colnames(FilledRose), "Mainsales")
L <- 84
s.totalmain <- ssa(wine.add.mainsales[, c("Mainsales", 
                                          "Total")],
                   L = L, kind = "mssa")
f.totalmain <- rforecast(s.totalmain, groups = list(1:14),
                         len = 11, only.new = TRUE)
plot(f.totalmain, main = "", xlab = NULL, oma = c(3, 1, 1, 1))

Produced output

Forecast ‘Total’ and ‘Mainsales’: Forecast to fill-in ‘Total’.

Fragment 4.4.4 (‘Australian wines’: Simultaneous decomposition by MSSA)

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))

Produced output

Trends ‘Australian wines’: Extraction of trends.

Trends ‘Australian wines’: Extraction of seasonality.