Singular spectrum analysis

The origins of SSA and, more generally, of subspace-based methods for signal processing, go back to the eighteenth century (Prony's method). A key development was the formulation of the spectral decomposition of the covariance operator of stochastic processes by Kari Karhunen and Michel Loève in the late 1940s (Loève, 1945; Karhunen, 1947). Broomhead and King (1986a, b) and Fraedrich (1986) proposed to use SSA and multichannel SSA (M-SSA) in the context of nonlinear dynamics for the purpose of reconstructing the attractor of a system from measured time series. These authors provided an extension and a more robust application of the idea of reconstructing dynamics from a single time series based on the embedding theorem. Several other authors had already applied simple versions of M-SSA to meteorological and ecological data sets (Colebrook, 1978; Barnett and Hasselmann, 1979; Weare and Nasstrom, 1982). Ghil, Vautard and their colleagues (Vautard and Ghil, 1989; Ghil and Vautard, 1991; Vautard et al., 1992; Ghil et al., 2002) noticed the analogy between the trajectory matrix of Broomhead and King, on the one hand, and the Karhunen–Loeve decomposition (Principal component analysis in the time domain), on the other. Thus, SSA can be used as a time-and-frequency domain method for time series analysis — independently from attractor reconstruction and including cases in which the latter may fail. The survey paper of Ghil et al. (2002) is the basis of the section of this article. A crucial result of the work of these authors is that SSA can robustly recover the "skeleton" of an attractor, including in the presence of noise. This skeleton is formed by the least unstable periodic orbits, which can be identified in the eigenvalue spectra of SSA and M-SSA. The identification and detailed description of these orbits can provide highly useful pointers to the underlying nonlinear dynamics. The so-called ‘Caterpillar’ methodology is a version of SSA that was developed in the former Soviet Union, independently of the mainstream SSA work in the West. This methodology became known in the rest of the world more recently (Danilov and Zhigljavsky, Eds., 1997; Golyandina et al., 2001; Zhigljavsky, Ed., 2010; Golyandina and Zhigljavsky, 2013; Golyandina et al., 2018). ‘Caterpillar-SSA’ emphasizes the concept of separability, a concept that leads, for example, to specific recommendations concerning the choice of SSA parameters. This method is thoroughly described in of this article. ==Methodology==

In practice, SSA is a nonparametric spectral estimation method based on embedding a time series \{X(t) : t=1,\ldots,N\} in a vector space of dimension M. SSA proceeds by diagonalizing the M\times M lag-covariance matrix {\textbf C}_X of X(t) to obtain spectral information on the time series, assumed to be stationary in the weak sense. The matrix {\textbf C}_X can be estimated directly from the data as a Toeplitz matrix with constant diagonals (Vautard and Ghil, 1989), i.e., its entries c_{ij} depend only on the lag |i-j|: : c_{ij} = \frac{1}{N-|i-j|} \sum_{t=1}^{N-|i-j|} X(t) X(t+|i-j|). An alternative way to compute {\textbf C}_X, is by using the N' \times M "trajectory matrix" {\textbf D} that is formed by M lag-shifted copies of {\it X(t)}, which are N' =N-M+1 long; then : {\textbf C}_X = \frac{1}{N'} {\textbf D}^{\rm t} {\textbf D}. The M eigenvectors {\textbf E}_k of the lag-covariance matrix {\textbf C}_ X are called temporal empirical orthogonal functions (EOFs). The eigenvalues \lambda_k of {\textbf C}_{X} account for the partial variance in the direction {\textbf E}_k and the sum of the eigenvalues, i.e., the trace of {\textbf C}_{X}, gives the total variance of the original time series X(t). The name of the method derives from the singular values \lambda^{1/2}_k of {\textbf C}_{X}. Decomposition and reconstruction Projecting the time series onto each EOF yields the corresponding temporal principal components (PCs) {\textbf A}_k: : A_k(t) = \sum_{j=1}^{M} X(t+j-1) E_k(j). An oscillatory mode is characterized by a pair of nearly equal SSA eigenvalues and associated PCs that are in approximate phase quadrature (Ghil et al., 2002). Such a pair can represent efficiently a nonlinear, anharmonic oscillation. This is due to the fact that a single pair of data-adaptive SSA eigenmodes often will capture better the basic periodicity of an oscillatory mode than methods with fixed basis functions, such as the sines and cosines used in the Fourier transform. The window width M determines the longest periodicity captured by SSA. Signal-to-noise separation can be obtained by merely inspecting the slope break in a "scree diagram" of eigenvalues \lambda_k or singular values \lambda^{1/2}_k vs. k. The point k^* = S at which this break occurs should not be confused with a "dimension" D of the underlying deterministic dynamics (Vautard and Ghil, 1989). A Monte-Carlo test (Allen and Smith, 1996; Allen and Robertson, 1996; Groth and Ghil, 2015) can be applied to ascertain the statistical significance of the oscillatory pairs detected by SSA. The entire time series or parts of it that correspond to trends, oscillatory modes or noise can be reconstructed by using linear combinations of the PCs and EOFs, which provide the reconstructed components (RCs) {\textbf R}_K: : R_{ K}(t) = \frac{1}{M_t} \sum_{k\in {\textit K}} \sum_{j={L_t}}^{U_t} A_k(t-j+1)E_k(j); here K is the set of EOFs on which the reconstruction is based. The values of the normalization factor M_t, as well as of the lower and upper bound of summation L_t and U_t, differ between the central part of the time series and the vicinity of its endpoints (Ghil et al., 2002). Multivariate extension Multi-channel SSA (or M-SSA) is a natural extension of SSA to an L-channel time series of vectors or maps with N data points \{X_{l}(t): l=1,\dots, L; t=1,\dots, N\}. In the meteorological literature, extended EOF (EEOF) analysis is often assumed to be synonymous with M-SSA. The two methods are both extensions of classical principal component analysis (PCA) but they differ in emphasis: EEOF analysis typically utilizes a number L of spatial channels much greater than the number M of temporal lags, thus limiting the temporal and spectral information. In M-SSA, on the other hand, one usually chooses L \leq M. Often M-SSA is applied to a few leading PCs of the spatial data, with M chosen large enough to extract detailed temporal and spectral information from the multivariate time series (Ghil et al., 2002). However, Groth and Ghil (2015) have demonstrated possible negative effects of this variance compression on the detection rate of weak signals when the number L of retained PCs becomes too small. This practice can further affect negatively the judicious reconstruction of the spatio-temporal patterns of such weak signals, and Groth et al. (2016) recommend retaining a maximum number of PCs, i.e., L=N. Groth and Ghil (2011) have demonstrated that a classical M-SSA analysis suffers from a degeneracy problem, namely the EOFs do not separate well between distinct oscillations when the corresponding eigenvalues are similar in size. This problem is a shortcoming of principal component analysis in general, not just of M-SSA in particular. In order to reduce mixture effects and to improve the physical interpretation, Groth and Ghil (2011) have proposed a subsequent VARIMAX rotation of the spatio-temporal EOFs (ST-EOFs) of the M-SSA. To avoid a loss of spectral properties (Plaut and Vautard 1994), they have introduced a slight modification of the common VARIMAX rotation that does take the spatio-temporal structure of ST-EOFs into account. Alternatively, a closed matrix formulation of the algorithm for the simultaneous rotation of the EOFs by iterative SVD decompositions has been proposed (Portes and Aguirre, 2016). M-SSA has two forecasting approaches known as recurrent and vector. The discrepancies between these two approaches are attributable to the organization of the single trajectory matrix {\textbf X} of each series into the block trajectory matrix in the multivariate case. Two trajectory matrices can be organized as either vertical (VMSSA) or horizontal (HMSSA) as was recently introduced in Hassani and Mahmoudvand (2013), and it was shown that these constructions lead to better forecasts. Accordingly, we have four different forecasting algorithms that can be exploited in this version of MSSA (Hassani and Mahmoudvand, 2013). Prediction In this subsection, we focus on phenomena that exhibit a significant oscillatory component: repetition increases understanding and hence confidence in a prediction method that is closely connected with such understanding. Singular spectrum analysis (SSA) and the maximum entropy method (MEM) have been combined to predict a variety of phenomena in meteorology, oceanography and climate dynamics (Ghil et al., 2002, and references therein). First, the “noise” is filtered out by projecting the time series onto a subset of leading EOFs obtained by SSA; the selected subset should include statistically significant, oscillatory modes. Experience shows that this approach works best when the partial variance associated with the pairs of RCs that capture these modes is large (Ghil and Jiang, 1998). The prefiltered RCs are then extrapolated by least-square fitting to an autoregressive model AR[p], whose coefficients give the MEM spectrum of the remaining “signal”. Finally, the extended RCs are used in the SSA reconstruction process to produce the forecast values. The reason why this approach – via SSA prefiltering, AR extrapolation of the RCs, and SSA reconstruction – works better than the customary AR-based prediction is explained by the fact that the individual RCs are narrow-band signals, unlike the original, noisy time series X(t) (Penland et al., 1991; Keppenne and Ghil, 1993). In fact, the optimal order p obtained for the individual RCs is considerably lower than the one given by the standard Akaike information criterion (AIC) or similar ones. Spatio-temporal gap filling The gap-filling version of SSA can be used to analyze data sets that are unevenly sampled or contain missing data (Kondrashov and Ghil, 2006; Kondrashov et al. 2010). For a univariate time series, the SSA gap filling procedure utilizes temporal correlations to fill in the missing points. For a multivariate data set, gap filling by M-SSA takes advantage of both spatial and temporal correlations. In either case: (i) estimates of missing data points are produced iteratively, and are then used to compute a self-consistent lag-covariance matrix {\textbf C}_X and its EOFs {\textbf E}_k; and (ii) cross-validation is used to optimize the window width M and the number of leading SSA modes to fill the gaps with the iteratively estimated "signal," while the noise is discarded. == As a model-free tool ==

The areas where SSA can be applied are very broad: climatology, marine science, geophysics, engineering, image processing, medicine, econometrics among them. Hence different modifications of SSA have been proposed and different methodologies of SSA are used in practical applications such as trend extraction, periodicity detection, seasonal adjustment, smoothing, noise reduction (Golyandina, et al, 2001). Basic SSA SSA can be used as a model-free technique so that it can be applied to arbitrary time series including non-stationary time series. The basic aim of SSA is to decompose the time series into the sum of interpretable components such as trend, periodic components and noise with no a-priori assumptions about the parametric form of these components. Consider a real-valued time series \mathbb{X}=(x_1,\ldots,x_{N}) of length N. Let L \ (1 be some integer called the window length and K=N-L+1. Main algorithm 1st step: Embedding. Form the trajectory matrix of the series \mathbb{X}, which is the L\!\times\! K matrix : \mathbf{X}=[X_1:\ldots:X_K]=(x_{ij})_{i,j=1}^{L,K}= \begin{bmatrix} x_1&x_2&x_3&\ldots&x_{K}\\ x_2&x_3&x_4&\ldots&x_{K+1}\\ x_3&x_4&x_5&\ldots&x_{K+2}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ x_{L}&x_{L+1}&x_{L+2}&\ldots&x_{N}\\ \end{bmatrix} where X_i=(x_{i},\ldots,x_{i+L-1})^\mathrm{T} \; \quad (1\leq i\leq K) are lagged vectors of size L. The matrix \mathbf{X} is a Hankel matrix which means that \mathbf{X} has equal elements x_{ij} on the anti-diagonals i+j =\,{\rm const}. 2nd step: Singular Value Decomposition (SVD). Perform the singular value decomposition (SVD) of the trajectory matrix \mathbf{X}. Set \mathbf{S}=\mathbf{X} \mathbf{X}^\mathrm{T} and denote by \lambda_1, \ldots,\lambda_L the eigenvalues of \mathbf{S} taken in the decreasing order of magnitude (\lambda_1\geq \ldots \geq \lambda_L\geq 0) and by U_1,\ldots,U_L the orthonormal system of the eigenvectors of the matrix \mathbf{S} corresponding to these eigenvalues. Set d= \mathop{\mathrm{rank}} \mathbf{X} =\max\{i,\ \mbox{such that}\ \lambda_i >0\} (note that d=L for a typical real-life series) and V_i=\mathbf{X}^\mathrm{T} U_i/\sqrt{\lambda_i} (i=1,\ldots,d). In this notation, the SVD of the trajectory matrix \mathbf{X} can be written as : \mathbf{X} = \mathbf{X}_1 + \ldots + \mathbf{X}_d, where :\mathbf{X}_i=\sqrt{\lambda_i}U_i V_i^\mathrm{T} are matrices having rank 1; these are called elementary matrices. The collection (\sqrt{\lambda_i},U_i,V_i) will be called the ith eigentriple (abbreviated as ET) of the SVD. Vectors U_i are the left singular vectors of the matrix \mathbf{X}, numbers \sqrt{\lambda_i} are the singular values and provide the singular spectrum of \mathbf{X}; this gives the name to SSA. Vectors \sqrt{\lambda_i}V_i=\mathbf{X}^\mathrm{T} U_i are called vectors of principal components (PCs). 3rd step: Eigentriple grouping. Partition the set of indices \{1,\ldots,d\} into m disjoint subsets I_1,\ldots,I_m. Let I=\{i_1,\ldots,i_p\}. Then the resultant matrix \mathbf{X}_I corresponding to the group I is defined as \mathbf{X}_I=\mathbf{X}_{i_1}+\ldots+\mathbf{X}_{i_p}. The resultant matrices are computed for the groups I=I_1, \ldots, I_m and the grouped SVD expansion of \mathbf{X} can now be written as : \mathbf{X}=\mathbf{X}_{I_1}+\ldots+\mathbf{X}_{I_m}. 4th step: Diagonal averaging. Each matrix \mathbf{X}_{I_j} of the grouped decomposition is hankelized and then the obtained Hankel matrix is transformed into a new series of length N using the one-to-one correspondence between Hankel matrices and time series. Diagonal averaging applied to a resultant matrix \mathbf{X}_{I_k} produces a reconstructed series \widetilde{\mathbb{X}}^{(k)}=(\widetilde{x}^{(k)}_1,\ldots,\widetilde{x}^{(k)}_N). In this way, the initial series x_1,\ldots,x_N is decomposed into a sum of m reconstructed subseries: : x_n = \sum\limits_{k=1}^m \widetilde{x}^{(k)}_n \ \ (n=1,2, \ldots, N). This decomposition is the main result of the SSA algorithm. The decomposition is meaningful if each reconstructed subseries could be classified as a part of either trend or some periodic component or noise. Theory of SSA separability The two main questions which the theory of SSA attempts to answer are: (a) what time series components can be separated by SSA, and (b) how to choose the window length L and make proper grouping for extraction of a desirable component. Many theoretical results can be found in Golyandina et al. (2001, Ch. 1 and 6). Trend (which is defined as a slowly varying component of the time series), periodic components and noise are asymptotically separable as N\rightarrow \infty. In practice N is fixed and one is interested in approximate separability between time series components. A number of indicators of approximate separability can be used, see Golyandina et al. (2001, Ch. 1). The window length L determines the resolution of the method: larger values of L provide more refined decomposition into elementary components and therefore better separability. The window length L determines the longest periodicity captured by SSA. Trends can be extracted by grouping of eigentriples with slowly varying eigenvectors. A sinusoid with frequency smaller than 0.5 produces two approximately equal eigenvalues and two sine-wave eigenvectors with the same frequencies and \pi/2-shifted phases. Separation of two time series components can be considered as extraction of one component in the presence of perturbation by the other component. SSA perturbation theory is developed in Nekrutkin (2010) and Hassani et al. (2011). Forecasting by SSA If for some series \mathbb{X} the SVD step in Basic SSA gives d, then this series is called time series of rank d (Golyandina et al., 2001, Ch.5). The subspace spanned by the d leading eigenvectors is called signal subspace. This subspace is used for estimating the signal parameters in signal processing, e.g. ESPRIT for high-resolution frequency estimation. Also, this subspace determines the linear homogeneous recurrence relation (LRR) governing the series, which can be used for forecasting. Continuation of the series by the LRR is similar to forward linear prediction in signal processing. Let the series be governed by the minimal LRR x_{n}=\sum_{k=1}^d b_k x_{n-k}. Let us choose L>d, U_1,\ldots,U_d be the eigenvectors (left singular vectors of the L-trajectory matrix), which are provided by the SVD step of SSA. Then this series is governed by an LRR x_{n}=\sum_{k=1}^{L-1} a_k x_{n-k}, where (a_{L-1},\ldots,a_1)^\mathrm{T} are expressed through U_1,\ldots,U_d (Golyandina et al., 2001, Ch.5), and can be continued by the same LRR. This provides the basis for SSA recurrent and vector forecasting algorithms (Golyandina et al., 2001, Ch.2). In practice, the signal is corrupted by a perturbation, e.g., by noise, and its subspace is estimated by SSA approximately. Thus, SSA forecasting can be applied for forecasting of a time series component that is approximately governed by an LRR and is approximately separated from the residual. Multivariate extension Multi-channel, Multivariate SSA (or M-SSA) is a natural extension of SSA to for analyzing multivariate time series, where the size of different univariate series does not have to be the same. The trajectory matrix of multi-channel time series consists of linked trajectory matrices of separate times series. The rest of the algorithm is the same as in the univariate case. System of series can be forecasted analogously to SSA recurrent and vector algorithms (Golyandina and Stepanov, 2005). MSSA has many applications. It is especially popular in analyzing and forecasting economic and financial time series with short and long series length (Patterson et al., 2011, Hassani et al., 2012, Hassani and Mahmoudvand, 2013). Other multivariate extension is 2D-SSA that can be applied to two-dimensional data like digital images (Golyandina and Usevich, 2010). The analogue of trajectory matrix is constructed by moving 2D windows of size L_x \times L_y. MSSA and causality A question that frequently arises in time series analysis is whether one economic variable can help in predicting another economic variable. One way to address this question was proposed by Granger (1969), in which he formalized the causality concept. A comprehensive causality test based on MSSA has recently introduced for causality measurement. The test is based on the forecasting accuracy and predictability of the direction of change of the MSSA algorithms (Hassani et al., 2011 and Hassani et al.,2012). MSSA and EMH The MSSA forecasting results can be used in examining the efficient-market hypothesis controversy (EMH). The EMH suggests that the information contained in the price series of an asset is reflected “instantly, fully, and perpetually” in the asset’s current price. Since the price series and the information contained in it are available to all market participants, no one can benefit by attempting to take advantage of the information contained in the price history of an asset by trading in the markets. This is evaluated using two series with different series length in a multivariate system in SSA analysis (Hassani et al. 2010). MSSA, SSA and business cycles Business cycles plays a key role in macroeconomics, and are interest for a variety of players in the economy, including central banks, policy-makers, and financial intermediaries. MSSA-based methods for tracking business cycles have been recently introduced, and have been shown to allow for a reliable assessment of the cyclical position of the economy in real-time (de Carvalho et al., 2012 and de Carvalho and Rua, 2017). MSSA, SSA and unit root SSA's applicability to any kind of stationary or deterministically trending series has been extended to the case of a series with a stochastic trend, also known as a series with a unit root. In Hassani and Thomakos (2010) and Thomakos (2010) the basic theory on the properties and application of SSA in the case of series of a unit root is given, along with several examples. It is shown that SSA in such series produces a special kind of filter, whose form and spectral properties are derived, and that forecasting the single reconstructed component reduces to a moving average. SSA in unit roots thus provides an `optimizing' non-parametric framework for smoothing series with a unit root. This line of work is also extended to the case of two series, both of which have a unit root but are cointegrated. The application of SSA in this bivariate framework produces a smoothed series of the common root component. Gap-filling The gap-filling versions of SSA can be used to analyze data sets that are unevenly sampled or contain missing data (Schoellhamer, 2001; Golyandina and Osipov, 2007). Schoellhamer (2001) shows that the straightforward idea to formally calculate approximate inner products omitting unknown terms is workable for long stationary time series. Golyandina and Osipov (2007) uses the idea of filling in missing entries in vectors taken from the given subspace. The recurrent and vector SSA forecasting can be considered as particular cases of filling in algorithms described in the paper. Detection of structural changes SSA can be effectively used as a non-parametric method of time series monitoring and change detection. To do that, SSA performs the subspace tracking in the following way. SSA is applied sequentially to the initial parts of the series, constructs the corresponding signal subspaces and checks the distances between these subspaces and the lagged vectors formed from the few most recent observations. If these distances become too large, a structural change is suspected to have occurred in the series (Golyandina et al., 2001, Ch.3; Moskvina and Zhigljavsky, 2003). In this way, SSA could be used for change detection not only in trends but also in the variability of the series, in the mechanism that determines dependence between different series and even in the noise structure. The method have proved to be useful in different engineering problems (e.g. Mohammad and Nishida (2011) in robotics), and has been extended to the multivariate case with corresponding analysis of detection delay and false positive rate. == Relation between SSA and other methods ==

; Autoregression: Typical model for SSA is x_n=s_n+e_n, where s_n=\sum_{k=1}^r a_k s_{n-k} (signal satisfying an LRR) and e_n is noise. The model of AR is x_n=\sum_{k=1}^r a_k x_{n-k}+ e_n. Despite these two models look similar they are very different. SSA considers AR as a noise component only. AR(1), which is red noise, is typical model of noise for Monte-Carlo SSA (Allen and Smith, 1996). ; Spectral Fourier Analysis: In contrast with Fourier analysis with fixed basis of sine and cosine functions, SSA uses an adaptive basis generated by the time series itself. As a result, the underlying model in SSA is more general and SSA can extract amplitude-modulated sine wave components with frequencies different from k/N. SSA-related methods like ESPRIT can estimate frequencies with higher resolution than spectral Fourier analysis. ; Linear Recurrence Relations: Let the signal be modeled by a series, which satisfies a linear recurrence relation s_{n}=\sum_{k=1}^r a_k s_{n-k}; that is, a series that can be represented as sums of products of exponential, polynomial and sine wave functions. This includes the sum of dumped sinusoids model whose complex-valued form is s_n=\sum_k C_k \rho_k^n e^{i 2\pi \omega_k n}. SSA-related methods allow estimation of frequencies \omega_k and exponential factors \rho_k (Golyandina and Zhigljavsky, 2013, Sect 3.8). Coefficients C_k can be estimated by the least squares method. Extension of the model, where C_k are replaced by polynomials of n, can be also considered within the SSA-related methods (Badeau et al., 2008). ; Signal Subspace methods: SSA can be considered as a subspace-based method, since it allows estimation of the signal subspace of dimension r by \mathop{\mathrm{span}}(U_1,\ldots,U_r). ; State Space Models: The main model behind SSA is x_n=s_n+e_n, where s_n=\sum_{k=1}^r a_k s_{n-k} and e_n is noise. Formally, this model belongs to the general class of state space models. The specifics of SSA is in the facts that parameter estimation is a problem of secondary importance in SSA and the data analysis procedures in SSA are nonlinear as they are based on the SVD of either trajectory or lag-covariance matrix. ; Regression: SSA is able to extract polynomial and exponential trends. However, unlike regression, SSA does not assume any parametric model which may give significant advantage when an exploratory data analysis is performed with no obvious model in hand (Golyandina et al., 2001, Ch.1). ; Linear Filters: The reconstruction of the series by SSA can be considered as adaptive linear filtration. If the window length L is small, then each eigenvector U_i=(u_1, \ldots, u_L)^\mathrm{T} generates a linear filter of width 2L-1 for reconstruction of the middle of the series \widetilde{x}_s, L\le s\le K. The filtration is non-causal. However, the so-called Last-point SSA can be used as a causal filter (Golyandina and Zhigljavsky 2013, Sect. 3.9). ; Density Estimation: Since SSA can be used as a method of data smoothing it can be used as a method of non-parametric density estimation (Golyandina et al., 2012). ==See also==