Basic Singular Spectrum Analysis and Forecasting with R

Nina Golyandina; Anton Korobeynikov

Source-linked AI summary

Basic Singular Spectrum Analysis and Forecasting with R

Nina Golyandina, Anton Korobeynikov

arXiv:1206.6910v2Jun 28, 2012stat.MEstat.AP

TL;DR

The paper addresses how SSA can be practically implemented for time-series analysis, forecasting, and parameter estimation. It describes the Rssa package and demonstrates its use through code examples and case studies, including decomposition, frequency estimation, and forecasting. The examples illustrate the methodology while also exposing parameter-selection and changing-seasonality constraints.

Problem
The paper addresses the practical implementation and use of SSA for time-series analysis, forecasting, and parameter estimation.
Method
The paper describes Rssa’s SSA implementation and demonstrates its workflows through typical R code, visualizations, and case studies.
Results
The examples demonstrate decomposition into trend, periodic components, and noise, parameter and frequency estimation, forecasting, and accuracy checking.
Takeaways & Limitations
Rssa provides a theoretically supported way to apply SSA methodology and inspect parameter choices and results through visual tools and runnable code.
Takeaways & Limitations
The package was recently created and remained under development, and a decomposition example leaves two additional sine-wave components hidden in the residuals.

Abstract

from arXiv · show

Singular Spectrum Analysis (SSA) as a tool for analysis and forecasting of time series is considered. The main features of the Rssa package, which implements the SSA algorithms and methodology in R, are described and examples of its use are presented. Analysis, forecasting and parameter estimation are demonstrated by means of case study with an accompanying code in R.

1. Introduction

The paper introduces SSA as a broad methodology for time-series analysis and forecasting, then shows how Rssa implements its analysis, forecasting, and parameter-estimation workflows through examples and case studies.

SSA and its scope: SSA comprises interlinked methods spanning decomposition, filtration, parameter estimation, and forecasting across diverse applications.Its scope includes exploratory analysis and signal-processing problems beyond conventional classical time-series analysis.
Rssa package: Rssa provides an extensively developed implementation of major SSA procedures with visual tools for parameter choice and result examination.The package is closely tied to established SSA theory and methodology, supporting its use in analysis.
Rssa package: The paper presents Rssa as the only SSA implementation then available from CRAN and describes it as almost certainly the fastest implementation.The authors qualify this assessment by noting that the package was recently created and remained under development.
Paper purpose: The paper’s aim is to demonstrate SSA analysis, forecasting, and parameter estimation using the Rssa package.The authors direct readers to specialized SSA references for fuller methodological and theoretical treatment.
Examples and case studies: Examples and case studies show decomposition into trend, periodic components, and noise, alongside parameter selection, frequency estimation, forecasting, and accuracy checking.The accompanying code also demonstrates Rssa calls and plotting commands needed to inspect parameter choices and justify results.

2. SSA algorithms and methodology

SSA decomposes a time series through trajectory-matrix embedding, eigentriple decomposition, grouping, and diagonal averaging, with separability guiding parameter choices and forecasting extending selected components.

Core SSA algorithm: SSA decomposes an observed series into interpretable components without a priori structural information, using complementary decomposition and reconstruction stages.The basic setup uses a window length L and K = N − L + 1 lagged vectors.
Embedding: Embedding forms lagged vectors and a Hankel trajectory matrix whose rows and columns are subseries of the original time series.Equal elements on anti-diagonals give the trajectory matrix its Hankel structure.
Decomposition: Decomposition represents the trajectory matrix through ordered eigenvectors and eigentriples, with Case A corresponding to the singular value decomposition of X.Eigentriples contain √λ_i, P_i, and Q_i, while eigenvalues are ordered decreasingly.
Grouping and reconstruction: Eigentriple grouping partitions nonzero components into disjoint subsets, and diagonal averaging converts each grouped matrix back into reconstructed series.The reconstructed series from grouped matrices sum to the initial series.
Separability and parameters: Separability determines whether components such as signal, noise, sine waves, trend, and seasonality can be approximately extracted and guides window-length selection.The w-correlation matrix identifies component groups: small weighted correlations indicate better separation, whereas correlated components should not be placed in different groups.
Separability and parameters: A w-correlation matrix contains weighted correlations between reconstructed components, helping identify correlated groups for subsequent eigentriple grouping.The weights reflect how often time-series terms enter the trajectory matrix.
Sequential SSA: Sequential SSA extracts trend with a small window and periodic components from residuals using L ∼ N/2 for complex series.For exploratory, model-free use, decomposition justification is based on separability theory and interpretability; statistical testing requires a concrete model.
Forecasting: Recurrent forecasting extends a series by applying a linear recurrence relation, while vector forecasting constructs lagged vectors, diagonal-averages them, and is typically more stable but costlier.The recurrent forecast produces M future terms directly from recurrence coefficients; when a separated component obeys an LRR, the two forecasts coincide and provide exact continuation.

2.4. Linear recurrence relations, time series of finite rank and roots

SSA forecasting models signals that are approximately governed by linear recurrence relations, whose characteristic roots encode signal structure and parameters. Forecasting and confidence assessment depend on appropriate model assumptions, root selection, and parameter choices.

Finite-rank series: Signals admitting exact SSA continuation are governed by linear recurrence relations, and their trajectory matrices are rank-deficient.These series provide a natural signal model for SSA forecasting.
Linear recurrence relations: The minimal LRR is unique, while the recurrent SSA forecast uses the best LRR among non-minimal relations of order L−1.The selected relation is discussed as preferable among the available non-minimal alternatives.
Characteristic roots: For real-valued series governed by an LRR, the resulting class consists of sums of polynomial, exponential, and sinusoidal terms.Characteristic-root multiplicities determine the polynomial factors, while complex roots produce oscillatory components.
Parameter estimation: Characteristic roots support parameter estimation: conjugate-root arguments yield sinusoidal frequencies, while root modulus gives α = ln ρ.For an exponentially modulated sinusoid, conjugate roots are written as ρe±i2πω.
Frequency estimation: ESPRIT estimates signal roots from eigenvalues of a pseudo-inverse matrix formed by deleting the first and last rows of selected eigenvectors.The estimated frequencies are the arguments of the resulting eigenvalues.
Forecast uncertainty: Bootstrap intervals for reconstructed signals and forecasts assume a rank-r signal with white-noise residuals and repeatedly resample the empirical residual distribution.The intervals use sample quantiles, while the sample mean is reported as the average bootstrap forecast.

3. Rssa package

Rssa implements the main SSA workflow in R, combining decomposition, reconstruction, forecasting, visualization, caching, and multiple SVD strategies. Its computational design automatically selects or reuses procedures while supporting large-window analyses, subject to documented data and memory constraints.

Core interface: The ssa function performs embedding and optionally decomposition, accepts standard numeric or time-series inputs, and returns an SSA object for downstream functions.For one-dimensional series, the paper focuses on kind="1d-ssa" and also discusses kind="toeplitz-ssa".
Core interface: reconstruct uses eigentriple groupings to return reconstructed series while preserving input time-series attributes unless drop is requested.The grouping is supplied as a list of index vectors for elementary series.
Object management: Rssa automatically calculates missing eigentriples when needed and reuses previous results through internal caching to save computation time.Previously calculated elementary series are stored inside the SSA object.
Object management: SSA objects use reference-like memory management, so deep copies require clone rather than the standard assignment operator, and cached data can be released with cleanup.This design minimizes memory copying but adds semantic constraints to object handling.
SVD implementations: Rssa provides truncated Lanczos-based SVD methods, propack and nutrlan, alongside full eigen and svd decompositions.The full methods ignore neig, whereas truncated methods target selected leading eigentriples.
Computational strategy: Lanczos SVD methods are most effective for large window lengths and non-small series, enabling larger windows that can improve separability.For a series of length 87000 and window length 43500, propack took 13 seconds versus 16 seconds for nutrlan.

4. Basic SSA with R

The Rssa package implements SSA analysis through decomposition, reconstruction, diagnostics, and visualization, illustrated first with “co2” and then with a sequential SSA analysis of “MotorVehicle”.

Basic SSA with R: Rssa decomposes time series by grouping reconstructed elementary components, with diagnostic plots supporting parameter and grouping choices.Eigenvalue, eigenvector, paired-eigenvector, and w-correlation plots are used to inspect components and separability.
Basic SSA with R: The “co2” example reconstructs trend, seasonality, and noise from the observed series using grouped eigentriples.Residual and original-series options allow the decomposition to be plotted as a sum of trend, seasonality, and noise.
Basic SSA with R: The “co2” result depends on choosing the window length and grouping appropriately, because additional sine-wave components remain hidden in the residuals.The paper presents the example primarily for code demonstration and refers readers elsewhere for detailed parameter-selection principles.
Case study: For “MotorVehicle”, Sequential SSA first extracts a complex trend with L = 12 before analyzing seasonality in the residual.The first-stage window length is chosen for smoothing, with divisibility by the periodic component’s period used as the rationale.
Case study: The first “MotorVehicle” eigentriple represents the trend, while later eigentriples contain high-frequency components unrelated to the trend.The leading eigenvector is nearly constant and corresponds to smoothing by the Bartlett filter.
Case study: The residual periodogram identifies seasonal periods 12, 6, 4, 3, and 2.4, which are extracted using a second SSA stage with L = 264.Eigenvalue steps, eigenvector scatterplots, and w-correlation diagnostics support pairing and separation of the seasonal components.

Fragment 4.8: “MotorVehicle”, 2nd stage: reconstruction and plotting of the results

Sequential SSA reconstructs the extracted seasonality and combines both stages to decompose MotorVehicle into trend, seasonality, and residuals. The resulting seasonally adjusted series is also produced.

The second-stage reconstruction extracts seasonality from eigentriples ET1–10 and plots it alongside the MotorVehicle series.The extracted component reflects complex periodic behavior caused by slowly changing sine-wave phases.
The resultant decomposition combines the first- and second-stage Sequential SSA results.Figure 10 presents the decomposition produced from both stages.
The reconstruction procedure also produces a seasonally adjusted MotorVehicle series.This auxiliary output is shown in Figure 11.
The residuals from the decomposition are heterogeneous rather than uniform.The residual behavior motivates subsequent estimation of a heterogeneous-noise envelope.

Fragment 4.9: “MotorVehicle”: finding noise envelope

The noise-envelope procedure estimates time-varying residual variability by extracting the trend of squared residuals with SSA, while Toeplitz SSA is unsuitable here because seasonality changes over time.

The residual envelope is computed by applying SSA to squared residuals and taking the square root of the reconstructed leading component.The resulting positive and negative envelope curves are plotted around the residual series.
Toeplitz SSA produces a worse decomposition for the MotorVehicle seasonality because the seasonal behavior changes over time.The paper notes that Toeplitz SSA is appropriate for stationary series but not for trend extraction in this example.
The ordering of eigentriples by matrix eigenvalues can differ from their contribution ordering, making the eigenvalue plot non-monotonic.The values of s$lambda equal contribution values, while ordering is performed by eigenvalues.

5. SSA forecasting with R

Rssa provides recurrent and vector SSA forecasts, tools for inspecting forecasting recurrence relations, and parameter-estimation workflows. In the co2 example, estimated signal roots support sine-wave and exponential-trend forecasts.

Rssa implements recurrent and vector forecasting after SSA decomposition.The package also provides confidence-interval forecasting and an all-in-one forecast wrapper compatible with the forecast package.
Forecasting routines attempt to preserve initial-series attributes, but time-index imputation requires separate handling for supported time-series classes.Standard classes such as ts and zooreg are handled explicitly.
The co2 forecast is represented by half-year and annual sine waves with nearly constant amplitude plus a trend approximated by two exponentials.This interpretation follows analysis of the forecasting linear-recurrence roots and ESPRIT estimates.
The forecasting workflow evaluates different window lengths and selected eigentriple groups using recurrent or vector forecasts.Forecast accuracy can be tested with sliding forecasts before choosing method parameters.

Fragment 5.4: Function for sliding forecasts

The sliding-forecast function repeatedly trains SSA on moving subseries, forecasts held-out points, computes mean squared errors, and averages them across windows.

The function selects recurrent or vector forecasting through the type argument and iterates over supplied component groups.Each sliding window produces a forecast for every requested group.
Each iteration fits SSA to a training subseries and forecasts the following forecast.len observations.The forecast is compared with the corresponding held-out segment.
For each group and window, the function records mean squared error and returns its average across sliding windows.This provides a parameter-comparison criterion for the selected forecasting setup.
With a 468-point series and 360-point sliding subseries, fl = 1 yields 108 one-step forecasts, whereas fl = 108 yields one 108-step forecast.These settings represent short-term and long-term forecasting evaluations, respectively.

Fragment 5.5: “co2”: dependence of forecast accuracy on choice of parameters

The examples use Rssa code to evaluate forecasting choices for time-series components, including window, grouping, and the number of removed starting points. For MotorVehicle, seasonality forecasts are stable across recurrent and vector methods, while the whole series with ET1–8 or ET1–10 is preferred.

co2 parameter choice: The co2 example varies the SSA window through Lseq and evaluates one-step forecasts across sliding windows.The code applies forecast.check for each candidate L and compares the resulting accuracy measures graphically.
MotorVehicle seasonality: MotorVehicle seasonality is forecast separately from the trend because the trend has a complex structure.The workflow reconstructs the trend, computes residual seasonality, and forecasts it with groups 1:10.
MotorVehicle seasonality: Recurrent and vector forecasts are stable enough that their difference is very small.Bootstrap confidence intervals can be calculated for the forecasted component to assess forecast error.
MotorVehicle seasonality: Forecasting based on the whole time series performs better with ET1–8 or ET1–10.The comparison is made after varying the number of removed starting points, using forecast.check on the residual trend.
MotorVehicle trend: A changing trend structure makes using the whole trend unreasonable for forecasting.The analysis therefore seeks the last point of structural change before selecting the forecasting interval.

Fragment 5.9: “MotorVehicle” trend: forecasting of last 3 years behavior

The short-horizon MotorVehicle trend forecast uses the final 36 points, corresponding to the last three years. Its structural-change boundary remains uncertain and requires additional study.

Forecasting the last 3 years: The resulting short-period trend forecast differs from the forecast obtained from a longer time period.The paper presents the longer-period comparison in Fig. 16.
Forecasting the last 3 years: The proposed structural-change point needs additional study using an SSA heterogeneity matrix or other change-point methods.The 2009 date is described as probable rather than established.

Fragment 5.10: “MotorVehicle” trend: forecasting of last 22 years behavior

The long-horizon MotorVehicle trend forecast uses a 22-year segment, producing a substantially different forecast from the last-three-years analysis. Choosing between them requires additional macroeconomic analysis.

Forecasting the last 22 years: The 22-year analysis constructs trend.end2 from trend[270:541] and forecasts it with SSA components 1:4.This is the longer-period forecast contrasted with the last-three-years result in Fig. 16.
Forecasting the last 22 years: Selecting the proper forecast requires additional macroeconomic analysis.The SSA comparison alone does not determine which forecast should be used.