3. TOOLKIT DEMONSTRATION

Singular-Spectrum Analysis

Having selected the data vector to be analyzed (here `soi', the SOI time-series) and the units of sampling interval, the principal SSA options to be specified are the Window Length, the type of Significance Test, and Covariance Estimator. The Default button is provided as a guide to select input parameters based on the length of the data time series.
The number of SSA Components specifies how many leading components in terms of captured variance will be retained for future analyses. After computation more information about these retained componentscan be obtained by using Advanced options (see below). Results of SSA are stored in matrix with a name specified in Spectrum field. In addition, T-EOFs and T-PCs (see below) are stored in matrices with names obtained by prefixing "eof_" and "pc_" to a Spectrum name, and can be accessed in Data I/O tool. If results from several SSA runs have been stored in different matrices, the parameters used in a particular SSA run will be restored in GUI by simply selecting correspondent matrix from a Spectrum pop-up list.

Window Length

We will set the Window Length to 60, which is a good choice for out time series (690 data points with a one month sampling rate) and the oscillatory periods (2 and 4 years) under investigation.

Covariance Estimation

As indicated in the theoretical section, Burg estimation is an iterative process based on fitting an AR model with a number of AR components equal to the SSA window length. The Burg approach in principal should involve less "power leakage" due to the finite length of the time series and should therefore improve resolution (Penland et al., 1991). However, Vautard et al. (1992) found that the Burg estimate can induce significant biases when nonstationarities and very low-frequency variations are present in the series. Thus in some cases it will be worthwhile to try the Vautard et al. method. Also, for long series (N on the order of 5,000), the Vautard -Ghil estimate is less computationally burdensome and thus is completed more quickly.

In practice, we have found that the Vautard-Ghil method is more prone than is the Burg method to numerical instabilities (yielding negative eigenvalues) when pure oscillations are analyzed. Both the Burg and Vautard-Ghil methods impose a Toeplitz structure upon the autocovariance matrix whereas the Broomhead & King method does not. The Broomhead & King method is thus somewhat less prone to problems with nonstationary time series (Allen et al., 1992b), although the Vautard et al. method seems untroubled by all but the most extreme nonstationarities. The Toeplitz methods also impose symmetries on the EOF shapes whereas the Broomhead & King method does not. However the Toeplitz methods do appear to yield more stable results under the influence of minor perturbations of the time series (Allen and Smith,1996). 

NB: The Burg covariance matrix estimation option is not supported by Monte-Carlo signficance test which defaults to Vautard-Ghil if Burg is selected on the main SSA panel.

We used Burg covariance method for analyzing SOI time series.

Advanced Options

By selectinfg rows of the table, user can plot correspondent T-EOFs or T-PCs, and perform reconstruction (discussed later). In addition, time series prediction can be done with settings defined in Prediction Options. See SSA Prediction chapter for more information.

Significance tests

• Heuristic
• Chi-Squared
• Monte Carlo
  
  

Heuristic

dl_i = (2kL/N)^1/2 l_i

where L is a typical decorrelation time for the time series, k is a user supplied Decorrelation weight between 1 and 2 (1.5 by default), and l_i is the i-th eigenvalue in the spectrum. The Toolkit estimates L to be the inverse of the logarithm of the lag-one autocorrelation of the time series. Weights close to 1 are closer--but not generally equal to--the liberal error bars of Vautard and Ghil (1989); weights approaching 2 are closer to the conservative error bars of Ghil and Mo (1991a).

After we click the Compute and Plot buttons a graphics window will open to display the eigenvalue spectrum of the specified SSA. The units of abscissa are SSA component number (eigenvalue rank), and with the variance contributed by each SSA component on the ordinate.

The Toolkit provides three pairing criteria to identify oscillatory pairs and trend components in clusters of eigenvalues whose ad-hoc error bars overlap; however, it is not necessary for the 'Heuristic' test to be selected on the main SSA panel. The results of these tests are written to the log file, that can be accessed in Log tab. These tests can be activated concurrently using checkboxes in Advanced panel

• `Same Frequency' --Following Vautard et al. (1992, section 4.2), the T-EOFs associated with potential pairs or clusters are subjected to a simple Fourier transform to identify the their dominant frequency. A pair (or cluster) is identified when the associated T-EOFs have the same dominant frequency, within a fraction of the SSA bandwidth. 
  
• `Strong FFT' --Also following Vautard et al. (1992), the same Fourier transform is used to determine how much of a given signal the potential oscillatory pair accounts for at their dominant frequency. This variance fraction must exceed 95% for a pair (or cluster) to be identified.
  
• `Do trend Test' --Two tests that help to identify trending and low-frequency SSA components:
  • Kendall's tau nonparametric trend tests on the T-PCs (and RCs if present?)
  • Counts of the numbers of zero crossings within the T-EOFs.

In the SOI example, the log file shows that pair 1-2 of SSA components passes both `Same Frequency' and `Strong FFT' tests for being oscillation, while the 5th component is a trend (see here for more on SSA detrending):

We will return to investigate the oscillations after considering some more of the controls on the main SSA panel.

Monte Carlo

Two fixed T-EOF bases are used using Basis popup menu in Advanced Options:

• Data: The data, together with many AR(1) noise realizations are projected onto the EOFs of the data covariance matrix
  
  
• AR(1): The data, together with many AR(1) noise realizations are projected onto the EOFs of the expected covariance matrix of pure noise, plus any selected "signal" EOFs of the data covariance matrix which are specified as "included" (see below) in the null-hypothesis, with the rest as for pure noise, after orthogonalizing. See Allen and Smith (1996) section 4.4. This is called the "null-hypothesis basis", and it is either pure noise, or a "hybrid basis" of noise plus certain data EOFs. The noise covariance matrix is computed analytically from the estimated lag-1 autocorrelation coefficient of the test time series. 
  
  
If certain data EOFs are "included", the resulting "hybrid basis" is the most reliable if you have already identified a strong signal in your data (e.g. an annual cycle), which would otherwise overwhelm the noise parameter estimation.

`Confidence': in percentiles, e.g. 95 will mean 95% confidence level with respect to red noise null-hypothesis.

`No. Surrogates:' Number of Monte Carlo noise realizations.

Chi-squared

Using SSA for identifiying oscillatory components

A typical testing procedure for Monte Carlo or Chi-Squared tests would be as follows (i) begin by testing against pure noise and identify frequencies at which the data displays anomalously high power against this null-hypothesis; (ii) find the data EOFs which correspond to these frequencies (being data-adaptive, these will often provide a better filter for the signal than the corresponding EOFs of the noise); finally (iii) re-test including these EOFs in the null-hypothesis to check if there are other features in the spectrum which were concealed by the dominant signal . It is, of course, consistent to iterate this procedure, but beware of the problem of test multiplicity discussed in Allen and Smith (1996) section 4.2: iteration makes the test increasingly liberal, so results which are still on the margins of acceptable significance after a couple of iterations including successively more EOFs into the "signal" should be ignored. 

Accordingly, to test against a pure red-noise null-hypothesis, we choose Chi-Squared as a significance test (leaving the Include EOFs field blank). After clicking the Compute and Plot buttons, and adjsuting the X-axis limits in Graph Controls window, the following graphics window will open:

For the Data Eofs basis we observe that eigenvalues corresponding to EOFS 1-2 and 3-4 appear almost superimposed on each other at frequencies equal to ~ 0.02 and 0.034 cycle/month, and lie outside the null-hypothesis error bars. Thus they are relatively unlikely ( at the 95% level set in Confidence levels of Advanced Options) to be merely due to selected null-hypothesis process, and represent two significant oscillations. The results of the Chi-squared test should always be checked using the Monte Carlo approach, which is also essential with more complex noise models.  We leave it as an exercise to the reader to check that similar results are obtained for our SOI time series with the Monte Carlo SSA test.

One can also choose 'basis' EOFs to be computed from the expected covariance matrix of the 'Null-hypothesis' by selecting AR(1) from the basis popup menu in Advanced SSA options. Projection of data does not give SSA eigenvalues here, but the interpretation of the graph is the same. Note how the EOFs of the pure noise NH are almost regularly spaced. According to Allen and Smith ('96) the NH basis avoids artificial variance-compression inherited in SSA and therefore it has a lower probability of false-positive results, i.e. identifying the noise components as significant. After doing calculations again the resulting spectrum is shown below:

The NH basis also confirms the significance of the identified significant pairs, and we conclude that they correspond to the QB and QQ components of El Niño.

Reconstruction and Prediction.

The SSA components table on the Advanced panel allows the reconstruction and prediction of the original time series from selected SSA components, as well as plotting associated T-EOFs and T-PCs by using Plot EOFs and Plot PCs buttons. The name of vector with reconstruction/prediction is specified by the user in Result field. Individual RC components are stored in matrix with name obtained by prefixing "mat_" to Result name, and can be accessed in Data I/O tool. The Result field is shared with the SSA gap-filling feature.

If results from several reconstructions have been stored in different vectors, the parameters for a particular reconstruction will be restored in GUI by simply selecting correspondent vector from Result pop-up list. For prediction user need to specify lead time in Lead field, and order of AR model to advance in time selected SSA components in Prediction Options panel; cross-validation is available as well for finding optimal parameters for prediction, see SSA prediction for more details. Checking "fliter out" box will filter out the selected components from the orginal data (prediction feature is disabled then); this can be quite useful for detrending.

Selecting the 1-4 table rows of identified 'significant' components , and after hitting Compute in Reconstruction/Prediction box, followed by Plot, next figure will be obtained (after adjusting its parameters in Graph Controls window):