Set up the SSA object and perform the decomposition, if necessary.
1 2 3 4 5 6 7 8 9 10 11 12 13 14
ssa(x, L = (N + 1) %/% 2, neig = NULL, mask = NULL, wmask = NULL, column.projector = "none", row.projector = "none", column.oblique = "identity", row.oblique = "identity", ..., kind = c("1d-ssa", "2d-ssa", "nd-ssa", "toeplitz-ssa", "mssa", "cssa"), circular = FALSE, svd.method = c("auto", "nutrlan", "propack", "svd", "eigen"), force.decompose = TRUE)
object to be decomposed. See
integer, window length. Fixed to half of the series length by default. Should be vector of length 2 for 2d SSA
integer, number of desired eigentriples. If 'NULL', then sane default value will be used, see 'Details'
for shaped 2d SSA case only. Logical matrix with same dimension as
for shaped 2d SSA case only. Logical matrix which specifies window form. See ‘Details’ for more information about the window shape selection
further arguments passed to
SSA method. This includes ordinary 1d SSA, 2d SSA, Toeplitz variant of 1d SSA, multichannel variant of SSA and complex SSA
logical vector of one or two elements, describes series topology for 1d SSA and Toeplitz SSA or field topology for 2d SSA. 'TRUE' means series circularity for 1d case or circularity by a corresponding coordinate for 2d case. See (Shlemov, 2014) for more info
singular value decomposition method. See 'Details' for more info
column and row signal subspaces projectors for SSA with projection. See ‘Details’ for information about methods of projectors specification
column and row matrix weights for Weighted Oblique SSA. See ‘Details’ for information about how to use this feature
logical, if 'TRUE' then the decomposition is performed before return.
This is the main entry point to the package. This routine constructs the SSA object filling all necessary internal structures and performing the decomposition if necessary.
The following implementations of the SSA method are supported
(corresponds to different values of
Basic 1d SSA as described in Chapter 1 of (Golyandina et al, 2001). This is also known as Broomhead-King variant of SSA or BK-SSA, see (Broomhead and King, 1986).
Toeplitz variant of 1d SSA. See Section 1.7.2 in (Golyandina et al, 2001). This is also knows as Vatuard-Gill variant of SSA or VG-SSA for analysis of stationary time series, see (Vautard and Ghil, 1989).
Multichannel SSA for simultaneous decomposition of several time series (possible of unequal length). See (Golyandina and Stepanov, 2005).
Complex variant of 1d SSA.
2d SSA for decomposition of images and arrays. See (Golyandina and Usevich, 2009, and Golyandina et.al, 2015) for more information.
Multidimensional SSA decomposition for arrays (tensors).
Window shape may be specified by argument
wmask is 'NULL',
then standard rectangular window (specified by
L) will be used.
wmask one may use following functions:
circular mask of radius
mask in form of isosceles right-angled triangle with
side. Right angle lay on topleft corner of container square
These functions are not exported, they defined only for
If one has objects with the same names and wants to use them rather than these functions,
one should use special wrapper function
I() (see 'Examples').
Projectors are specified by means of
row.projector arguments (see Golyandina and Shlemov, 2017).
Each may be a matrix of orthonormal
(otherwise QR orthonormalization process will be perfomed) basis of
projection subspace, or single integer, which will be interpreted as
dimension of orthogonal polynomial basis (note that the dimension
equals to degree plus 1, e.g. quadratic basis has dimension 3), or
one of following character strings (or unique prefix): 'none',
'constant' (or 'centering'), 'linear', 'quadratic' or 'qubic' for
orthonormal bases of the corresponding functions.
Here is the the list of the most used options
corresponds to ordinary 1D SSA,
corresponds to 1D SSA with centering,
corresponds to 1D SSA with double centering.
SSA with centering and double centering may improve the separation of linear trend (see (Golyandina et.al, 2001) for more information).
Corresponding matrix norm weights may be specified for ordinary 1D SSA case
by means of
row.oblique arguments. These
arguments should be either 'identical' or positive numeric vectors of length
N - L + 1 for
Weighted Oblique SSA inside
may improve finite-rank estimation of signal
(see e.g. Cadzow(alpha) iterations in (Zvonarev and Golyandina, 2017)
for more information).
The main step of the SSA method is the singular decomposition of the
so-called series trajectory matrix. Package provides several
implementations of this procedure (corresponds to different values of
Automatic method selection depending on the series length, window length, SSA kind and number of eigenvalues requested.
Thick-restart Lanczos eigensolver which operates on cross-product matrix. This methods exploits the Hankel structure of the trajectory matrix efficiently and is really fast. The method allows the truncated SVD (only specifid amount of eigentriples to be computed) and the continuation of the decomposition. See (Korobeynikov, 2010) for more information.
SVD via implicitly restarted Lanczos bidiagonalization with partial reothogonalization. This methods exploits the Hankel structure of the trajectory matrix efficiently and is really fast. This is the 'proper' SVD implementation (the matrix of factor vectors are calculated), thus the memory requirements of the methods are higher than for nu-TRLAN. Usually the method is slightly faster that nu-TRLAN and more numerically stable. The method allows the truncated SVD (only specifid amount of eigentriples to be computed). See (Korobeynikov, 2010) for more information.
Full SVD as provided by LAPACK DGESDD routine. Neither continuation of the decomposition nor the truncated SVD is supported. The method does not assume anything special about the trajectory matrix and thus is slow.
Full SVD via eigendecompsition of the cross-product matrix. In many cases faster than previous method, but still really slow for more or less non-trivial matrix sizes.
ssa function tries to provide the best SVD
implementation for given series length and the window size. In
particular, for small series and window sizes it is better to use
generic black-box routines (as provided by 'svd' and 'eigen'
methods). For long series special-purpose routines are to be used.
Object of class ‘ssa’. The precise layout of the object is mostly
meant opaque and subject to change in different version of the
ssa-object for details.
Broomhead, D.S., and King, G.P. (1986a): Extracting qualitative dynamics from experimental data, Physica D, 20, 217–236.
Vautard, R., and Ghil, M. (1989): Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series, Physica D, 35, 395–424.
Golyandina, N., Nekrutkin, V. and Zhigljavsky, A. (2001): Analysis of Time Series Structure: SSA and related techniques. Chapman and Hall/CRC. ISBN 1584881941
Golyandina, N. and Stepanov, D. (2005): SSA-based approaches to analysis and forecast of multidimensional time series. In Proceedings of the 5th St.Petersburg Workshop on Simulation, June 26-July 2, 2005, St. Petersburg State University, St. Petersburg, 293–298. http://www.gistatgroup.com/gus/mssa2.pdf
Golyandina, N. and Usevich, K. (2009): 2D-extensions of singular spectrum analysis: algorithm and elements of theory. In Matrix Methods: Theory, Algorithms, Applications. World Scientific Publishing, 450-474.
Korobeynikov, A. (2010): Computation- and space-efficient implementation of SSA. Statistics and Its Interface, Vol. 3, No. 3, Pp. 257-268
Golyandina, N., Korobeynikov, A. (2012, 2014): Basic Singular Spectrum Analysis and Forecasting with R. Computational Statistics and Data Analysis, Vol. 71, Pp. 934-954. http://arxiv.org/abs/1206.6910
Golyandina, N., Zhigljavsky, A. (2013): Singular Spectrum Analysis for time series. Springer Briefs in Statistics. Springer.
Golyandina, N., Korobeynikov, A., Shlemov, A. and Usevich, K. (2015): Multivariate and 2D Extensions of Singular Spectrum Analysis with the Rssa Package. Journal of Statistical Software, Vol. 67, Issue 2. https://www.jstatsoft.org/article/view/v067i02
Golyandina, N. and Shlemov, A. (2017): Semi-nonparametric singular spectrum analysis with projection. Statistics and its Interface. http://arxiv.org/abs/1401.4980
Shlemov, A. and Golyandina, N. (2014): Shaped extensions of singular spectrum analysis. 21st International Symposium on Mathematical Theory of Networks and Systems, July 7-11, 2014. Groningen, The Netherlands. p.1813-1820. http://arxiv.org/abs/1507.05286
Zvonarev, N. and Golyandina, N. (2017): Iterative algorithms for weighted and unweighted finite-rank time-series approximations. Statistics and its Interface. http://arxiv.org/abs/1507.02751
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97
# Decompose 'co2' series with default parameters s <- ssa(co2) # Show the summary summary(s) # Reconstruct the series, with suitable grouping r <- reconstruct(s, groups = list(c(1, 4), c(2, 3), c(5, 6))) plot(r) # Decompose 'EuStockMarkets' series with default parameters s <- ssa(EuStockMarkets, kind = "mssa") r <- reconstruct(s, groups = list(Trend = 1:2)) # Plot original series, trend and residuals superimposed plot(r, plot.method = "xyplot", superpose = TRUE, auto.key = list(columns = 3), col = c("blue", "green", "red", "violet"), lty = c(rep(1, 4), rep(2, 4), rep(3, 4))) # Artificial image for 2dSSA mx <- outer(1:50, 1:50, function(i, j) sin(2*pi * i/17) * cos(2*pi * j/7) + exp(i/25 - j/20)) + rnorm(50^2, sd = 0.1) # Decompose 'mx' with circular window s <- ssa(mx, kind = "2d-ssa", wmask = circle(5), neig = 10) # Reconstruct r <- reconstruct(s, groups = list(1, 2:5)) # Plot components, original image and residuals plot(r) # Real example: Mars photo data(Mars) # Decompose only Mars image (without backgroud) s <- ssa(Mars, mask = Mars != 0, wmask = circle(50), kind = "2d-ssa") # Plot eigenarrays plot(s, type = "vectors", idx = 1:25) # Plot factor arrays plot(s, type = "vectors", vectors = "factor", idx = 1:25) # Reconstruct and plot trend plot(reconstruct(s, 1), fill.uncovered = "original") # Reconstruct and plot texture pattern plot(reconstruct(s, groups = list(c(13,14, 17, 18)))) # I()-wrapper demo circle <- 50 s <- ssa(Mars, wmask = circle(R = I(circle))) # CSSA-based trend extraction s <- ssa(EuStockMarkets[, 1] + 1.0i*EuStockMarkets[, 2], kind = "cssa") r <- reconstruct(s, groups = list(Trend = 1:2)) plot(r) # `co2' decomposition with double projection to linear functions s <- ssa(co2, column.projector = "centering", row.projector = "centering") plot(reconstruct(s, groups = list(trend = seq_len(nspecial(s))))) # Artificial 2d example with double projection ii <- matrix(1:100, 100, 100); jj <- t(ii) x <- ii + 2 * jj s <- ssa(x, column.projector = "centering", row.projector = "centering") plot(s) plot(reconstruct(s, groups = list(trend = seq_len(nspecial(s))))) # 3D-SSA example (2D-MSSA) data(Barbara) Barbara.noised <- Barbara # Corrupt image by regular noise noise <- outer(seq_len(dim(Barbara)), seq_len(dim(Barbara)), function(i, j) sin(2*pi * (i/13 + j/23))) Barbara.noised[,, 1] <- Barbara.noised[,, 1] + 10 * noise Barbara.noised[,, 2] <- Barbara.noised[,, 2] + 30 * noise Barbara.noised[,, 3] <- Barbara.noised[,, 3] + 5 * noise # Normalize image for plotting Barbara.noised <- (Barbara.noised - min(Barbara.noised)) / diff(range(Barbara.noised)) plot(c(0, 1), c(0, 1), type = "n", xlab = "", ylab = "") rasterImage(Barbara.noised, 0, 0, 1, 1, interpolate = FALSE) # Multichannel 2D-SSA ss <- ssa(Barbara.noised, L = c(50, 50, 1)) plot(ss, type = "series", groups = 1:18, slice = list(k = 1)) plot(ss, type = "vectors", idx = 1:12, slice = list(k = 1)) plot(ss, type = "vectors", vectors = "factor", idx = 1:12, slice = list(k = 3)) # Denoise image Barbara.rec <- residuals(ss, groups = 5:6) plot(c(0, 1), c(0, 1), type = "n", xlab = "", ylab = "") rasterImage(Barbara.rec, 0, 0, 1, 1, interpolate = FALSE)