Our motivation to focus on advanced spatio-temporal data analysis is to better understand the decadal climate variability in the Earth system and illuminate the capabilities of prediction tools to capture the associated signals (Meehl et al., 2014). Inter- and intra-decadal climate variability is inherent to the ocean–atmosphere system and is further coupled to other Earth system components, such as sea-ice and land surface (Meehl et al., 2009). The variability appears as complex four-dimensional (or spatio-temporal) structures in Earth system variables, such as wind, temperature and precipitation (Solomon et al., 2011).
These structures are embedded in extremely large-dimensional data sets gathered and generated in reanalysis of atmospheric and oceanic observations, and in massive simulation endeavours using Earth system models world-wide. Applicability of advanced data analysis tools is severely hampered by the very large dimensionality of the climate data.
Many common analysis methods, such as principal component analysis (PCA; Von Storch and Zwiers, 1999), involve eigen-problems, which become impossible to solve with increasing data dimension. Earlier we illustrated the use of random projections (RP) as a tool to tackle high-dimensional problems (Seitola et al., 2014). We demonstrated how PCA can be applied in three-dimensions to problems that are beyond practical computational limits without efficient dimension reduction. PCA is not an ideal tool, however, to extract and illustrate four-dimensional eigen-features in climate data. In this respect, the multichannel singular spectrum analysis (MSSA; Broomhead and King, 1986a, b) is a much more appealing method since the MSSA eigen-problem inherently contains the auto-covariance in the lagged copies of the original data vectors. The computational burden is, however, even larger than in PCA. We overcome this burden by a novel randomised version of MSSA, called RMSSA. To our knowledge, this approach has not been suggested before. We note that Oropeza and Sacchi (2011) incorporate a randomising operator into MSSA for noise attenuation in seismic data, but their algorithm is not aimed directly at large-dimensional problems. In RMSSA, RPs are used essentially to enable analysis of extremely large-dimensional data sets.
In this article, we present the RMSSA algorithm in detail and also include a test for the statistical significance of the results (Monte-Carlo MSSA; Allen and Robertson, 1996) in the algorithm. We demonstrate the use of RMSSA by decomposing the 20th century global monthly mean near-surface temperature variability into its low-frequency components. The data sources are described in Section 2.3.
MSSA was introduced into the study of dynamical systems by Broomhead and King (1986a, b). The method is equivalent to extended empirical orthogonal function (EEOF) analysis (Weare and Nasstrom, 1982), but there are differences in the choice of some important parameters and in the interpretation of the results (Plaut and Vautard, 1994).
In traditional PCA or EOF analysis (e.g. Rinne and Karhila, 1979), spatial correlations (in case of climatic data sets) are used in determining the patterns that explain most of the variability in the data set, but MSSA differs from this traditional method by also taking into account the temporal correlations. In other words, standard PCA decomposes a spatio-temporal field into spatial PC loading patterns (EOFs) and corresponding PC score time series (PCs), whereas MSSA also adds a temporal dimension to EOFs. MSSA PCs and EOFs are often called space-time PCs (ST-PCs) and space-time EOFs (ST-EOFs), and we have adopted this notation here. A more detailed description of MSSA is presented in Ghil et al. (2002) and in Appendix A.1 here.
The idea of MSSA, in brief, is to find the patterns that maximise the lagged covariance of the data set X_{N×L} within M lags. In case of a gridded climate data set, N represents the time steps and L is the number of grid points. The columns of the data matrix X are often called channels. The length of the lag window M is a user choice. For example, Elsner and Tsonis (1996) suggest that the results of MSSA do not change significantly with varying M as long as N>>M and they recommend using M=N/4. Vautard and Ghil (1989) recommend to choose M no larger than approximately N/3. Clearly, if the number of channels L is large in the beginning, choosing large M would result in a very high-dimensional data matrix with M×L columns, including lagged copies of each channel in X.
Determining the length of the lag window M is a trade-off between spectral resolution and statistical significance of the obtained components. The larger M is chosen, the more temporal information can be extracted but at the same time the variance is distributed on a larger set of components. If M is small, the statistical significance of the obtained components is enhanced. In this study, we used several values of M in order to test its effects on the results.
ST-PCs/ST-EOFs often appear in pairs ('sinusoidal’) that explain approximately the same variance and are π/2 out of phase with each other. These pairs are said to present stationary or propagating oscillatory modes of the data set (Plaut and Vautard, 1994). Modes with period less than or equal to M can be only presented by such pairs. However, existence of such a pair does not guarantee any physical oscillation, and according to Allen and Robertson (1996) such pairs can also be generated by non-oscillatory processes, such as first-order autoregressive noise.
This finding led Allen and Robertson (1996) to formulate a test for the statistical significance of MSSA components. The identified components are tested against a null-hypothesis of the data being generated by independent AR(1) processes (i.e. red noise) with the same variance and lag-1 autocorrelation as the original input time series. This procedure is called Monte-Carlo MSSA (MC-MSSA), and it is described in more detail in the original study of Allen and Robertson (1996) as well as in Appendix A.1 of this article.
ST-PCs cannot be compared to the original time series as such; instead, they can be represented in the original coordinate system by their reconstructed components, RCs (Plaut and Vautard, 1994; Ghil et al., 2002). In the reconstruction, the ST-PCs are projected back onto the eigenvectors (ST-EOFs) and each RC is a kind of filtered version of the original multivariate time series. Construction of RCs is illustrated in Fig. 1. Several ST-PCs/ST-EOFs can be used in the reconstructions, and if there is an oscillation that appears as a sinusoidal pair, both of these ST-PCs/ST-EOFs should be included in the reconstruction of that certain oscillatory mode. This is done by summing up the corresponding RCs. No information is lost in the reconstruction, and the original time series is a sum of all individual RCs.
As mentioned earlier, the computational burden of MSSA becomes soon prohibitively high if the original data set is high-dimensional and M is chosen to be large. This is typically the situation in studies of low-frequency variability in climate data sets. Traditionally, the dimensionality reduction has been obtained by calculating first a conventional PCA and retaining a set of dominant PCs for the MSSA (see chapter 2.2.3). However, in this article we apply a different approach to dimensionality reduction. That is, we use RP to reduce the dimensionality of the original data set before performing MSSA.
In Halko et al. (2011), it is stated that randomised methods provide a powerful tool for constructing approximate matrix factorisations. Compared with standard deterministic algorithms, the randomised methods are often faster and more robust. Halko et al. (2011) present also numerical evidence that these algorithms succeed for real computational problems.
In our approach, RP is applied to reduce the dimension of the original data matrix X after which the traditional MSSA calculation is performed in the lower-dimensional subspace. Finally, we reconstruct the ST-EOFs and RCs in the original space. We call this algorithm randomised multichannel singular spectrum analysis (RMSSA).
In RP, the original data set is projected onto a matrix R of Gaussian distributed (zero mean and unit variance) random numbers in order to construct a lower-dimensional representation P of the data set:
In other words, we are projecting our data set onto k random directions determined by the column vectors of R. From these projections a lower-dimensional representation of the original data set can be constructed. Due to the simplicity of RP, involving only matrix multiplication, it can be applied to a wide range of data sets, even very high-dimensional ones.
RP has already been applied to climate data in Seitola et al. (2014) and it has been shown to preserve structures of the original data very well. In that article, the theoretical background of RP is presented in more detail with additional references.
The projected lower-dimensional data set P can be processed through MSSA where instead of original L channels we have now only k channels. This implies substantial computational savings (see Appendix A.2, algorithm 1). In the literature, there are some estimates of an appropriate value for k (e.g. Frankl and Maehara, 1988; Dasgupta and Gupta, 2003). However, these theoretical lower bounds for k are the worst case estimates and usually much lower values for k still give good results, retaining most of the information of the original data set (see e.g. Bingham and Mannila, 2001; Seitola et al., 2014). In practice, the value for k is usually chosen adaptively keeping the desired size for lower-dimensional approximation in mind.
A final step of the algorithm is to calculate the reconstructed components. This requires the recovery of the MSSA eigenmodes from the reduced space back to the original space, allowing the reconstruction of the original time series. This means that the eigenvectors (ST-EOFs) should be calculated in the original space instead of the reduced one. This part of the algorithm is also presented in Appendix A.2. Furthermore, in Appendix A.4 we explain how RP preserves the lagged covariance structure of the original data set.
To our knowledge, the proposed RMSSA algorithm is unique. Some published work comes close to our approach but RMSSA has some important differences to the randomised MSSA algorithms used in seismic data processing (Oropeza and Sacchi, 2011; Chiu, 2013). The aim of Chiu (2013) was to introduce a new rank-based-reduction denoising algorithm to perform coherent and random noise filtering concurrently. Chiu (2013) named this algorithm, or rather filter, MSSARD (MSSA in the randomised domain). In MSSARD, the randomising operator randomly rearranges the order of the input data and reorganises the coherent noise into incoherent noise. The most important difference to our algorithm is in the randomising operator: In our case, we are using RP to reduce the dimensionality of the input data whereas Chiu's (2013) approach is to randomly rearrange the input data.
The technique of Oropeza and Sacchi (2011) was to embed a spatial data at a given temporal frequency into a block Hankel matrix after which a randomised singular value decomposition (SVD) was adopted to accelerate the rank reduction stage of the algorithm. Construction of a Hankel matrix corresponds to the construction of an augmented data matrix A in our algorithm (see Appendix A.1). Our algorithm is different in the sense that we apply RP on the original input data before construction of the augmented (or Hankel) matrix. This notably reduces the computational burden of MSSA because we are processing a much smaller data set already in the augmentation phase of the algorithm (see algorithm 1 in Appendix A.2).
In addition to these main differences, the above-mentioned seismological applications involve handling a data set where each time/frequency slice of spatial (x-y) data is processed separately through the algorithm. In our case, we are processing the whole time–longitude–latitude data set at once through the RMSSA algorithm.
In many studies, where MSSA is used as an analysis method (e.g. Plaut and Vautard, 1994; Moron et al., 2012), the dimension of the original data matrix has been reduced by calculating a conventional PCA of the original data matrix and then limiting MSSA into the dominant PCs. One has to bear in mind that the problem dimension may be prohibitive to contemplate solving even PCA, let alone MSSA. Nevertheless, the number of retained PCs is a somewhat arbitrary choice, but can be estimated by inspecting the eigenvalue spectrum and choosing the PCs that account for the majority of the variance and are separated from the rest of the spectrum. In geophysical datasets, however, the eigenvalue spectrum often decreases monotonically and it is difficult to distinguish the appropriate cut-off point. The aim of the study does also affect the choice of the PCs. For example, if the focus is on large-scale patterns, it might be more convenient to choose the low-frequency PCs for further analysis. Performing the calculations with different number of PCs and comparing the results can also help in finding the appropriate number of PCs. Importantly, RMSSA (Appendix A.2, algorithm 1) does not suffer from this problem because the lower-dimensional data set has essentially the same structure as the original high-dimensional data set.
PCA-based dimensionality reduction is, however, a preferred method if the oscillatory modes identified with MSSA are tested against a red noise null-hypothesis through Monte-Carlo simulation. According to Allen and Robertson (1996) the test is only useful if the channels in the data matrix are orthogonal or at least very nearly so. The PCs fulfil the orthogonality condition exactly. The randomised method can still accelerate – and in the case of a very-high-dimensional data set even enable – the calculation of the PCs (see Appendix A.2, algorithm 2).
This also raises the question as to whether the projected data set [i.e. matrix P in eq. (1)] could be used directly in MC-MSSA. Like the PCs, RP is also an orthogonal projection and the columns of P are also nearly orthogonal. However, this question is beyond the scope of this study and will not be discussed here any further.
As an illustration of applying the RMSSA algorithm, we analysed the monthly mean near-surface air temperature field from the 20th Century Reanalysis V2 data, hereafter 20CR, provided by the NOAA/OAR/ESRL PSD (Compo et al., 2011). In addition, we repeated the analysis for the historical 20th century simulations by Hadley Global Environment Model 2 – Earth System HadGEM2-ES (Collins et al., 2011), hereafter HadGEM2, and MPI Earth System Model (ESM) running on a medium resolution grid MPI-ESM-MR (Stevens et al., 2013), hereafter MPI-ESM. We extended the historical simulations (1901–2005) until 2012 using the rcp45 simulations. The historical and rcp45 simulations were extracted from the CMIP5 data archive and they follow the CMIP5 experimental protocol (Taylor et al., 2012). In the 20th century simulations, the historical record of climate forcing factors such as greenhouse gases, aerosols and natural forcings such as solar and volcanic changes is used. Rcp45 simulations follow the RCP4.5 greenhouse gas scenario. We used a single ensemble member of each model: r2i1p1 in case of HadGEM2 and r1i1p1 in case of MPI-ESM.
The 20CR data set is produced using an ensemble of perturbed reanalyses, and the final data set corresponds to the ensemble mean. Only surface pressure observations are assimilated, and the observed monthly sea-surface temperature and sea-ice distributions are used as boundary conditions to generate full three-dimensional estimates of the state of the troposphere (Compo et al., 2011). The 20CR data set is available from 1871 to 2012 but to be consistent with HadGEM2 and MPI-ESM simulations, the time sequence analysed here is 1901–2012 (1344 time steps). 20CR has ~2.0 degree horizontal resolution and we have used Gaussian gridded (192×94) data from 3-hour forecast values. HadGEM2 and MPI-ESM have both a global grid of 144×73 points. Thus, we have original data sets X_{N×L} with N=1344, L=18048 (20CR) and L=10512 (HadGEM2 and MPI-ESM).
As an illustrative example of the high-dimensionality of the MSSA problem, let's choose a lag window of M=240 (months). In the case of the 20CR data set, this would result in an augmented matrix with M×L=4331520 columns. Clearly some kind of dimensionality reduction is needed in order to make the computations more efficient or even make them possible.
In the previous section, we have introduced the RMSSA algorithm and the data sets to be analysed. Next we will proceed to the applications of the proposed method and discuss the results.
First, the original data sets were mean centred and RMSSA (algorithm 1 in Appendix A.2) was applied with k=500. The first 1–30 ST-PCs of 20CR are shown in Fig. 2. In order to find the most powerful frequencies associated with the ST-PCs, the Multitaper spectral analysis method (Thomson, 1982; Mann and Lees, 1996) was applied. The power spectra of the ST-PCs are shown on the right in Fig. 2. The first pair of ST-PCs is clearly related to the annual cycle and this pair together explains the majority of the variance of the data set (almost 90%). The pairs of ST-PCs 3–4, 7–8 and 12–13 are the subharmonic frequencies of the annual cycle. The periods of ST-PCs 5, 6 and 11 as well as of ST-PCs 14, 17 and 18 fall outside the lag window length M and are the so-called trend components. ST-PC5 may be related to a centennial scale warming trend and ST-PC11 has a multi-decadal scale variability. ST-PCs 22 and 24 have clear spectral peaks on a 5–6 yr period and ST-PCs 29 and 30 are oscillating on a period of 3–4 yr. Those ST-PCs might be related to the El Niño-Southern Oscillation (ENSO) which is a prominent phenomenon on those time scales. ST-PCs 19–21 are related to a decadal scale variability, but the spectra of those components are quite broad on a 10–20 yr time scale.
The above analysis was also performed for the HadGEM2 and MPI-ESM data sets (figures not shown). As the annual cycle is too dominating in each data set, the analysis in the following sections will be repeated without the annual cycle. We also integrate a MC-MSSA step in the RMSSA algorithm (Appendix A.2, algorithm 2) in order to study the statistical significance of the obtained components.
Some pre-processing of the original data sets was crucial in order to assess the statistical significance of the low-frequency variability using MC-MSSA. First of all, the original data sets were standardised (i.e. the time series of each grid point was mean centred and divided by its standard deviation) in order to avoid overweighting the grid points with higher variance. Furthermore, the annual cycle of the time series of each grid point was estimated by STL (Loess based Seasonal-Trend Decomposition) and removed from the original data set. The STL method is a filtering procedure for decomposing a time series into trend, seasonal and remainder components. It includes some parameter choices controlling, for instance, how rapidly the trend and seasonal components can change. The method is described in detail in Cleveland et al. (1990) and we have followed their guidelines in choosing the related parameters. Without this procedure the annual cycle would dominate the results and starve the lower ranked MSSA components of power when tested against the red noise null-hypothesis. Linear trends were also fitted and removed from the data sets in order to avoid the dominance of the centennial scale trend.
For the sake of comparison, the annual cycle was also estimated by calculating the mean values of each calendar month and those values were subtracted from the data to get monthly anomaly time series. However, determining the base for the anomaly calculation is not that straightforward and the choice of a base period may have severe impacts on the results (Kawale et al., 2011). Furthermore, the average annual cycle is only removed and if the annual cycle varies in the time series, the anomalies still contain a residual annual cycle.
The dimensions of the original data sets were reduced by applying RP with k=500 to have a lower-dimensional approximation P_{N×k} of each data set. To be able to perform MC-MSSA, we further calculated SVD of P and retained 30 first PCs of each data set, explaining approximately 72% (20CR), 67% (HadGEM2) and 64% (MPI-ESM) of the variance. Those 30 PCs were used as input channels in the MC-MSSA-step.
The ST-PCs 1–30 of each data set and their spectra are presented in Figs. (3–5). These figures show the results after applying the steps 1–8 of algorithm 2 in Appendix A.2 (note that the annual cycle and linear trend were removed from the original data sets). In 20CR (Fig. 3), the ST-PCs 1–2 are so-called trend components explaining together almost 9% of the variability of the data set. Pairs of ST-PCs 3–4 and 5–6 in 20CR have clear peaks in frequencies corresponding to 3–4 yr and over 5 yr periods. In addition, 2–3 yr periodicities are distributed on several ST-PCs beginning from the 14th one. When the model simulations are compared to the 20CR components, the main differences are the prominent decadal scale components of HadGEM2 (ST-PCs 2–3, 9.3% of explained variance) and the 2–7 yr variability of MPI-ESM that is distributed on a large set of successive components. For more details, the readers are advised to study Figs. (3–5).
In MC-MSSA step, in total of 1000 realisations of red-noise surrogates were generated and the red-noise basis was used to estimate the 90, 95 and 99% confidence intervals for the eigenvalues generated by the noise model that consists of independent first-order autoregressive processes. Figure 6 shows the results of the Monte-Carlo significance test of 20CR, HadGEM2 and MPI-ESM with a 20 yr lag window (M=240 months). In that figure, the data eigenvalues and 2.5th and 97.5th percentiles of the distribution of the surrogate eigenvalues are plotted against the dominant frequencies of the corresponding red-noise basis vectors (noise ST-EOFs). The dominant frequencies are estimated using fast Fourier transform (FFT). It should be noted, that the estimate of the dominant frequency of the noise ST-EOFs may not be exactly the same as the dominant frequency of the data ST-EOFs which may cause some small inaccuracies in the results.
The significant signals (at 5% significance level) in Fig. 6 are those whose data eigenvalues lie above the 97.5th percentiles of the surrogate eigenvalues. According to the test these signals have more variance than would be expected to have from a noise process. According to Plaut and Vautard (1994) the use of a lag window length M typically allows the distinction of oscillations with periods in the range [M/5, M]. Therefore we only show the significance test of the periodicities that are covered by the 20 yr lag window used in this example. From Fig. 6, we can see that in 20CR data set there are some significant periodicities (at 5% level) between 1.7 and 5.5 yr. HadGEM2 has somewhat more significant periodicities compared to 20CR, especially on 10 yr time scales, but MPI-ESM has hardly any eigenvalues lying above the 97.5th percentile.
As noted earlier, the Monte-Carlo simulations were performed with varying lag window M to estimate its effect on the statistical significance of the oscillations. Spectral resolution increases with lag window length and oscillatory pairs with longer periodicity can be identified. However, at the same time the statistical significance of the identified oscillations may decline. We used the following values of M: 5 yr (M=60 months), 10 yr (M=120), 20 yr (M=240), approx. 28 yr [M=340≈N/4, following the recommendation of Elsner and Tsonis (1996)] and approx. 38 yr [M=450≈N/3, following Vautard and Ghil (1989)].
The identified periodicities and their significance levels with increasing lag window are presented in Fig. 7. The numbers in Fig. 7 show the dominant periods associated with the oscillations. These dominant periods are estimated using FFT. From Fig. 7 we can see that in 20CR the significant periodicities are consistently found at 3.6, 2.3 and 1.7 yr, depending to some extent on M. Those periods are more or less visible in HadGEM2 and to a lesser extent in MPI-ESM. Significant 5–6 yr oscillations are identified in all the data sets and especially a ~5.5 yr variability is found consistently.
2–6 yr oscillations are usually attributed to ENSO which is a globally dominating form of variability on annual to decadal time scales (e.g. Kleeman, 2008). It is a broadband phenomenon with several spectral peaks and the highest peak is around 4 yr. This can also be seen in our analysis of 20CR, HadGEM2 and MPI-ESM data sets because most of the significant oscillations are concentrated on 2–6 yr time scales. However, the spectra of MPI-ESM (Fig. 5) differs distinctly from the spectra of the other two data sets: the power on 2–8 yr time scales is distributed on a large set of components (especially ST-PCs 4–18) which also decreases the statistical significance of oscillations on those time scales.
In HadGEM2, significant decadal scale oscillations are identified with all lag window lengths. Dominant peak on the decadal time scales has been noted by Collins et al. (2008) and one of the possible reasons for this is in deficiencies of simulation of the ENSO phase-changing process in HadGEM2 (Martin et al., 2010).
There are also significant multi-decadal components in 20CR data set, but their period decreases with increasing lag window M. The time series to be analysed become shorter with increasing M and this may have an effect on the identified period length. We did not find significant multi-decadal components in HadGEM2 and MPI-ESM data sets, although 27 yr and 26 yr periods are identified on 10% significance levels, but only with a single lag window length. However, the use of a lag window M typically allows only the distinction of oscillations with periods ≤M and thus the interpretation of those multi-decadal components remains uncertain.
The final step of our analysis is to reconstruct the decomposed signals in the original space. As an illustration, we have chosen to reconstruct the signal corresponding to approximately 5.5 yr variability, which was identified in all the data sets.
In order to see the time evolution of the ~5.5 yr variability, we have reconstructed the time series in each gridpoint of the original data set with the ST-PCs corresponding to the signal of interest. I.e., in the reconstruction we have projected the original (centred) data set onto ST-PCs (calculated in the reduced space) to obtain ST-EOFs in the original space and then projected the ST-PCs onto those ST-EOFs (see Appendices A.2 and A.4 for more details). In order to see the global effects of the ~5.5 yr cycle, the time series of each grid point has its original variance. The above calculations were completed for each data set using their own ~5.5 yr patterns. ST-PCs 5 and 6 of the 20CR data set (Fig. 3), ST-PCs 6 and 7 of HadGEM2 (Fig. 4) and ST-PCs 4 and 5 of MPI-ESM (Fig. 5) were used in the reconstruction.
Once we have reconstructed the time series in each gridpoint we can plot the anomalies related to the signal month by month. These plots are presented as animations of each data set (20CR, HadGEM2 and MPI-ESM) for a time period of 1901–2012 (the animations are available at www.youtube.com/channel/UCRjwc6cI-TzbvtShONYZ7cg). In Fig. 8, we also show the global patterns of the ~5.5 yr variability of near-surface temperature anomaly. The patterns are composites of eight cases, when the oscillation is in its positive phase in the equatorial Pacific. Positive events are defined as an average of winter months (November–March).
The temperature anomalies of 20CR have many similarities to global El Niño effects, such as above average temperatures in the central and eastern equatorial Pacific Ocean, in the western and northern parts of North-America and South-America as well as in South-East Asia, Australia and southern Africa. Below average anomalies are found in the south-east parts of North-America, in the north-west and south-west Pacific as well as in northern parts of Eurasia. In 20CR, a typical north-south wave train is also seen, but the east-west patterns are weaker, except for the anomalies at the Amazonas.
HadGEM2 and MPI-ESM show similarities to 20CR, but differences can be seen, for example, in the Pacific forcing patterns. Especially in MPI-ESM the centre of the forcing pattern seems to be more western. In the model simulations, the negative anomaly near the west-coast of North-America extends to the continent, which is not detected in 20CR. The positive anomalies in HadGEM2 and MPI-ESM also extend into the northern Eurasia and there are anomaly patterns in the southern Indian Ocean which are absent in 20CR. MPI-ESM has a stronger positive anomaly in the coast of South-East Asia compared to the other two data sets. In addition, there is a strong anomaly near the Antarctic Peninsula in the Weddel Sea in the 20CR data set which is not detected in the model simulations. The anomaly patterns in the Atlantic Ocean are also weaker in 20CR compared to simulations.
The animations of the ~5.5 yr pattern (available at www.youtube.com/channel/UCRjwc6cI-TzbvtShONYZ7cg) show some more features in addition to the ones seen in Fig. 8. For instance, in 20CR animation there is a quite strong anomaly pattern to the west of Ural Mountains. This pattern is not usually associated with ENSO, and its maximum negative and positive phases seem to occur at different times compared to the ENSO-related anomaly patterns in the Pacific. However, this pattern to the west of Ural might also reflect some other phenomenon, mixed with the ENSO patterns.
The animations also show that the variability has a more propagating character in 20CR data set whereas the anomaly patterns in the model simulations are more stationary. In the northern and southern Pacific Ocean, for example, the anomalies seem to propagate eastward in the 20CR animation.
Compared to 20CR, HadGEM2 and MPI-ESM show a richer structure in Fig. 8 and in the animations. One has to remember that the reanalysis data set is an ensemble mean whereas the analysis of the climate model simulations is conducted on a single ensemble member of each model. This may also contribute to the structure seen in the model simulations. Different, more or less real, phenomena may also be mixed in the variability patterns of the simulations.
We have introduced an RMSSA algorithm, which allows the calculation of MSSA of extremely high-dimensional problems. The RMSSA algorithm first reduces the dimension of the original data set by RP, then decomposes the data set into components of different frequencies by calculating MSSA in a reduced space, and finally reconstructs the components in the original high-dimensional space.
We have applied the RMSSA algorithm to decompose the monthly mean near-surface air temperature of the 20th century reanalysis and the historical 20th century simulations of HadGEM2-ES and MPI-ESM-MR extracted from the CMIP5 data archives. We have also performed Monte-Carlo simulations in order to estimate the significance of the identified low-frequency components. Our analysis shows that 2–6 yr oscillations are present in all the data sets. Their statistical significance is highest in HadGEM2 while in MPI-ESM the power on those timescales is distributed on a large set of components decreasing their statistical significance.
2–6 yr oscillations are usually attributed to ENSO which is a globally dominating form of variability on annual to decadal time scales. Our global monthly animations of 5–6 yr near-surface temperature cycle match quite well with the known temperature anomalies related to ENSO. The reanalysis and the historical simulations have similar anomaly patterns in the central and eastern Pacific Ocean, around the northern part of Indian Ocean as well as in the north-west North-America, but also some notable differences in several areas, such as Eurasia. Also, our animations of the 5–6 yr cycle reveal a propagating structure in the near-surface temperature anomalies of 20CR, while the variability in HadGEM2 and MPI-ESM data sets is more stationary. The focus of this study was to introduce the RMSSA algorithm and the discussion on the possible causes for the differences in oscillatory patterns of the data sets is limited. However, this would be a subject for a further study with a larger set of climate model data sets included.
RMSSA algorithm is a powerful tool when the dimensions of the data sets become prohibitively large. It allows a computationally efficient way of decomposing a data set into its spatio-temporal patterns. Several climatic state variables can be incorporated in the RMSSA at the same time in order to find the co-varying signals and illustrate their propagation. RMSSA can also be used in studying the oscillations in three dimensions including data from several atmospheric levels in the analysis.
We thank Petteri Uotila and Mikko Alestalo from Finnish Meteorological Institute for interesting discussions at earlier phases of the work. We also thank Jouni Räisänen from University of Helsinki for preparing the HadGEM2-ES and MPI-ESM-MR data sets used in this study. In addition, we thank the reviewer for her/his valuable comments and suggestions that helped in improving the manuscript. This research has been funded by the Academy of Finland (project number 140771).
The aim of MSSA is to identify spatially and temporally coherent patterns in a multivariate data set. In MSSA terminology, the columns of the original data matrix X_{N×L} are called channels. In case of gridded data set, N represents the time steps and L is the number of grid points:
The next step is to construct an augmented data matrix A, which contains M lagged copies of each channel in X:
and
In MSSA, M represents the lag window. A has now ML columns and $N\prime =N-M+1$ rows. The singular value decomposition (SVD) of A can now be calculated:
The vectors of U_{A} are the eigenvectors of ${\mathbf{Z}}_{A}=\frac{1}{ML}\mathbf{A}{\mathbf{A}}^{T}$ and ${\mathbf{V}}_{A}^{T}$ contains the eigenvectors of ${\mathbf{C}}_{A}=\frac{1}{N\prime}{\mathbf{A}}^{T}\mathbf{A}$. These vectors are orthogonal and often called space-time principal components (ST-PCs) and space-time empirical orthogonal functions (ST-EOFs), respectively. Diagonal elements of D_{A} are the eigenvalues of C_{A} or Z_{A}.
Optionally the lag-covariance matrix C_{A} (or Z_{A}) and its eigendecomposition can be calculated to yield eigenvectors ${\mathbf{V}}_{A}^{T}$ (or U_{A}) and eigenvalues (diagonal elements of matrix ${\mathbf{D}}_{A}={\mathbf{V}}_{A}^{T}{\mathbf{C}}_{A}{\mathbf{V}}_{A}$or ${\mathbf{D}}_{A}={\mathbf{U}}_{A}^{T}{\mathbf{Z}}_{A}{\mathbf{U}}_{A}$). Matrix U_{A} (or ${\mathbf{V}}_{A}^{T}$) can be obtained by projecting A onto ${\mathbf{V}}_{A}^{T}$(or U_{A}). If $N\prime >ML$ (or $ML>N\prime $), it is more convenient to estimate C_{A} (or Z_{A}) because it is smaller. See Allen and Robertson (1996) for details.
The components obtained by MSSA can be tested against a null-hypothesis of the data being generated by independent AR(1) processes (i.e. red noise). The red noise model has the form:
where γ_{s} is the lag-1 autocorrelation of channel s (in the original data set), ${\alpha}_{s}=\sqrt{{c}_{s}(1-{\gamma}_{s}^{2})}$ (c_{s} is the variance of channel s) and W_{t,s} is Gaussian white noise. The data set generated by the model in (A5) is called the surrogate data set and it is subjected to MSSA in the same way as the original data set. Large number of surrogates are generated in order to estimate the confidence limits for the MSSA results of the original data set.
In the test of Allen and Robertson (1996), the lag-covariance matrices of the original data set and the red-noise surrogates are projected either onto the data-adaptive basis (i.e. U_{A} or ${\mathbf{V}}_{A}^{T}$) or the null-hypothesis basis. The null-hypothesis basis can be calculated from the expected lag-covariance matrix C_{N} of the red-noise surrogates. C_{N} can be estimated analytically by
Projection onto the red-noise basis is considered more reliable because the use of the data-adaptive basis assumes the existence of an oscillation even in a case where it is uncertain whether or not the oscillation is significant.
According to Allen and Robertson (1996), the input channels should be uncorrelated (or at least nearly) at zero lag for the test to be useful. In the case of a gridded data set, where all the grid point time series are used as input channels, the decorrelation condition is not valid. The test might still be useful if we are using grid points sufficiently far from each other as the input channels for the test (Ghil et al., 2002).
1: Original MSSA algorithm enhanced by RP
2: PC-based MSSA algorithm enhanced by RP
The method to back-project from the reduced space to the original space in the case of RP + SVD is explained in Seitola et al. (2014) (Appendix A.1) but we also present it briefly here:
The SVD of the original data matrix X_{N×L} is:
U contains the eigenvectors of Z=XX^{T}.
Random projection (RP) of X is P=XR, where R_{L×k} is the projection matrix and the row vectors of R are scaled to have unit length. Thus, we can write:
In the previous, we have assumed that $\mathbf{R}{\mathbf{R}}^{T}\approx \mathbf{I}$ because the row vectors of R are nearly orthonormal. It is also possible to make the vectors of R strictly orthonormal, in which case $\mathbf{R}{\mathbf{R}}^{T}=\mathbf{I}$. However, orthogonalisation is often not necessary, because the difference between the orthogonalised and non-orthogonalised random vectors is very small, especially in high-dimensions.
Let's rewrite (A7) as ${\mathbf{X}}_{N\times L}={\mathbf{U}}_{N\times r}{\mathbf{D}}_{r\times r}{\mathbf{V}}_{r\times L}^{T}$, where r=rank(X). Now we can manipulate (A7):
Because Z≈Z_{RP} we can approximate
In the previous, we have defined U_{RP} as N×k and D_{RP} as a k×k matrix, where k is the rank of matrix P_{N×k}.
In this appendix, we will explain how to get from the reduced space back to the original space in the case of RP + MSSA.
Let's write the original data matrix X_{N×L} as
where x_{i} are the row vectors of X.
The augmented matrix A of X is already defined in Appendix A.1. Now let's calculate AA^{T}.
After some algebra we get
Now let's calculate RP of X:
The augmented matrix of is A_{RP}:
Let's calculate ${\mathbf{A}}_{RP}{\mathbf{A}}_{RP}^{T}$:
Because $\mathbf{R}{\mathbf{R}}^{T}\approx \mathbf{I}$, the first element of ${\mathbf{A}}_{RP}{\mathbf{A}}_{RP}^{T}$ can be written as ${\mathbf{x}}_{\mathbf{1}}\mathbf{R}{\mathbf{R}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{1}}^{\mathbf{T}}+{\mathbf{x}}_{\mathbf{2}}\mathbf{R}{\mathbf{R}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{2}}^{\mathbf{T}}+\cdots +{\mathbf{x}}_{\mathbf{M}}\mathbf{R}{\mathbf{R}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{M}}^{\mathbf{T}}\approx {\mathbf{x}}_{\mathbf{1}}{\mathbf{x}}_{\mathbf{1}}^{\mathbf{T}}+{\mathbf{x}}_{\mathbf{2}}{\mathbf{x}}_{\mathbf{2}}^{\mathbf{T}}+\cdots +{\mathbf{x}}_{\mathbf{M}}{\mathbf{x}}_{\mathbf{M}}^{\mathbf{T}}$${\mathbf{x}}_{\mathbf{1}}\mathbf{R}{\mathbf{R}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{1}}^{\mathbf{T}}+{\mathbf{x}}_{\mathbf{2}}\mathbf{R}{\mathbf{R}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{2}}^{\mathbf{T}}+\cdots +{\mathbf{x}}_{\mathbf{M}}\mathbf{R}{\mathbf{R}}^{\mathbf{T}}{\mathbf{x}}_{\mathbf{M}}^{\mathbf{T}}\approx {\mathbf{x}}_{\mathbf{1}}{\mathbf{x}}_{\mathbf{1}}^{\mathbf{T}}+{\mathbf{x}}_{\mathbf{2}}{\mathbf{x}}_{\mathbf{2}}^{\mathbf{T}}+\cdots +{\mathbf{x}}_{\mathbf{M}}{\mathbf{x}}_{\mathbf{M}}^{\mathbf{T}}$
After calculating all the elements of ${\mathbf{A}}_{RP}{\mathbf{A}}_{RP}^{T}$ as above, we see that $\mathbf{A}{\mathbf{A}}^{T}\approx {\mathbf{A}}_{RP}{\mathbf{A}}_{RP}^{T}$. Therefore, as in Appendix A.3, we can approximate
Same kind of reasoning applies also when the PCs of the data set are used as channels in MSSA. We can write PCs as ${\mathbf{U}}_{N\times r}{\mathbf{D}}_{r\times r}={\mathbf{X}}_{N\times L}{\mathbf{V}}_{L\times r}$, where r=rank(X). Vectors of V are orthonormal, so in the above calculations we can replace R with V.
Allen M. R. , Robertson A. W . Distinguishing modulated oscillations from coloured noise in multivariate datasets . Clim. Dyn . 1996 ; 12 ( 11 ): 775 – 784 .
Bingham E. , Mannila H . Random projection in dimensionality reduction: applications to image and text data . Proceedings of the Seventh ACM SIGKDD International Conference on Knowledge Discovery and Data Mining . 2001 ; ACM : New York . 245 – 250 . KDD '01 .
Broomhead D. S. , King G. P . Extracting qualitative dynamics from experimental data . Physica D . 1986a ; 20 : 217 – 236 .
Broomhead D. S. , King G. P . Sarkar S . On the qualitative analysis of experimental dynamical systems . Nonlinear Phenomena and Chaos . 1986b ; Bristol : Adam Hilger . 113 – 144 .
Chiu S. K . Coherent and random noise attenuation via multichannel singular spectrum analysis in the randomized domain . Geophys. Prospect . 2013 ; 61 : 1 – 9 .
Cleveland R. B. , Cleveland W. S. , McRae J. E. , Terpenning I . STL: a seasonal-trend decomposition procedure based on loess . J Off. Stat . 1990 ; 6 : 373 .
Collins W. J. , Bellouin N. , Doutriaux-Boucher M. , Gedney N. , Halloran P. , co-authors . Development and evaluation of an Earth-System model HadGEM2 . Geosci. Model Dev . 2011 ; 4 : 1051 – 1075 .
Collins W. J. , Bellouin N. , Doutriaux-Boucher M. , Gedney N. , Hinton T. , co-authors . Evaluation of the HadGEM2 model. Met Office Hadley Centre Technical Note no. HCTN 74 . 2008 . available from Met Office, FitzRoy Road, Exeter EX1 3PB. Online at: http://www.metoffice.gov.uk/publications/HCTN/index.html .
Compo G. P. , Whitaker J. S. , Sardeshmukh P. D. , Matsui N. , Allan R. J. , co-authors . The twentieth century reanalysis project . Quart. J. Roy. Meteorol. Soc . 2011 ; 137 : 1 – 28 .
Dasgupta S. , Gupta A . An elementary proof of a theorem of Johnson and Lindenstrauss . Rand. Struct. Algo . 2003 ; 22 : 60 – 65 .
Elsner J. B. , Tsonis A. A . Singular Spectrum Analysis: A New Tool in Time Series Analysis . 1996 ; New York, NY, USA : Springer Science & Business Media .
Frankl P. , Maehara H . The Johnson-Lindenstrauss lemma and the sphericity of some graphs . J. Combin. Theory Ser. B . 1988 ; 44 : 355 – 362 .
Ghil M. , Allen M. R. , Dettinger M. D. , Ide K. , Kondrashov D. , co-authors . Advanced spectral methods for climatic time series . Rev. Geophys . 2002 ; 40 ( 1 ): 1 – 41 .
Halko N. , Martinsson P. G. , Tropp J. A . Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions . SIAM Rev . 2011 ; 53 ( 2 ): 217 – 288 .
Kawale J. , Chatterjee S. , Kumar A. , Liess S. , Steinbach M. , co-authors . Anomaly construction in climate data: issues and challenges . 2011 ; 189 – 203 . Proceedings of the 2011 Conference on Intelligent Data Understanding, CIDU, Mountain View, CA, USA .
Kleeman R . Stochastic theories for the irregularity of ENSO . Philos. Trans. Roy. Soc. A Math. Phys. Eng. Sci . 2008 ; 366 ( 1875 ): 2509 – 2524 .
Mann M. E. , Lees J. M . Robust estimation of background noise and signal detection in climatic time series . Clim. Change . 1996 ; 33 : 409 – 445 .
Martin G. M. , Milton S. F. , Senior C. A. , Brooks M. E. , Ineson S. , co-authors . Analysis and reduction of systematic errors through a seamless approach to modeling weather and climate . J. Clim . 2010 ; 23 ( 22 ): 5933 – 5957 .
Meehl G. A. , Goddard L. , Boer G. , Burgman R. , Branstator G. , co-authors . Decadal climate prediction: an update from the trenches . Bull. Amer. Meteor. Soc . 2014 ; 95 : 243 – 267 .
Meehl G. A. , Goddard L. , Murphy J. , Stouffer R. J. , Boer G. , co-authors . Decadal prediction: can it be skillful? . Bull. Amer. Meteor. Soc . 2009 ; 90 : 1467 – 1485 .
Moron V. , Robertson A. W. , Ghil M . Impact of the modulated annual cycle and intraseasonal oscillation on daily-to-interannual rainfall variability across monsoonal India . Clim. Dyn . 2012 ; 38 : 2409 – 2435 .
Oropeza V. , Sacchi M . Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis . Geophysics . 2011 ; 76 : V25 – V32 .
Plaut G. , Vautard R . Spells of low-frequency oscillations and weather regimes in the Northern Hemisphere . J. Atmos. Sci . 1994 ; 51 ( 2 ): 210 – 236 .
Rinne J. , Karhila V . Empirical orthogonal functions of 500 mb height in the northern hemisphere determined from a large data sample . Quart. J. Roy. Meteorol. Soc . 1979 ; 105 : 873 – 884 .
Seitola T. , Mikkola V. , Silen J. , Järvinen H . Random projections in reducing the dimensionality of climate simulation data . Tellus A . 2014 ; 66 25274. DOI: http://dx.doi.org/10.3402/tellusa.v66.25274 .
Solomon A. , Goddard L. , Kumar A. , Carton J. , Deser C. , co-authors . Distinguishing the roles of natural and anthropogenically forced decadal climate variability . Bull. Amer. Meteor. Soc . 2011 ; 92 : 141 – 156 .
Stevens B. , Giorgetta M. , Esch M. , Mauritsen T. , Crueger T. , co-authors . Atmospheric component of the MPIM Earth System Model: ECHAM6 . J. Adv. Model. Earth Syst . 2013 ; 5 ( 2 ): 146 – 172 .
Taylor K. E. , Stouffer R. J. , Meehl G. A . An overview of CMIP5 and the experiment design . Bull. Amer. Meteor. Soc . 2012 ; 93 : 485 – 498 .
Thomson D. J . Spectrum estimation and harmonic analysis . Proc. IEEE . 1982 ; 70 : 1055 – 1096 .
Vautard R. , Ghil M . Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series . Phys. D Nonlin. Phenom . 1989 ; 35 ( 3 ): 395 – 424 .
Von Storch H. , Zwiers F. W . Statistical Analysis in Climate Research . 1999 ; Cambridge : Cambridge University Press .
Weare B. C. , Nasstrom J. S . Examples of extended empirical orthogonal function analysis . Mon. Weath. Rev . 1982 ; 110 : 481 – 485 .