Chapter 2: SSA analysis of one-dimensional time series, Sections 2.1 — 2.7
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 2: SSA analysis of one-dimensional time series, Sections 2.1 — 2.7
- Fragments 2.1.1 (‘Australian Wines’: Input) and 2.1.2 (‘FORT’: Reconstruction)
- Fragment 2.1.3 (‘FORT’: Identification)
- Fragment 2.2.1 (Noisy sinusoid: Toeplitz SSA)
- Fragment 2.2.2 (Simulation: comparison of Toeplitz and Basic SSA)
- Fragment 2.3.1 (‘CO2’: SSA with projection)
- Fragment 2.3.2 (Polynomial trend: SSA with projection)
- Fragment 2.4.1 (Noisy sum of three sinusoids: Iterative O-SSA)
- Fragment 2.4.2 (Noisy sum of three sinusoids: Iterative O-SSA, summary)
- Fragment 2.4.3 (Dependence of iossa error on difference between frequencies)
- Fragment 2.5.1 (Separation of two sine waves with equal amplitudes)
- Fragment 2.6.1 (Decomposition for series with a gap)
- Fragment 2.6.2 (Incomplete decomposition for a series with a gap)
- Fragment 2.7.1 (‘White dwarf’: Auto grouping by clustering)
- Fragment 2.7.2 (‘Production’: Auto grouping by frequency analysis)
Fragments 2.1.1 (‘Australian Wines’: Input) and 2.1.2 (‘FORT’: Reconstruction)
library("Rssa")
data("AustralianWine", package = "Rssa")
wine <- window(AustralianWine, end = time(AustralianWine)[174])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))Produced output
‘FORT’: Decomposition.
Fragment 2.1.3 (‘FORT’: Identification)
plot(s.fort, type = "vectors", idx = 1:8)
plot(s.fort, type = "paired", idx = 2:11, plot.contrib = FALSE)
print(parestimate(s.fort, groups = list(2:3, 4:5),
method = "pairs"))
plot(wcor(s.fort, groups = 1:30),
scales = list(at = c(10, 20, 30)))
plot(reconstruct(s.fort, groups = list(G12 = 2:3, G4 = 4:5,
G6 = 6:7, G2.4 = 8:9)),
plot.method = "xyplot", layout = c(2, 2),
add.residuals = FALSE, add.original = FALSE)Produced output
$F1
period rate | Mod Arg | Re Im
11.971 0.000000 | 1.00000 0.52 | 0.86540 0.50109
$F2
period rate | Mod Arg | Re Im
4.005 0.000000 | 1.00000 1.57 | 0.00177 1.00000
‘FORT’: 1D graphs of eigenvectors.
‘FORT’: 2D scatterplots of eigenvectors.
‘FORT’: Weighted correlations.
‘FORT’: Reconstructed sine waves.
Fragment 2.2.1 (Noisy sinusoid: Toeplitz SSA)
N <- 100
sigma <- 0.5
set.seed(1)
F <- sin (2 * pi * (1:N) / 7) + sigma * rnorm(N)
Fcenter <- F - mean(F)
st <- ssa(Fcenter, L = 50, kind = "toeplitz-ssa")
s <- ssa(F, L = 50, kind = "1d-ssa")
p <- plot(s, type = "vectors", idx = 1:4, layout = c(4, 1))
pt <- plot(st, type = "vectors", idx = 1:4, layout = c(4, 1))
plot(pt, split = c(1, 1, 1, 2), more = TRUE)
plot(p, split = c(1, 2, 1, 2), more = FALSE)
pt <- plot(reconstruct(st, groups = list(1:2)),
plot.method = "xyplot", layout = c(3, 1))
p <- plot(reconstruct(s, groups = list(1:2)),
plot.method = "xyplot", layout = c(3, 1))
plot(pt, split = c(1, 1, 1, 2), more = TRUE)
plot(p, split = c(1, 2, 1, 2), more = FALSE)Produced output
Noisy sinusoid: 1D graphs of eigenvectors (top: Toeplitz SSA, bottom: Basic SSA).
Noisy sinusoid: Reconstruction (top: Toeplitz SSA, bottom: Basic SSA).
Fragment 2.2.2 (Simulation: comparison of Toeplitz and Basic SSA)
Warning: this example takes a lot of computational time.
SIMUL <- FALSE
N <- 100
sigma <- 0.5
set.seed(1)
alpha <- seq(0.0, 0.01, 0.001)
L <- 50
Q <- 1000
if (SIMUL) {
RMSE <-
sapply(alpha,
function(a) {
sqrt(rowMeans(replicate(Q, {
S <- exp(a * (1:N)) * sin(2 * pi * (1:N) / 7)
F <- S + sigma * rnorm(N)
Fcenter <- F - mean(F)
st <- ssa(Fcenter, L = L, kind = "toeplitz-ssa")
s <- ssa(F, L = L, kind = "1d-ssa")
rec <- reconstruct(s, groups = list(1:2))$F1
rec.t <- reconstruct(st, groups = list(1:2))$F1
c("1d-ssa" = mean((rec - S)^2),
"toeplitz" = mean((rec.t - S)^2))
})))
})
toeplitz.sim <- as.data.frame(t(RMSE))
} else {
data("toeplitz.sim", package = "ssabook")
}
matplot(alpha, toeplitz.sim, type = "l", ylim = c(0, 0.25))Produced output
Simulation: Reconstruction accuracy
of Toeplitz (dash red line) and Basic SSA (solid black line).
Fragment 2.3.1 (‘CO2’: SSA with projection)
s2 <- ssa(co2, column.projector = "centering",
row.projector = "centering")
plot(reconstruct(s2, groups =
list(Linear.trend = seq_len(nspecial(s2)))),
add.residuals = FALSE, plot.method = "matplot")
s4 <- ssa(co2, column.projector = 2, row.projector = 2)
plot(reconstruct(s4, groups =
list(Trend = seq_len(nspecial(s4)))),
add.residuals = FALSE, plot.method = "matplot")
plot(s4, type = "vectors", idx = 1:12)
r <- reconstruct(s4,
groups =
list(Signal = c(seq_len(nspecial(s4)), 5:8)))
plot(r, plot.method = "xyplot")Produced output
‘CO2’: Reconstruction of linear trend.
‘CO2’: Reconstruction of the cubic trend.
‘CO2’: 1D graphs of eigenvectors.
‘CO2’: Reconstruction of signal.
Fragment 2.3.2 (Polynomial trend: SSA with projection)
N <- 199
tt <- (1:N) / N
r <- 5
F0 <- 10 * (tt - 0.5)^r
F <- F0 + sin(2 * pi * (1:N) / 10)
L <- 100
dec <- ssa(F, L = L, column.projector = 3, row.projector = 3)
rec1 <- reconstruct(dec, groups =
list(Trend = seq_len(nspecial(dec))))
fit1 <- rec1$Trend
fit1_3b <- lm(fit1 ~ poly((1:N), r, raw = TRUE))
fit3b <- lm(F ~ poly((1:N), r, raw = TRUE))
li <- 1:199
d <- data.frame(Initial = F[li],
dproj = fit1[li],
dproj_reg = predict(fit1_3b)[li],
regr = predict(fit3b)[li], trend = F0[li])
xyplot(as.formula(paste(paste(colnames(d), collapse = "+"),
"~", "1:nrow(d)")),
data = d,
type = "l", ylab = "", xlab = "",
lty = c(1, 1, 1, 1, 1), lwd = c(1, 2, 2, 2, 2),
auto.key = list(columns = 3,
lines = TRUE, points = FALSE))Produced output
Polynomial trend: Comparison of trend reconstructions.
Fragment 2.4.1 (Noisy sum of three sinusoids: Iterative O-SSA)
N <- 100
L <- 50
omega1 <- 0.07
omega2 <- 0.065
omega3 <- 0.15
sigma <- 0.1
set.seed(3)
F <- 2 * sin(2 * pi * omega1 * (1:N)) +
sin(2 * pi * omega2 * (1:N)) +
3 * sin(2 * pi * omega3 * (1:N)) + sigma * rnorm(N)
xyplot(F ~ 1:N, type = "l")
s <- ssa(F, L)
plot(s, type = "vectors", idx = 1:8, layout = c(4, 2))
plot(wcor(s, groups = 1:20), scales = list(at = seq(1,20,2)))
ios <- iossa(s, nested.groups = list(3:4, 5:6), maxiter = 1000)
plot(ios, type = "vectors", idx = 1:8, layout = c(4, 2))
ior <- reconstruct(ios, groups = c(list(1:2), ios$iossa.groups))
plot(ior, plot.method = "xyplot", add.original = FALSE,
add.residuals = FALSE)Produced output
Noisy sum of three sinusoids: The original series.
Noisy sum of three sinusoids, Basic SSA: w-Correlation matrix.
Noisy sum of three sinusoids, Basic SSA: Eigenvectors.
Noisy sum of three sinusoids, Iterative O-SSA: Eigenvectors.
Noisy sum of three sinusoids, Iterative O-SSA: Reconstruction.
Fragment 2.4.2 (Noisy sum of three sinusoids: Iterative O-SSA, summary)
print(ios$iossa.groups)
summary(ios)Produced output
[[1]]
[1] 3 4
[[2]]
[1] 5 6
Call:
iossa.ssa(x = s, nested.groups = list(3:4, 5:6), maxiter = 1000)
Series length: 100, Window length: 50, SVD method: eigen
Special triples: 0
Computed:
Eigenvalues: 50, Eigenvectors: 50, Factor vectors: 6
Precached: 0 elementary series (0 MiB)
Overall memory consumption (estimate): 0.0352 MiB
Iterative O-SSA result:
Converged: yes
Iterations: 243
Initial mean(tau): 0.1032
Initial tau: 0.0007976, 0.2055299
I. O-SSA mean(tau): 0.0004452
I. O-SSA tau: 0.0006709, 0.0002196
Initial max wcor: 0.02442
I. O-SSA max wcor: 0.06986
I. O-SSA max owcor: 0.0732
Fragment 2.4.3 (Dependence of iossa error on difference between frequencies)
Warning: this example takes a lot of computational time.
rowMeansQuantile <- function(x, level = 0.05) {
apply(x, 1,
function(x) {
q <- quantile(x, c(level / 2, 1 - level / 2))
x[x < q[1]] <- q[1]
x[x > q[2]] <- q[2]
mean(x)
})
}
data("iossa.err", package = "ssabook")
lseq <- c(seq(0.03, 0.058, 0.002), seq(0.062, 0.1, 0.002))
iter.real <- rowMeansQuantile(iossa.err$iter.real)
iter.est <- iossa.err$iter.est
err1 <- sqrt(rowMeansQuantile(iossa.err$err1))
err2 <- sqrt(rowMeansQuantile(iossa.err$err2))
xlab <- expression(omega[1])
ylab1 <- expression(N[plain(iter)])
ylab2 <- expression(RMSE)
p1 <- xyplot(iter.real + iter.est ~ lseq,
type = "l", ylab = ylab1, xlab = xlab)
p2 <- xyplot(err1 + err2 ~ lseq,
type = "l", ylab = ylab2, xlab = xlab)
print(p1, split = c(1, 1, 1, 2), more = TRUE)
print(p2, split = c(1, 2, 1, 2), more = FALSE)Produced output
Dependence of number of iterations (top) and RMSE errors of
frequency estimations (bottom) on 𝝎1 for 𝝎2=0.6.
Fragment 2.5.1 (Separation of two sine waves with equal amplitudes)
N <- 100
L <- 50
omega1 <- 0.03
omega2 <- 0.06
sigma <- 0.1
set.seed(3)
F <- sin(2 * pi * omega1 * (1:N)) +
sin(2 * pi * omega2 * (1:N)) +
sigma * rnorm(N)
s <- ssa(F, L = L, neig = min(L, N - L + 1)) #full decomposition
plot(s)
p1 <- plot(s, type = "vectors", idx = 1:4, layout = c(4, 1),
main = "Eigenvectors, Basic SSA")
fos <- fossa(s, nested.groups = list(1:2, 3:4), gamma = 10,
normalize = FALSE)
#the total percent is equal to 100%
print(sum(fos$sigma^2) / sum(s$sigma^2) * 100)
p2 <- plot(fos, type = "vectors", idx = 1:4, layout = c(4, 1),
main = "Eigenvectors, SSA with derivatives")
ios1 <- iossa(s, nested.groups = list(1:2, 3:4), maxiter = 1)
#the total percent is not equal to 100%
print(sum(ios1$sigma^2) / sum(s$sigma^2) * 100)
p3 <- plot(ios1, type = "vectors", idx = 1:4, layout = c(4, 1),
main = "Eigenvectors, Iterative O-SSA, 1 iter")
ios2 <- iossa(ios1, nested.groups = list(1:2, 3:4), maxiter = 1)
#the total percent is not equal to 100%
print(sum(ios2$sigma^2) / sum(s$sigma^2) * 100)
p4 <- plot(ios2, type = "vectors", idx = 1:4, layout = c(4, 1),
main = "Eigenvectors, Iterative O-SSA, 2 iter")
plot(p1, split = c(1, 1, 1, 4), more = TRUE)
plot(p2, split = c(1, 2, 1, 4), more = TRUE)
plot(p3, split = c(1, 3, 1, 4), more = TRUE)
plot(p4, split = c(1, 4, 1, 4), more = FALSE)
fo.rec <- reconstruct(fos, groups = list(1:2, 3:4))
pr1 <- plot(fo.rec, plot.method = "xyplot",
main = "SSA with derivatives", xlab = "")
io.rec <- reconstruct(ios2, groups = ios2$iossa.groups)
pr2 <- plot(io.rec, plot.method = "xyplot",
main = "Iterative O-SSA", xlab = "")
plot(pr1, split = c(1, 1, 1, 2), more = TRUE)
plot(pr2, split = c(1, 2, 1, 2), more = FALSE)Produced output
Noisy sum of sinusoids: Graph of eigenvalues for Basic SSA.
Noisy sum of sinusoids: 1D graphs of eigenvectors for Basic SSA (top),
DerivSSA (middle) and Iterative O-SSA, 1 iteration (second from bottom) and 2 iterations (bottom).
Noisy sum of sinusoids: Reconstructions for DerivSSA (top) and
Iterative O-SSA, 2 iterations (bottom).
[1] 100
[1] 99.62939
[1] 101.7544
In .contribution(x, idx, ...):
Elementary matrices are not F-orthogonal (max F-cor is -0.016).
Contributions can be irrelevant
Fragment 2.6.1 (Decomposition for series with a gap)
F <- co2; F[100:200] <- NA
#prompt for the choice of window length
clplot(F)
# Perform shaped SSA
s1 <- ssa(F, L = 72)
plot(s1, type = "vectors", idx = 1:12)
plot(s1, type = "series", groups = 1:6, layout = c(2, 3))
plot(wcor(s1, groups = 1:20), scales = list(at = seq(1,20,2)))
plot(reconstruct(s1, groups = list(c(1, 4, 7))),
add.residuals = FALSE,
plot.method = "xyplot", superpose = TRUE)Produced output
‘CO2’ with gaps, Shaped SSA: Dependence of proportion of complete vectors on window length.
‘CO2’ with gaps, Shaped SSA: Eigenvectors, L=72.
‘CO2’ with gaps, Shaped SSA: Elementary reconstructed series, L=72.
‘CO2’ with gaps, Shaped SSA: w-Correlation matrix, L=72.
Fragment 2.6.2 (Incomplete decomposition for a series with a gap)
s2 <- ssa(F, L = 120)
#plot(s2, type = "vectors")
#plot(wcor(s2, groups = 1:20))
#plot of reconstruction
plot(reconstruct(s2, groups = list(c(1, 6, 7))),
add.residuals = FALSE,
plot.method = "xyplot", superpose = TRUE)Produced output
‘CO2’ with gaps, Shaped SSA: Trend reconstruction, L=72.
‘CO2’ with gaps, Shaped SSA: Incomplete trend reconstruction, L=120.
Fragment 2.7.1 (‘White dwarf’: Auto grouping by clustering)
data("dwarfst", package = "ssabook")
s <- ssa(dwarfst, L = 100)
g <- grouping.auto(s, grouping.method = "wcor",
method = "average", nclust = 2)
print(g[[1]])
plot(wcor(s, groups = 1:30), scales = list(at = c(1, 11, 30)))
plot(reconstruct(s, groups = g),
add.residuals = FALSE,
plot.method = "xyplot", superpose = FALSE)Produced output
‘White dwarf’: w-Correlation matrix, L=100.
‘White dwarf’: Decomposition with automatic
grouping performed by clustering.
[1] 1 2 3 4 5 6 7 8 9 10 11
Fragment 2.7.2 (‘Production’: Auto grouping by frequency analysis)
data("oilproduction", package = "ssabook")
s <- ssa(oilproduction, L = 120)
plot(s, type = "vectors", vectors = "factor", idx = 1:12)
g0 <- grouping.auto(s, base = "series",
freq.bins = list(Tendency = 1/240,
Trend = 1/24),
threshold = 0.1)
plot(g0, order = TRUE, type = "b")
contrib <- attr(g0, "contributions")[, 2]
print(thr <- sort(contrib, decreasing = TRUE)[9])
g <- grouping.auto(s, base = "series",
freq.bins = list(Tendency = 1/240,
Trend = 1/24),
threshold = thr)
print(g[[1]])
print(g[[2]])
plot(reconstruct(s, groups = g),
add.residuals = FALSE,
plot.method = "xyplot", superpose = TRUE)Produced output
‘Production’: Ordered frequency contributions of factor vectors, L=120.
‘Production’: Factor vectors, L=120.
‘Production’: Two extracted trends of different resolution,
automatic grouping by frequencies.
8
0.861955
[1] 1 2
[1] 1 2 3 6 8 11 12 17 18