Nonlinear multichannel singular spectrum analysis of the tropical Pacific climate variability using a neural network approach William W. Hsieh and Aiming Wu Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, British Columbia, Canada Received 4 May 2001; revised 3 January 2001; accepted 11 January 2001; published 18 July 2002. [1] Singular spectrum analysis (SSA), a linear (univariate and multivariate) time series technique, performs principal component analysis (PCA) on an augmented data set containing the original data and time-lagged copies of the data. Neural network theory has meanwhile allowed PCA to be generalized to nonlinear PCA (NLPCA). In this paper, NLPCA is further extended to perform nonlinear SSA (NLSSA): First, SSA is applied to reduce the dimension of the data set; the leading principal components (PCs) of the SSA then become inputs to an NLPCA network (with a circular node at the bottleneck). This network performs the NLSSA by nonlinearly combining all the input SSA PCs. The NLSSA is applied to the tropical Pacific sea surface temperature anomaly (SSTA) field and to the sea level pressure anomaly (SLPA) field for the 1950–2000 period. Unlike SSA modes, which display warm and cool periods of similar duration and intensity, NLSSA mode 1 shows the warm periods to be shorter and more intense than the cool periods, as observed for the El Niño-Southern Oscillation. Also, in SSTA NLSSA mode 1 the peak warm event is centered in the eastern equatorial Pacific, while the peak cool event is located around the central equatorial Pacific, an asymmetry not found in the individual SSA modes. A quasi-triennial wave of about a 39 month period is found in NLSSA mode 2 of the SSTA and of the SLPA. INDEX TERMS: 3339 Meteorology and Atmospheric Dynamics: Ocean/atmosphere interactions (0312, 4504); 3220 Mathematical Geophysics: Nonlinear dynamics; 4215 Oceanography: General: Climate and interannual variability (3309); 4522 Oceanography: Physical: El Niño; KEYWORDS: El Niño, Southern Oscillation, neural networks, singular spectrum analysis, tropical Pacific 1. Introduction [2] Principal component analysis (PCA), also known as empirical orthogonal function (EOF) analysis, is a classical multivariate statistical technique used to analyze a set of variables {xi} [von Storch and Zwiers, 1999]. It is com- monly used to reduce the dimensionality of the data set {xi} and to extract features (or recognize patterns). [3] By the 1980s, interests in chaos theory and dynam- ical systems led to further extension of PCA to singular spectrum analysis (SSA): The information contained in a continuous variable and its first to (L  1)th derivatives can be approximated by a discrete time series and the same time series lagged by 1,. . ., L  1 time steps [Elsner and Tsonis, 1996, chap. 4]. Thus a given time series and its lagged versions can be regarded as a set of variables {xi} (i = 1,. . ., L), an augmented data set, which can be analyzed by the PCA. This resulting method is the SSA with window L. In the multivariate case where there is more than one time series, one can again make lagged copies of the time series, treat the lagged copies as extra variables, and apply the PCA to this augmented data set, resulting in the multichannel SSA (MSSA) method, also called the space-time PCA (ST-PCA) method, or the extended EOF (EEOF) method (though in typical EEOF applications, only a small number of lags are used). For brevity we will use the term SSA to denote both SSA and MSSA. A recent review of the SSA method is given by Ghil et al. [2002]. [4] Neural network (NN) models, which first became popular in the late 1980s, have significantly advanced nonlinear empirical modeling. Each member of the hierachy of classical multivariate methods (multiple linear regression, PCA, and canonical correlation analysis (CCA)) has been nonlinearly generalized by NN models: nonlinear multiple regression [Rumelhart et al., 1986], nonlinear PCA (NLPCA) [Kramer, 1991], and nonlinear CCA (NLCCA) [Hsieh, 2000]. The tropical Pacific climate variability has been studied by the NLPCA method [Mon- ahan, 2001; Hsieh, 2001a] and by the NLCCA method [Hsieh, 2001b]. [5] In PCA a straight line approximation to the data set is sought that accounts for the maximum amount of variance in the data. In NLPCA the straight line is replaced by an open continuous curve for approximating the data [Kramer, 1991]. Kirby and Miranda [1996] introduced a NLPCA with a circular node at the network bottleneck (henceforth NLPCA.cir), so that the nonlinear principal component (NLPC) as represented by the circular node is an angular variable q and the NLPCA.cir is capable of approximating the data by a closed continuous curve. JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 107, NO. C7, 3076, 10.1029/2001JC000957, 2002 Copyright 2002 by the American Geophysical Union. 0148-0227/02/2001JC000957$09.00 13 - 1 [6] This paper belongs to a series of papers on nonlinear PCA and its extensions. In Hsieh [2001a] both NLPCA and NLPCA.cir were applied to study the tropical Pacific sea surface temperature anomalies (SSTA). W. W. Hsieh and A. Wu (Nonlinear singular spectrum analysis by neural networks, submitted to Neural Networks, 2002, hereinafter referred to as Hsieh and Wu, submitted manuscript, 2002). developed the nonlinear SSA (NLSSA) method from the NLPCA.cir network, resulting in a new nonlinear time series technique. This nonlinear spectral technique allows the detection of highly anharmonic oscillations, as was illustrated by a stretched square wave imbedded in white noise, which showed NLSSA to be superior to SSA and classical Fourier spectral analysis (Hsieh and Wu, submitted manuscript, 2002). In this present paper the NLSSA method is used to analyze the tropical Pacific SSTA field, and the sea level pressure anomaly (SLPA) field, for the period from 1950 to 2000. 2. Theory: From NLPCA to NLSSA [7] Since NLPCA.cir has been explained fully by Hsieh [2001a], we only outline its main properties here. The input data are in the form x(t) = [x1,. . .,xl], where each variable xi(i = 1,. . .,l ) is a time series containing n observations. The model (Figure 1) is a standard feed forward NN (i.e., multilayer perceptron), with 3 ‘‘hidden’’ layers of variables or ‘‘neurons’’ (denoted by circles) sandwiched between the input layer x on the left and the output layer x0 on the right. Next to the input layer (with l neurons) is the encoding layer (with m neurons), followed by the ‘‘bottleneck’’ layer, then the decoding layer (with m neurons), and finally, the output layer (with l neurons). The jth neuron vj (k) in the kth layer receives its value from the neurons v(k1) in the preceding layer; that is, vj (k) = fk(wj (k)  v(k1) + bj(k) ), where wj(k) is a vector of weight parameters, bj (k) is a bias parameter, and the transfer functions f1 and f3 are the hyperbolic tangent functions, while f2 and f4 are simply the identity functions. Hence a total of four successive layers of transfer functions are needed to map from the inputs x to the outputs x0. In NLPCA.cir the bottleneck contains two neurons p and q confined to lie on a unit circle, i.e., only 1 degree of freedom, as represented by the angle q. Effectively, a non- linear function q = F(x) maps from the higher dimension input space to the one-dimensional (1-D) bottleneck space, followed by an inverse transform x0 = G(q) mapping from the bottleneck space back to the original space, as repre- sented by the outputs. To make the outputs as close to the inputs as possible, the cost function J = hkx  x0k2i (i.e., the mean square error (MSE)) is minimized (where hi denotes a sample or time mean). Through the optimization the values of the weight and bias parameters are solved. Data compression is achieved by the bottleneck, yielding the nonlinear principal component q. [8] Hsieh [2001a] noted that the NLPCA.cir, with its ability to extract closed curve solutions, is particularly ideal for extracting periodic or wave modes in the data. In SSA it is common to encounter periodic modes, each of which had to be split into a pair of modes [Elsner and Tsonis, 1996], as the underlying PCA technique is not capable of modeling a periodic mode (a closed curve) by a single mode (a straight line). Thus, two SSA modes can easily be combined into one NLPCA.cir mode. In principle, the original variables and their lagged versions can be input to the NLPCA.cir to extract the NLSSA solution. [9] The problem is that usually one cannot afford to use many lags before the large number of input variables to the network results in more model parameters than samples (i.e., measurements in time). Even in NLPCA or NLPCA.cir, the number of input variables can be so large that there are more model parameters than samples. To avoid this sit- uation, the data are usually first condensed by the classical PCA method where only the first few leading PCs (i.e., the time coefficients from the PCA) are retained, resulting in far fewer input variables to the NLPCA or NLPCA.cir network. [10] An analogous approach to reduce greatly the number of input variables to the network can be used for the NLSSA, except that instead of PCA, SSA is used to prefilter the data. First, the original data are analyzed by the SSA with window L. Only the first few leading SSA modes are retained, and their PCs, [x1,. . .,xl], are then served as input variables to the NLPCA.cir network. The NLPCA.cir finds a continuous curve solution by nonli- nearly relating the PCs, thereby giving NLSSA mode 1. Hsieh [2001a] pointed out that the general configuration of the NLPCA.cir can model not only closed curve solutions but also open curve solutions like the Kramer [1991] NLPCA. Because the NLPCA.cir is more general than the NLPCA (with one bottleneck neuron), it is used here to perform the NLSSA. More details of the NLSSA theory are given by Hsieh and Wu (submitted manuscript, 2002). [11] A small amount of weight penalty was added to the cost function as by Hsieh [2001a] and Hsieh and Wu (submitted manuscript, 2002), with the penalty parameter P = 0.02. Because of local minima in the cost function, there is no guarantee that the optimization algorithm reaches the global minimum. Hence an ensemble of 30 NNs with random initial weights and bias parameters was run. Also, 20% of the data were randomly selected as validation data and withheld from the training of the NNs. Runs where the MSE was larger for the validation data set than for the Figure 1. A schematic diagram of the NN model for calculating the NLPCA with a circular node at the bottleneck (NLPCA.cir). 13 - 2 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS Figure 2. The SSA (i.e., ST-PCA or EEOF) modes (a) 1, (b) 2, (c) 3, (d) 4, (e) 5, and (f ) 6 for the tropical Pacific SSTA. The contour plots display the space-time eigenvectors (loading patterns), showing the SSTA along the equator as a function of the lag. Solid contours indicate positive anomalies, and dashed contours indicate negative anomalies, with the zero contour indicated by the thick solid curve. In a separate panel beneath each contour plot, the PC of each SSA mode is also plotted as a time series (where each tick mark on the abscissa indicates the start of a year). The time of the PC is synchronized to the lag time of 0 month in the space-time eigenvector. HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 3 training data set were rejected to avoid overfitted solu- tions. Then the NN with the smallest overall MSE (i.e., the MSE computed over both training and validation data) was selected as the solution. Separate runs were made using m = 2,. . .,6. The overall MSE dropped with increasing m but eventually levelled off. Following the principle of parsimony, we chose the smallest m when the levelling occurred as the solution. For the SSTA NLSSA modes, m = 4, while for SLPA modes, m = 6. 3. NLSSA of the Tropical Pacific SSTA [12] The Smith et al. [1996] monthly SST for the domain 21S–21N, 123E–69W during January 1950 to Decem- ber 2000 were used. The original 2  2 data were averaged to 2 latitude  4 longitude resolution (to reduce the memory requirement in the subsequent SSA calcula- tion), then the climatological seasonal cycle was subtracted and the linear trend was removed from the SST to yield the SSTA data. [13] As we want to resolve the El Niño-Southern Oscillation (ENSO) variability, a window of 73 months was chosen. With a lag interval of 3 months, the original plus 24 lagged copies of the SSTA data formed the augmented SSTA data set. The first eight SSA modes explain 12.4, 11.7, 7.1, 6.7, 5.4, 4.4, 3.5, and 2.8%, respectively, of the total variance of the augmented data set (with the first six modes shown in Figure 2). The first two modes have space-time eigenvectors (i.e., loading patterns) showing an oscillatory timescale of about 48 months, comparable to the ENSO timescale, with the mode 1 anomaly pattern occurring about 12 months before a very similar mode 2 pattern; that is, the two patterns are in quadrature. The PC time series also show similar timescales for modes 1 and 2. Modes 3 and 5 show longer timescale fluctuations, while modes 4 and 6 show shorter timescale fluctuations, around the 30 month time- scale. [14] With the eight PCs as input x1,. . ., x8 to the NLPCA.cir network the resulting NLSSA mode 1 is a closed curve in the 8-D PC space and is plotted in the x1-x2-x3 space in Figure 3. NLSSA mode 1 is basically a large loop aligned parallel to the x1-x2 plane, thereby combin- ing the first two SSA modes. The solution also shows some modest variations in the x3 direction. This NLSSA mode 1 explains 24.0% of the total variance of the augmented data set, essentially that explained by the first two SSA modes together. The linear PCA solution is shown as a straight line in Figure 3, which is, of course, simply SSA mode 1. Of interest is r, the ratio of the MSE of the nonlinear solution to the MSE of the linear solution. Here r = 0.71. [15] The NLSSA mode 1 space-time loading pattern for a given value of q can be obtained by mapping from q to the outputs x0, which are the eight PC values corresponding to the given q. Multiplying each PC value by its corresponding SSA eigenvector and summing over the eight modes, we obtain the NLSSA mode 1 pattern corresponding to the given q. [16] The NLSSA mode 1 loading patterns for various q values are shown in Figure 4. Comparing the patterns in Figure 4 with the patterns from the first two SSA modes in Figure 2, we find three notable differences: (1) The presence of warm anomalies for 24 months followed by cool anomalies for 24 months in the first two SSA modes is replaced by warm anomalies for 18 months followed by cool anomalies for about 33 months in NLSSA mode 1. Although the cool anomalies can be quite mild for long periods, they can develop into full La Niña cool events (Figure 4c). (2) The El Niño warm events are strongest near the eastern boundary, while the La Niña events are strongest near the central equatorial Pacific in NLSSA mode 1, an asymmetry not found in the individual SSA modes. (3) The magnitude of the peak positive anomalies is significantly larger than that of the peak negative anomalies in NLSSA mode 1 (Figure 4c), again an asymmetry not found in the individual SSA modes. All three differences indicate that NLSSA mode 1 is much closer to the observed ENSO properties than the first two SSA modes are. [17] The NLPC q of NLSSA mode 1, plotted as a time series (cyclically bounded between p and p radians) in Figure 5, reveals q generally increasing with time. Thus the patterns in Figure 4 generally evolve from Figures 4a–4f with time. [18] Next, we try to reconstruct the SSTA field from NLSSA mode 1 for the whole record. The PC at a given time gives not only the SSTA field at that time but a total of 25 SSTA fields (spread 3 months apart) covering the 73 months in the lag window. Thus (except for the final 72 months of the record) the SSTA field at any time is taken to be the average of the 25 SSTA fields produced by the PC of the current month and the Figure 3. NLSSA mode 1 for the tropical Pacific SSTA. The PCs of SSA modes 1–8 were used as inputs x1,. . ., x8 to the NLPCA.cir network, with the resulting NLSSA mode 1 shown as (densely overlapping) crosses in the x1–x2–x3 3-D space. The projections of this mode onto the x1–x2, x1–x3, and x2–x3 planes are denoted by the (densely overlapping) small circles, and the projected data are denoted by dots. For comparison, linear SSA mode 1 is shown by the dashed line in the 3-D space and by the projected solid lines on the 2-D planes. 13 - 4 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS preceding 72 months. The SSTA field reconstructed from NLSSA mode 1 (called the reconstructed component (NLRC1)) during 1971–2000 is shown in Figure 6, which compares well with the observed field, except for the weaker amplitude (Figure 7). [19] To extract NLSSA mode 2, we removed the NLSSA mode 1 x0 from the data x and input the residuals (x  x0) into the same NLPCA.cir network. Since this NLSSA mode 2 mainly combines SSA modes 3 and 4, it is plotted in the x2–x3–x4 space (Figure 8), where the most prominent feature is the large loop in the x3–x4 plane. There are also nonlinear relations manifested in the x2–x3 and x2–x4 planes, as well as in the x1–x3 and x1–x4 planes (not shown). NLSSA mode 2 explains 12.6% of the total variance of the residuals. In contrast, if the residuals were input into a PCA model, the leading mode (shown by the straight lines in Figure 8) would only explain 7.3% of the variance. For mode 2, r = 0.77. [20] The space-time loading patterns associated with this NLSSA mode 2 (Figure 9) reveal oscillations at around the triennial period. The corresponding NLPC q is also plotted in Figure 5. NLRC2, the SSTA field reconstructed from Figure 4. The SSTA NLSSA mode 1 space-time loading patterns for (a) q = 0, (b) q = 60, (c) q = 120, (d) q = 180, (e) q = 240, and (f ) q = 300. The contour plots display the SSTA anomalies along the equator as a function of the lag time. The contour interval is 0.2C. HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 5 NLSSA mode 2 during 1971–2000 (Figure 10), manifests oscillations at around the 39 month period. [21] What is striking about this NLSSA mode 2 when compared with NLSSA mode 1 is its great regularity. Here the warm and cool anomaly patterns are approximately mirror images of each other, unlike the large differences (in duration and intensity) between warm and cool anomaly patterns found in NLSSA mode 1 (Figure 6). The oscilla- tions at around the 39 month period, except for a few years when the oscillations faded away, are far more regular (Figure 10) than the fluctuations generated by NLSSA mode 1 (Figure 6). The great regularity strongly suggests that this quasi-triennial oscillation of NLSSA mode 2 is a much more linear phenomenon than NLSSA mode 1. [22] We next plot the spatial anomaly patterns during peak warming and during peak cooling in the Nino 3 region (150–90W, 5S–5N) in the eastern equatorial Pacific. Searching over all lags and q from 0 to 360, we found the lag and q when the maximum value of the loading in the Nino 3 region occurred, as well as for the minimum value. The spatial anomaly patterns for NLSSA modes 1 and 2 during maximum and minimum SSTA in the Nino 3 region are shown in Figure 11. For NLSSA mode 1 the cool pattern (Figure 11b) is centered much farther offshore from Peru than the warm pattern (Figure 11a), while for NLSSA mode 2 such a difference is not found between the warm pattern (Figure 11c) and the cool pattern (Figure 11d). The anomalies are more closely confined to the equator for mode 1 than those for mode 2. [23] For the period 1971–2000, comparing the recon- structed SSTA from NLSSA mode 1 (Figure 6), that from NLSSA mode 2 (Figure 10), and the observed SSTA (Figure 7), we note that during the major El Niño events (e.g., 1972, 1982–1983, and 1997) both NLSSA mode 1 and mode 2 developed warm anomalies, while during La Niña events (e.g., 1973–1974 and 1988) both NLSSA mode 1 and mode 2 developed cool anomalies. Intrigu- ingly, in Figures 6, 7, and 10 the positive and negative anomalies in Figures 6a, 7a, and 10a for 1971–1985 match quite well with the anomalies in Figures 6b, 7b, and 10b for 1986–2000, as though history nearly repeated itself after 15 years. The sensitivity of these results to the choice of the window L = 73 months was checked by repeating the calculations with L = 97 months, which yielded essentially the same results. 4. NLSSA of the Tropical Pacific SLPA [24] The tropical Pacific monthly SLP data from the Comprehensive Ocean-Atmosphere Data Set [Woodruff et al., 1987] for the domain 21S–21N, 121E–71W during January 1950 to December 2000 were used. The 2  2 resolution SLP data were interpolated for missing values and averaged into 2 latitude  4 longitude gridded data. The seasonal cycle and the linear trend were then removed from the SLP to give the SLPA data. [25] The first eight SSA modes of the SLPA accounted for 7.9, 7.1, 5.0, 4.9, 4.0, 3.1, 2.5, and 1.9%, respectively, of the total variance of the augmented data. The first two modes displayed the Southern Oscillation (SO), the east- west seesaw of SLPA at around the 50 month period, while the higher modes displayed fluctuations at around the quasi-biennial oscillation (QBO) [Hamilton, 1998] average period of 28 months (not shown). [26] The eight leading PCs of the SSA were then used as inputs, x1,. . .,x8, to the NLPCA.cir network, yielding NLSSA mode 1 for the SLPA. This mode accounts for 17.1% of the variance of the augmented data, significantly more than the variance explained by the first two SSAmodes (15.0%). This is not surprising as the NLSSA mode did more than just combine SSA modes 1 and 2; it also nonlinearly connects SSA mode 3 to SSA modes 1 and 2 (Figure 12). In Figure 5. The NLPC q time series for SSTA NLSSA mode 1 (bottom curve) and mode 2 (second curve from bottom) and for the SLPA NLSSA mode 1 (second curve from top) and mode 2 (top curve). The time series have been vertically displaced by multiples of 2p radians for better visualization. 13 - 6 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS the x1-x3 plane the bowl-shaped projected solution implies that PC3 tends to be positive when PC1 takes on either large positive or large negative values. Similarly, in the x2-x3 plane the hill-shaped projected solution indicates that PC3 tends to be negative when PC2 takes on large positive or negative values. These curves reveal nonlinear interactions between the longer timescale SSA modes 1 and 2, and the shorter timescale SSA mode 3. Here r = 0.68, smaller than that Figure 6. The reconstructed SSTA of NLSSA mode 1 from January 1971 to December 2000. The contour plots display the SSTA along the equator as a function of time, with a contour interval of 0.2C. HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 7 found for SSTA mode 1, indicating somewhat stronger nonlinearity in the SLPA than in the SSTA. [27] The NLPC q (Figure 5) shows q generally increasing with time. The space-time loading patterns for NLSSA mode 1 at various q reveal that the negative phase of the SO is much shorter and more intense than the positive phase (not shown), in agreement with observations and in contrast to SSA modes 1 and 2, where the negative and positive Figure 7. The observed SSTA along the equator from January 1971 to December 2000. The data have been low-pass filtered, with periods of 12 months or less removed. Note, the contour interval is 0.5C instead of the 0.2C in Figure 6. 13 - 8 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS phases of the SO are about equal in duration and magnitude. The SLPA field reconstructed from NLSSA mode 1 (NLRC1) during 1971–2000 is shown in Figure 13, which also reveals the negative phase of the SO to be more intense but of shorter duration than the positive phase. [28] After NLSSA mode 1 had been removed from the data, the residual data were then input into the NLPCA.cir network to extract NLSSA mode 2. This mode is not shown in the PC space as it resembles Figure 8 for SSTA NLSSA mode 2. SLPA NLSSA mode 2 accounts for 8.2% of the variance, versus 4.8% accounted for by the leading PCA mode of the residual data. Here r = 0.77, the same as that for SSTA mode 2. [29] The space-time loading patterns of NLSSA mode 2 reveal east-west seesaw oscillations in the SLPA, but at around the 39 month period (not shown). The SLPA field reconstructed from NLSSA mode 2 (NLRC2) during 1971– 2000 (Figure 14) also manifests oscillations at around the 39 month period. [30] Again, comparing SLPA NLSSA mode 2 (Figure 14) with SLPA NLSSA mode 1 (Figure 13), we note that the positive and negative anomaly patterns of mode 2 are approximately mirror images of each other, in contrast to the large differences (in duration and intensity) between the positive and negative anomaly patterns found in NLSSA mode 1. The oscillations at around the 39 month period, except for a few years when the oscillations faded, are far more regular (Figure 14) than the fluctuations generated by NLSSA mode 1 (Figure 13). The great regularity strongly suggests that this quasi-triennial oscillation of NLSSA mode 2 is a much more linear phenomenon than NLSSA mode 1. The quasi-triennial oscillations in the SLPA (Figure 14) coincide very well with those in the SSTA (Figure 10). [31] The spatial SLPA patterns for NLSSA modes 1 and 2 during maximum and minimum SLPA in the Nino 3 region were also calculated. Again, the strongest SLPA in NLSSA mode 2 tend to occur farther away from the equator than the strongest SLPA in NLSSA mode 1, analogous to what was found for the SSTA (Figure 11). [32] For the period 1971–2000, comparing NLRC1 (Figure 13) and NLRC2 (Figure 14), we note that during the major El Niño events (e.g., 1972, 1982–1983, and 1997) both NLSSA mode 1 and mode 2 developed negative SLPA in the eastern equatorial Pacific (and positive anomalies in the west), while during La Niña events (e.g., 1973–1974 and 1988) both NLSSA mode 1 and mode 2 developed positive SLPA in the eastern equa- torial Pacific. Again, the positive and negative anomalies in Figures 13a and 14a for 1971–1985 match quite well with the anomalies in Figure 13b and 14b for 1986–2000, as though history nearly repeated itself after 15 years. [33] To explore relations between the NLRC1 for the SSTA and that for the SLPA, the NLRC1 of the SSTA averaged over the Nino 3 region was lagged correlated with that of the SLPA over Nino 3. The peak correlation of 0.56 occurred when SLPA leads SSTA by 1 month. In the Nino 3.4 region (170–120W, 5S–5N) the peak correlation was 0.54 with SLPA leading SSTA by 2 months. [34] Between the NLRC2 of the SSTA and the NLRC2 of the SLPA the peak correlation in Nino 3 was 0.96 with SLPA leading SSTA by 1 month, while in Nino 3.4 the peak correlation was 0.92 with SLPA leading by 2 months. The peak correlations here are much higher than those involving the NLRC1, suggesting that the coupling between NLSSA mode 2 of the SSTA and that of the SLPA may be stronger than the NLSSA mode 1 coupling. 5. Summary and Discussion [35] As NN modeling has generalized PCA to NLPCA, it is natural to similarly extend SSA to NLSSA. SSA is simply PCA applied to an augmented data set (containing the original data set and copies of the data set lagged by a range of time steps). In most NLPCA applications the data set is first condensed by the PCA method, and the first few PCs are used as inputs to the NLPCA. Similarly, with NLSSA the data set is first condensed by the SSA, and the first few PCs from the SSA are chosen as inputs to the NLPCA.cir network (the NLPCAwith a circular node at the bottleneck), which extracts the NLSSA mode by nonli- nearly combining the various SSA modes. [36] In general, NLSSA has several advantages over SSA: (1) The PCs from different SSA modes are linearly uncorre- lated; however, they may have nonlinear relationships that can be detected by the NLSSA. (2) Although the SSA modes are not restricted to sinusoidal oscillations in time like the Fourier spectral components, in practice, they are inefficient in modeling strongly anharmonic signals (e.g., the stretched square wave by Hsieh and Wu (submitted manuscript, 2002)), scattering the signal energy into many SSA modes. The NLSSA recombines the SSA modes to extract the anharmonic signal. (3) As different SSA modes are associ- ated with different timescales (e.g., timescales of the ENSO, QBO, and decadal oscillations), the nonlinear relations found by the NLSSA reveal the timescales among which Figure 8. NLSSA mode 2 for the tropical Pacific SSTA shown in the x2–x3–x4 3-D PC space. After removing NLSSA mode 1 from the data we input the residuals into the same NLPCA.cir network to extract NLSSA mode 2. The dots indicate the residuals projected onto the 2-D planes, and the straight lines indicate the PCA approximation to the residuals. HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 9 Figure 9. The SSTA NLSSA mode 2 space-time loading patterns for (a) q = 0, (b) q = 90, (c) q = 180, and (d) q = 270. The contour plots display the SSTA along the equator as a function of the lag time. The contour interval is 0.1C. 13 - 10 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS there are nonlinear relations, thereby disclosing nonlinear relations between seemingly separate phenomena. [37] For the tropical Pacific SSTA the NLSSA mode 1 combines mainly the SSA modes 1 and 2, characterizing the ENSO phenomenon. However, in the SSA modes the warm and cool periods are of similar duration and inten- sity, whereas in the NLSSA mode the warm period is of shorter duration and greater intensity. Furthermore, in the Figure 10. The NLSSA mode 2 reconstructed SSTA along the equator from January 1971 to December 2000. The contour interval is 0.1C. HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 11 NLSSA mode the peak warm event is centered in the eastern equatorial Pacific, while the peak cool event is centered much farther west, a contrast not found in SSA modes 1 and 2. In short, the NLSSA mode gives a much more accurate picture of the ENSO phenomenon than the leading pair of SSA modes. [38] For the tropical Pacific SLPA, NLSSA mode 1 not only combines SSA modes 1 and 2 with ENSO timescales but also reveals nonlinear interactions between these two and mode 3, which has a QBO timescale. Again, the symmetry between the positive and negative phases of the SO found in SSA modes 1 and 2 disappeared in NLSSA mode 1, in agreement with observations. [39] An interesting comparison can be made between NLSSA mode 1 of the SSTA and that of the SLPA, which shows that the SLPA is more nonlinear than the SSTA at the ENSO timescale. This is in contrast to the NLPCA study by Monahan [2001] and the NLCCA study by Hsieh [2001b], where the SSTA was found to be more nonlinear than the SLPA. The difference is that in the two earlier studies, time lags were not incorporated, so only non- linearity in space was being measured. In fact, NLPCA can be regarded as a limiting case of NLSSA, with the window L = 1 month. By incorporating time lags, NLSSA measures the nonlinearity in both space and time and reveals SLPA to be the more nonlinear field because of the stronger nonlinear interactions between ENSO and QBO timescales. [40] Perhaps the most intriguing achievement of the NLSSA is the isolation of a quasi-triennial oscillation of about 39 month period in both SSTA NLSSA mode 2 (Figures 9 and 10) and SLPA NLSSA mode 2 (Figure 14). In fact, the correlation between the SSTA and the SLPA in this 39 month oscillation is much higher than that found in the main ENSO signal (NLSSA mode 1) at the 4–5 year period. While a 3 year oscillation has been known to exist Figure 11. The NLSSA mode 1 SSTA when the eastern equatorial Pacific (Nino 3 region) is (a) warmest and (b) coolest. For comparison, the NLSSA mode 2 reconstructed SSTA when the Nino 3 region is (c) warmest and (d) coolest. The contour interval is 0.3C in Figures 11a and 11b and 0.15C in Figures 11c and 11d. Figure 12. The NLSSA mode 1 for the tropical Pacific SLPA plotted in the x1–x2–x3 3-D PC space and its 2-D planes. The straight lines indicate the linear SSA mode 1. 13 - 12 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS in the SO for a long time and, in fact, has been known to dominate over the 5 year oscillation during 1882–1926 [Troup, 1965], its close proximity to the main ENSO signal at the 4–5 year period renders a clean extraction difficult. NLSSA reveals this 39 month signal to be, in fact, considerably more regular, without the east-west shift of the SSTA between warm and cold events as found in NLSSA mode 1, and hence likely to be more linear than Figure 13. The reconstructed SLPA of NLSSA mode 1 along the equator from January 1971 to December 2000. The contour interval is 0.2 mbar. HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 13 the main signal (NLSSA mode 1). The 39 month period is also longer than most of the QBO periods reported by previous studies using linear techniques [Ghil et al., 2002]. Here the longer 39 month period appears to arise from nonlinearly combining the SSA modes of QBO periods and the longer-period SSA modes 1 and 2. This 39 month signal also appears less closely confined to the equator than the main signal. Figure 14. The NLSSA mode 2 reconstructed SLPA along the equator from January 1971 to December 2000. The contour interval is 0.1 mbar. 13 - 14 HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS [41] Acknowledgments. Support through strategic and research grants from the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. References Elsner, J. B., and A. A. Tsonis, Singular Spectrum Analysis, Plenum, New York, 1996. Ghil, M., et al., Advanced spectral methods for climatic time series, Rev. Geophys., 10.1029/2000RG000092, in press, 2002. Hamilton, K., Dynamics of the tropical middle atmosphere: A tutorial review, Atmos.-Ocean, 36, 319–354, 1998. Hsieh, W. W., Nonlinear canonical correlation analysis by neural networks, Neural Networks, 13, 1095–1105, 2000. Hsieh, W. W., Nonlinear principal component analysis by neural networks, Tellus, Ser. A, 53, 599–615, 2001a. Hsieh, W. W., Nonlinear canonical correlation analysis of the tropical Pacific climate variability using a neural network approach, J. Clim., 14, 2528–2539, 2001b. Kirby, M. J., and R. Miranda, Circular nodes in neural networks, Neural Comp., 8, 390–402, 1996. Kramer, M. A., Nonlinear principal component analysis using autoassocia- tive neural networks, AIChE J., 37, 233–243, 1991. Monahan, A. H., Nonlinear principal component analysis: Tropical Indo- Pacific sea surface temperature and sea level pressure, J. Clim., 14, 219–233, 2001. Rumelhart, D. E., G. E. Hinton, and R. J. Williams, Learning internal representations by error propagation, in Parallel Distributed Processing, vol. 1, edited by D. E. Rumelhart, J. L. McClelland, and P. R. Group, pp. 318–362, MIT Press, Cambridge, Mass., 1986. Smith, T. M., R. W. Reynolds, R. E. Livezey, and D. C. Stokes, Recon- struction of historical sea surface temperatures using empirical orthogonal functions, J. Clim., 9, 1403–1420, 1996. Troup, A. J., The ‘‘Southern Oscillation,’’ Q. J. R. Meteorol. Soc., 91, 490– 506, 1965. von Storch, H., and F. W. Zwiers, Statistical Analysis in Climate Research, Cambridge Univ. Press, New York, 1999. Woodruff, S. D., R. J. Slutz, R. L. Jenne, and P. M. Steurer, A compre- hensive ocean-atmosphere data set, Bull. Am. Meteorol. Soc., 68, 1239– 1250, 1987.  W. W. Hsieh and A. Wu, Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, B.C., Canada V6T 1Z4. (whsieh@eos.ubc.ca) HSIEH AND WU: NONLINEAR MULTICHANNEL SINGULAR SPECTRUM ANALYSIS 13 - 15