Create a new SSA object
Description
Set up the SSA object and perform the decomposition, if necessary.
Usage
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("1dssa", "2dssa", "ndssa", "toeplitzssa", "mssa", "cssa"),
circular = FALSE,
svd.method = c("auto", "nutrlan", "propack", "svd", "eigen"),
force.decompose = TRUE)

Arguments
x 
object to be decomposed. See

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

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 
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, 2014) for more info 
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. 
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.
Variants of SSA
The following implementations of the SSA method are supported
(corresponds to different values of kind
argument):
 1dssa
Basic 1d SSA as described in Chapter 1 of (Golyandina et al, 2001). This is also known as BroomheadKing variant of SSA or BKSSA, see (Broomhead and King, 1986).
 toeplitzssa
Toeplitz variant of 1d SSA. See Section 1.7.2 in (Golyandina et al, 2001). This is also knows as VatuardGill variant of SSA or VGSSA 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.
 2dssa
2d SSA for decomposition of images and arrays. See (Golyandina and Usevich, 2009, and Golyandina et.al, 2015) for more information.
 ndssa
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 rightangled 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 finiterank 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
socalled 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
Thickrestart Lanczos eigensolver which operates on crossproduct 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 nuTRLAN. Usually the method is slightly faster that nuTRLAN 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 crossproduct matrix. In many cases faster than previous method, but still really slow for more or less nontrivial matrix sizes.
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 blackbox routines (as provided by 'svd' and 'eigen'
methods). For long series specialpurpose routines are to be used.
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 ssaobject
for details.
References
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): SSAbased approaches to analysis and forecast of multidimensional time series. In Proceedings of the 5th St.Petersburg Workshop on Simulation, June 26July 2, 2005, St. Petersburg State University, St. Petersburg, 293–298. http://www.gistatgroup.com/gus/mssa2.pdf
Golyandina, N. and Usevich, K. (2009): 2Dextensions of singular spectrum analysis: algorithm and elements of theory. In Matrix Methods: Theory, Algorithms, Applications. World Scientific Publishing, 450474.
Korobeynikov, A. (2010): Computation and spaceefficient implementation of SSA. Statistics and Its Interface, Vol. 3, No. 3, Pp. 257268
Golyandina, N., Korobeynikov, A. (2012, 2014): Basic Singular Spectrum Analysis and Forecasting with R. Computational Statistics and Data Analysis, Vol. 71, Pp. 934954. 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): Seminonparametric 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 711, 2014. Groningen, The Netherlands. p.18131820. http://arxiv.org/abs/1507.05286
Zvonarev, N. and Golyandina, N. (2017): Iterative algorithms for weighted and unweighted finiterank timeseries approximations. Statistics and its Interface. http://arxiv.org/abs/1507.02751
See Also
svd
,
ssaobject
,
ssainput
,
decompose
,
reconstruct
,
plot
,
forecast
,
Examples
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 = "2dssa", 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 = "2dssa")
# 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)))
# CSSAbased 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)))))
# 3DSSA example (2DMSSA)
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)
# Multichannel 2DSSA
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)
