g. Multi-Channel Singular-Spectrum Analysis

As an example, we apply M-SSA to the near-global data set of monthly sea-surface temperatures from the Global Sea-Ice and Sea-Surface Temperature (GISST) data set Rayner et al.-1995 for 1950-94, from 30$^{\circ}$S to 60$^{\circ}$N, on a 4$^{\circ}$-latitude by 5$^{\circ}$-longitude grid.

This translates into a dataset with 1295 channels, 540 months long. To read the GISST data we use the Read Matrix function from the File/Data menu on the main panel, which places the data into a matrix with a default name "mat" of 540x1295 size.

Selecting the `MSSA' option from the Analysis Tools menu on the main panel launches the following window (showing its state after pressing Get Default Values button, and specifying the parameters as described below):

Having specified the data to be analyzed (here `mat', the GISST time-series) and the sampling interval, the principal MSSA options to be specified are the Window Length, the type of Significance Test, Covariance, preprocessing PCA switch, and number of PCA channels(if necessary). The number of MSSA Components specifies how many components will be retained for future analyses. The Get Default Values button is provided as a guide.

Window length

For comparison with previous ENSO studies and to demonstrate the usage of N'-windows vs. M-windows [see Eq. (11) and Fig.1, for the definition of ``complementary windows''], we set M=480, which yields an N'-window of 5 years, and M=270, which yields the M- and the N'-windows equal both to N/2.

Note!: The value in the Window Length is taken to equal to N' if `Reduced' option is chosen for Covariance, and M for other Covariance choices.

Preprocessing with PCA

We set the pre-PCA switch to Yes, and theNo. of PCA channels to 10 to prefilter data with standard PCA Preisendorfer-1988 to retain the 10 leading spatial PCs that describe 55.2% of the variance, as we can see from the Log file:

This favors the association of larger decorrelation times with larger spatial scales, as expected for climatic Fraedrich and Boettger-1978 and other geophysical fields, and the channels are uncorrelated at zero lag. After PCA the input data for MSSA has L=10 channels and length of N=540 months, and it is stored in the list of matrices as SPC; correspondent spatial EOFs are stored in a matrix SEOF

Covariance Estimation

The method for estimating the lag-covariance matrix that is decomposed (diagonalized) by M-SSA is chosen by selecting either `Reduced' , `Vautard-Ghil' , or ` Broomhead &King ' from the `Covariance' menu on the main M-SSA control panel (Fig. 11).
We choose Reduced for the Covariance estimator, since in all cases ML > N' so that it is more efficient to diagonalize the reduced ( $N'\times N'$) covariance matrix with elements given by Eq. (11), rather than the ( $ML \times ML$) matrix whose elements are given by Eq. (6) (` Broomhead &King), or Eq. (1) `Vautard-Ghil'.

We set the Window length to 60 and 270 (N' for such option, see above).

Significance test

There are three choices for Significance test

described below. If None is selected, the eigenspectrum is displayed in order of eigenvalue rank. The leading oscillatory pair over the entire domain has a quasi-quadrennial period for all values of M. The quasi-quadrennial pair (modes 2 and 3 in this case ), accounts for less then 30% vof ariance for N'=60, as we can see by opening the Log File shown in the next figure:

The smaller one of M and N' determines the approximate spectral resolution 1/N' or 1/M. Choosing M=N'=270 yields the maximum spectral resolution.

The dominant frequencies of MSSA modes are computed only in Monte-Carlo or apporximate, but much faster Chi-squared MSSA test. Choosing Chi-squared in the Significance tests menu, and M=N'=270 we obtain the following plot:

Here the red-noise surrogate projections are plotted against the dominant frequencies associated with each MSSA mode, and we have zoomed in on the 0-0.05 cy/month frequency interval of interest. Using M=N'=270 captures less variance - 12.1% - for the significant oscillatory quasi-quadriennial pair formed in this case by modes 3 and 4, as seen in the Log File:

Monte-Carlo Test Options

The quasi-biennial pair (modes 14 and 15) does not pass the significance test when data eigenmodes are used as a basis for projections. However since the Monte-Carlo test is biased if we project onto the data eigenmodes, we project onto the eigenmodes provided by the covariance matrix of the AR(1) noise. So, we select AR(1) instead of Data option in Test Options pull-down menu:

to obtain the following result:

Here, the projections are plotted against the dominant frequencies associated with each noise eigenvector, and we have zoomed in on the 0-0.05 cy/month frequency interval of interest using the xmgr plotting tool. Since the latter are near-sinusoidal in this case, the resulting spectrum is closely related to a traditional Fourier power spectrum. Both the quasi-quadrennial(~0.023 cycle/month) and the quasi-biennial(~0.038 cycle/month) modes pass the test at the 95% level (as specified in Test Options). They are well separated in frequency by about 1/(20 months), which far exceeds the spectral resolution of 1/M = 1/N' = $1/(270\;\;{\rm months}) \approx 1/(22\;\; {\rm years})$. The two modes are thus significantly distinct from each other spectrally, in agreement with the univariate SOI results using MEM and MTM, respectively.

We leave to a user to verify the above results with the Monte-Carlo test with 100 surrogates for M=N'=270, the maximum effective resolution. Computation make take a while, so be patient!

Plot Options

We can check that ST-EOFs of oscillatory MSSA pairs are indeed in phase quadrature using Plot Options pull-down menu:

The following plot is for a quasi-quadriennial pair 3 and 4 and channel 1:

We leave as an exersise to check that a phase quadrature for a quasi-biennial pair (modes 14 and 15) is a mostly prominent for channel 2.


We can reconstruct the original multi-channel data or PCA components from selected MSSA modes using the Reconstruction pull-down menu in the MSSA window.

We can reconstruct either spatial PCs if PCA-prepocessing has been used or original data: choose PCA or Grid options, respectively. Here we show reconstruction in PCA space quasi-quadriennial pair for channel 1: