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
- Chapter 4: SSA for multivariate time series
- Fragment 4.1.1 (‘Stocks’: Reconstruction)
- Fragment 4.2.1 (‘FORT’ and `DRY’: Reconstruction)
- Fragment 4.2.2 (‘FORT’ and ‘DRY’: Identification)
- Fragment 4.3.1 (‘FORT’ and ‘DRY’: Forecast)
- Fragment 4.3.2 (Simulation for accuracy estimation)
- Fragment 4.4.1 (‘FORT’ and ‘ROSE’: Influence of series scales)
- Fragment 4.4.2 (‘FORT’ and ‘ROSE’: Filling-in the missing data in ‘ROSE’)
- Fragment 4.4.3 (‘Total’ and ‘Mainsales’: Forecast to fill-in ‘Total’)
- Fragment 4.4.4 (‘Australian wines’: Simultaneous decomposition by MSSA)
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
‘Stocks’: Reconstructed trends.
‘Stocks’: Eigenvectors, real and imaginary parts.
[1] 6156.061+8492.425i 6169.006+8507.808i
Fragment 4.2.1 (‘FORT’ and `DRY’: Reconstruction)
library("Rssa")
data("AustralianWine", package = "Rssa")
wine <- window(AustralianWine, end = time(AustralianWine)[174])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
‘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
‘FORT’ and ‘DRY’: 1D graphs of eigenvectors.
‘FORT’ and ‘DRY’: 2D scatterplots of 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
‘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
‘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
‘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
‘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
‘Australian wines’: Extraction of trends.
‘Australian wines’: Extraction of seasonality.