Rssa (version 1.0.2)

ssa: Create a new SSA object

Description

Set up the SSA object and perform the decomposition, if necessary.

Usage

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", "rspectra", "primme"),
    force.decompose = TRUE)

Arguments

x

object to be decomposed. See ssa-input for more information

L

integer, window length. Fixed to half of the series length by default. Should be vector of length 2 for 2d SSA

neig

integer, number of desired eigentriples. If 'NULL', then sane default value will be used, see 'Details'

mask

for shaped 2d SSA case only. Logical matrix with same dimension as x. Specifies form of decomposed array. If 'NULL', then all non-NA elements will be used

wmask

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 decompose

kind

SSA method. This includes ordinary 1d SSA, 2d SSA, Toeplitz variant of 1d SSA, multichannel variant of SSA and complex SSA

circular

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 and Golyandina (2014) for more information

svd.method

singular value decomposition method. See 'Details' for more info

column.projector, row.projector

column and row signal subspaces projectors for SSA with projection. See `Details' for information about methods of projectors specification

column.oblique, row.oblique

column and row matrix weights for Weighted Oblique SSA. See `Details' for information about how to use this feature

force.decompose

logical, if 'TRUE' then the decomposition is performed before return.

Value

Object of class `ssa'. The precise layout of the object is mostly meant opaque and subject to change in different version of the package. See ssa-object for details.

Details

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. For the comprehensive description of SSA modifications and their algorithms see Golyandina et al (2018).

Variants of SSA

The following implementations of the SSA method are supported (corresponds to different values of kind argument):

1d-ssa

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

Toeplitz variant of 1d SSA. See Section 1.7.2 in Golyandina et al (2001). This is also known as Vautard-Ghil variant of SSA or VG-SSA for analysis of stationary time series, see Vautard and Ghil (1989).

mssa

Multichannel SSA for simultaneous decomposition of several time series (possible of unequal length). See Golyandina and Stepanov (2005).

cssa

Complex variant of 1d SSA.

2d-ssa

2d SSA for decomposition of images and arrays. See Golyandina and Usevich (2009) and Golyandina et.al (2015) for more information.

nd-ssa

Multidimensional SSA decomposition for arrays (tensors).

Window shape selection (for shaped 2d SSA)

Window shape may be specified by argument wmask. If wmask is 'NULL', then standard rectangular window (specified by L) will be used.

Also in wmask one may use following functions:

circle(R)

circular mask of radius R

triangle(side)

mask in form of isosceles right-angled triangle with cathetus side. Right angle lay on topleft corner of container square matrix

These functions are not exported, they defined only for wmask expression. 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 specification for SSA with projection

Projectors are specified by means of column.projector and 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

both projectors are 'none'

corresponds to ordinary 1D SSA,

column.projector='centering'

corresponds to 1D SSA with centering,

column.projector='centering' and row.projector='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).

Weighted Oblique SSA

Corresponding matrix norm weights may be specified for ordinary 1D SSA case by means of column.oblique and row.oblique arguments. These arguments should be either 'identical' or positive numeric vectors of length L and N - L + 1 for column.oblique and row.oblique respectively.

Weighted Oblique SSA inside Cadzow iterations may improve finite-rank estimation of signal (see e.g. Cadzow(alpha) iterations in Zvonarev and Golyandina (2017) for more information).

SVD methods

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 svd.method) argument:

auto

Automatic method selection depending on the series length, window length, SSA kind and number of eigenvalues requested.

nutrlan

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.

propack

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.

svd

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.

eigen

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.

rspectra

SVD via svds function from Rspectra package (if installed)

primme

SVD via svds function from PRIMME package (if installed)

Usually the 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.

References

Golyandina N., Korobeynikov A., Zhigljavsky A. (2018): Singular Spectrum Analysis with R. Use R!. Springer, Berlin, Heidelberg.

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.

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

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, p.1-78. https://www.jstatsoft.org/article/view/v067i02

Golyandina, N. and Shlemov, A. (2017): Semi-nonparametric singular spectrum analysis with projection. Statistics and its Interface, Vol. 10, Issue 1, p. 47-57. http://arxiv.org/abs/1401.4980

Zvonarev, N. and Golyandina, N. (2017): Iterative algorithms for weighted and unweighted finite-rank time-series approximations. Statistics and its Interface, Vol. 10, Issue 1, p. 5-18. http://arxiv.org/abs/1507.02751

See Also

svd, ssa-object, ssa-input, ssa.capabilities, decompose, reconstruct, plot, forecast,

Examples

# NOT RUN {
# 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)))
# }
# NOT RUN {
plot(r)
# }
# NOT RUN {
# 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
# }
# NOT RUN {
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)))
# }
# NOT RUN {
# 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))
# }
# NOT RUN {
# Plot components, original image and residuals
plot(r)
# }
# NOT RUN {
# Real example: Mars photo
data(Mars)
# Decompose only Mars image (without backgroud)
s <- ssa(Mars, mask = Mars != 0, wmask = circle(50), kind = "2d-ssa")
# }
# NOT RUN {
# 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)))))
# }
# NOT RUN {
# 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)[1]),
               seq_len(dim(Barbara)[2]),
               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)

# }
# NOT RUN {
# 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)
# }