## Abstract

A new statistical time series prediction method based on temporal empirical orthogonal function (T-EOF) is introduced in this study. This method first applies singular spectrum analysis (SSA) to extract dominant T-EOFs from historical data. Then, the most recent data are projected onto an optimal subset of the T-EOFs to estimate the corresponding temporal principal components (T-PCs). Finally, a forecast is constructed from these T-EOFs and T-PCs. Results from forecast experiments on the El Niño sea surface temperature (SST) indices from 1993 to 2000 showed that this method consistently yielded better correlation skill than autoregressive models for a lead time longer than 6 months. Furthermore, the correlation skills of this method in predicting Niño-3 index remained above 0.5 for a lead time up to 36 months during this period. However, this method still encountered the “spring barrier” problem. Because the 1990s exhibited relatively weak spring barrier, these results indicate that the T-EOF based prediction method has certain extended forecasting capability in the period when the spring barrier is weak. They also suggest that the potential predictability of ENSO in a certain period may be longer than previously thought.

## 1. Introduction

The El Niño–Southern Oscillation (ENSO), manifested from a complex air–sea interaction in the equatorial Pacific, is the strongest climate signal in the interannual timescales (Rasmusson and Carpenter 1982; Enfield 1989; Philander 1990; Harrison and Larkin 1998). It has a quasiperiodic behavior with dominant periods of around 2–7 yr. Because of this quasiperiodicity, it is thought that ENSO may be predictable with a lead time of more than a year. However, most of the current prediction models suffer a marked drop in forecast skill when predictions are made across the boreal spring (the so-called spring barrier; Latif et al. 1994; Balmaseda et al. 1995). Therefore, these models generally exhibit useful skill only up to a lead time of 6–12 months (Barnston et al. 1994; Latif et al. 1998; Kim and North 1999). Nevertheless, the predictability of ENSO had been shown to vary on a decadal scale (Balmaseda et al. 1995). In decades when the seasonal dependence in the autocorrelation of the tropical Pacific sea surface temperature anomalies (SSTAs) was weak, such as in the 1980s, ENSO was more predictable. On the other hand, in decades when the autocorrelation was strongly phase-locked with the annual cycle, such as in the 1970s, ENSO was less predictable. Furthermore, Chen et al. (1995), using the physically based Cane–Zebiak model (Cane et al. 1986) with an improved initialization procedure, indicated that the spring barrier might be nonphysical and suggested that ENSO might be predictable up to 2 yr in advance. These results suggest that a well-designed statistical prediction model may also have the capability to predict ENSO beyond 1 yr in decades when the spring barrier is weak. In this study, I introduced a new T-EOF based prediction method that, by optimizing the use of historical data, showed extended forecasting capability for the Niño-3 SSTA index in such a situation (i.e., in the 1990s).

The fundamental concepts behind the development of this method are as follows. Let **s** = [*a*_{1}, … , *a*_{n}] be a time series of *n* samples and 𝗔 be a (*n* + *M* − 1)*M* matrix in which each of its columns consist of **s** plus (*M* − 1) empty elements. Either for the purpose of reconstructing the phase space trajectory (e.g., Broomhead and King 1986) or singular spectrum analysis (SSA; Vautard and Ghil 1989; Vautard et al. 1992), the 𝗔 entries are arranged as

Because of this kind of arrangement, one can further divide 𝗔 into three submatrices, 𝗔_{1}, 𝗔_{2}, and 𝗔_{3} as shown in Eq. (1). Among them, only 𝗔_{2} contains no empty elements. Therefore, 𝗔_{2} can be decomposed, by using the principal component analysis (PCA; Preisendorfer and Mobley 1988), as

where

are the *k*th row vector, the *j*th principal component (PC), and the *j*th empirical orthogonal function (EOF) of 𝗔_{2}, respectively. The superscript T in (2) denotes the matrix transposition. Because 𝗔_{2} is constructed solely from the time series **s**, its associated EOFs correspond to temporal patterns instead of spatial patterns as with a conventional PCA. This special PCA application to 𝗔_{2} is usually called SSA and the resultant EOFs and PCs are commonly referred to as temporal empirical orthogonal functions (T-EOFs) and temporal principal components (T-PCs). Thus, SSA is a PCA performed in the time domain with window length *M.* The resultant T-EOFs can also be seen as the phase space basis of a dynamic system reconstructed from **s**, while the resultant T-PCs describe the reconstructed phase space trajectory of the system.

However, the application of SSA to 𝗔_{2} alone does not extract all of the information embedded in **s**. For instance, the nonempty entries of 𝗔_{3} contain the most recent evolution history of **s**. Nevertheless, due to the existence of empty entries, SSA cannot be used to extract the information embedded in 𝗔_{3}. Because the most recent evolutionary history in **s** is likely to have the most influence on its own future behavior, the ability to use the information embedded in 𝗔_{3} may have a crucial effect on predicting its future evolution. To use this information, one observes that there are many entry overlaps between 𝗔_{2} and 𝗔_{3}. Since 𝗔_{2} includes all of the available data in **s**, the dominant T-EOFs extracted from it are likely to correspond to the dominant time series signals. Therefore, it is reasonable to expect that these T-EOFs are also the dominant T-EOFs of 𝗔_{3}. By projecting the 𝗔_{3} entries onto the T-EOFs from 𝗔_{2} to derive their corresponding T-PCs, one can construct a forecast value from these T-PCs and T-EOFs. Designed in this way, the T-EOF prediction method not only uses the whole range of data to extract the dominant T-EOFs, but also uses the most recent evolution history of the time series embedded in 𝗔_{3} as a predictor to produce forecasts. Because the most recent data are used in extracting the dominant T-EOFs, as well as in initializing the prediction, this method naturally put more weight on data that are more critical to the future evolution of the time series. Therefore, the T-EOF prediction method appears to have the property of automatically optimizing the information embedded in a time series to produce forecasts.

The detailed methodology of this T-EOF prediction method is described in section 2. The remainder of this paper is organized as follows. In section 3, the design of retroactive real-time forecast experiments for Niño-3, Niño-4, and Niño-3.4 indices are described. The results from these forecast experiments are shown in section 4. Finally, a summary and a conclusion of this study are given in section 5.

## 2. The T-EOF based prediction method

To produce the first forecast value **a**^{*}_{n+1}, from the above discussion, one assumes that the first row from 𝗔_{3} can be rather faithfully represented by the dominant T-EOFs from 𝗔_{2} with unknown T-PCs, that is,

where the subscript *L* is the number of T-EOFs retained from 𝗔_{2}. Because *a*^{*}_{n+1} and **y**_{1} are both unknown, one must use nonempty entries from the first row of 𝗔_{3} to solve **y**_{1} first. Therefore, by deleting the last row of **z**^{*}_{1} and 𝗕*, the equation for **y**_{1} becomes

Then, **y**_{1} can be solved as

Once **y**_{1} has been determined, one can simply multiply **y**_{1} by the last row of 𝗕* to yield

To find *a*^{*}_{n+r} for *r* ≥ 2, 𝗕 is assumed to remain the same as in (4), while the **z**_{r} entries consist of nonempty entries from the *r*th row of 𝗔_{3} and the previously predicted *a*^{*}_{n+1}, … , *a*^{*}_{n+r−1}. Then **y**_{r} is then obtained by substituting **z**_{r} for **z**_{1} in Eq. (5). Finally, *a*^{*}_{n+r} can be estimated as

From Eqs. (3)–(7) one can see that, in this T-EOF prediction method, 𝗕* is the forecast model with *a*^{*}_{n+r} and **z**_{r} predictands and predictors, respectively. It should be noted that there are two free parameters, *M* (window length) and *L* (number of retained T-EOFs), in this method. The *M* parameter determines the temporal resolution of SSA in decomposing a given time series, and the *L* parameter determines how many of the resultant T-EOFs are considered as dominant signals. The different choices of *M* and *L* means that the prediction method applies a different filter to the time series. Therefore, the method can have quite different performances with different values for these two parameters. One generally needs to scan through various *M* and *L* parameters to determine the values that optimize the performance of this prediction method.

## 3. Forecast experiment design

To evaluate the forecasting capability of this new prediction method, Niño-3.4 (5°N–5°S, 120°–170°W), Niño-3 (5°N–5°S, 90°–150°W), and Niño-4 (5°N–5°S, 160°E–150°W) SST indices between January 1950 and October 2000 were used in the retroactive real-time forecast experiments. These indices of area-averaged monthly SST are available from the National Centers for Environmental Prediction/Climate Prediction Center (http://www.cpc.ncep.noaa.gov/products/monitoring_and_data). Because there are two free parameters in this T-EOF prediction method, one can build many forecast models using different *M* and *L* choices. To properly conduct the experiments, each SST index was divided into three periods: (i) base, (ii) selection, and (iii) validation. The base period consisted of 30 yr of data from January 1950 to December 1979. Data in this period were treated as “initial historical data” and were available for use in all forecast models. The selection period consisted of 9 yr of data from November 1983 to October 1992. Data in this period were used to evaluate the performances of forecast models using various *M* and *L* parameters to select the optimal forecast model for the time series. The validation period consisted of 8 yr of data from November 1992 to October 2000. Data in this period were used to validate the predicting skill of the selected optimal model.

The forecast experiments were conducted as follows. First, three annual-cycle-removed SSTA time series were constructed by subtracting the mean annual cycle of the base period from whole records of each Niño SST index. Then, a total of 800 forecast models, each corresponding to a specific *M* and *L* for *M* from 30 to 225 months, with an increment of 5 months and for *L* from 11 to 30 modes, were set up for each time series to prepare for retroactive real-time forecasts from January 1980 until the end of the selection period. During this period, a set of forecasts with lead times of 1–48 months was produced for each month, based on data prior to that month only, with every one of these models. Each model's forecast capability was evaluated based on the temporal correlation coefficients between its forecasts and the observations. Because forecasts at different leads corresponded to slightly different periods, to appropriately compare the forecasting skill among different leads, the respective temporal correlation coefficients were calculated from data only in the common period (i.e., the selection period). The optimal forecast model for each time series was then chosen by selecting the model with the highest averaged temporal correlation coefficient among lead times of 0–36 months. The optimal forecast model was used to produce forecasts to the end of the validation period.

## 4. Results of retroactive real-time forecast experiments

To select the optimal forecast model for each of Niño SSTA index, Fig. 1 shows the averaged temporal correlation coefficients among 0–36 months lead time for all 800 models in the selection period for (Fig. 1a) Niño-3.4, (Fig. 1b) Niño-3, and (Fig. 1c) Niño-4 as a function of window length. From this figure and the figure with the averaged temporal correlation coefficients as a function of retained modes (not shown), one can see that models with *M* = 195, 190, 190, and *L* = 23, 25, 15 yielded the highest averaged temporal correlation coefficients for Niño-3.4, Niño-3, and Niño-4 SSTA indices, respectively. Hence, they were chosen as the optimal models. The need of nearly 16 yr of window length for the optimal models may result from the quasi-quadrennial oscillation of ENSO that requires several cycles of evolution to properly describe its behavior.

To evaluate the predicting capability of these optimal models, Figs. 2–4 show the temporal correlation coefficients and the root mean square errors (rmse's) between the observations and the forecasts (plus sign) as a function of lead time in both selection and validation periods for Niño-3.4, Niño-3, and Niño-4 SSTA indices, respectively. In each panel of these figures, results from three reference models are also shown. The first one is the persistence forecast (thin line). The second one is the forecast from the autoregressive model of order 17 (AR-17; dotted line) that was also used as a reference model in Kim and North (1999). The third one is the forecast from the optimal AR model (i.e., the AR model of order that yielded the best correlation skill in the selection period; thick line). The optimal AR models for Niño-3.4, Niño-3, and Niño-4 SSTA indices are AR-46, AR-46, and AR-44, respectively. The orders of these models are consistent with the dominant quasi-quadrennial behavior of ENSO in the 1980s.

From the figures, one noted that the optimal AR models generally had higher correlation coefficients than others in the selection period. However, their skills deteriorated severely in the validation period. On the other hand, although the optimal T-EOF based models did not yield correlation skills as high as those of the optimal AR models in the selection period, they did not lose as much skill in the validation period. The optimal T-EOF based model for the Niño-3 SSTA index, as shown in Fig. 3, even had better skill in the validation period than in the selection period. Generally speaking, the optimal T-EOF based models had worse correlation skills than other models for a lead time less than 6 months in both periods, but had better correlation skills for a lead time longer than 6 months in the validation period. As for the prediction errors, rmse's of the optimal T-EOF based models were generally larger than those of the AR models in the selection period, but were more comparable in the validation period. The only noticeable exception was the prediction of Niño-3 SSTA index in the validation period, where the optimal T-EOF based model had lower rmse's than those of AR models for a lead time longer than 16 months. These results suggest that the performance of the optimal T-EOF based models is more consistent than AR models.

The difference in performances between the optimal AR model and the optimal T-EOF based model in different periods can be understood by comparing the basic designs between these two methods. One notes that an AR model of order *M*-1 and a T-EOF based model of window length *M* use the same historical data to build the models and use the same initial conditions to produce predictions. However, the AR method constructs the model by finding coefficients that optimally fit the historical data. Because of this optimal fitting strategy, both signal and noise components of the historical data contributed to determining the coefficients of the optimal AR models in the selection period. Hence, the high skills of these models shown in the selection period were unable to carry over to the validation period. On the other hand, the T-EOF prediction method has a built-in filter process (i.e., the method has the freedom to choose the number of retained modes *L*) to reduce the risk of bringing too much of the contribution from the noise component of the historical data into the model. The above results indicated that this built-in filter process caused the optimal T-EOF based models to have higher rmse's and lower correlation skills than the optimal AR models in the selection period. However, it also allowed the performance of these models to be more consistent in both periods.

From the figures, one also noted a marked change in behavior from persistence curve (which also corresponds to lag autocorrelations) comparisons between the selection and the validation periods for Niño-3.4 and Niño-4 SSTA indices. A strong return in positive correlation occurred near a lead time of 4 yr in the selection period, while the similar feature occurred near a lead time of 3 yr in the validation period. It is believed that hardly any statistical model can have the capability to cope with this kind of marked change in behavior. Therefore, the deterioration of the forecasting skill shown by both the AR and the T-EOF methods in the validation period for Niño-3.4 and Niño-4 SSTA indices might not be unexpected. The situation for the Niño-3 SSTA index was somewhat different. The persistence curve for the Niño-3 SSTA index in the validation period showed only a weak return in positive correlation near a lead time of 3 yr compared to that for Niño-3.4 or Niño-4. Because the change in lag autocorrelations was less dramatic, the performances of the optimal AR and T-EOF based models for the Niño-3 SSTA index were more consistent in both periods. These results suggest that the steadiness of the lag autocorrelations of a time series may have a very important effect on the performance of any statistical prediction model.

The above results indicate that the T-EOF prediction method indeed has certain extended forecasting capability for Niño SSTA indices, especially the Niño-3 index, in the 1990s. To examine the closeness of individual forecast to the observations, all available forecasts for the Niño-3 SSTA index starting from the beginning of validation period for each lead time were plotted as a function of time in Fig. 5. Generally speaking, the evolutions of short-lead forecasts more closely followed the evolution of the observations, while those of long-lead forecasts tended to be smoother than the observations. However, none deviated much from the observations, except in the period of 1995–97 where forecasts of almost all leads were out of phase with the observations. The forecast failure in this period was not a unique property of this prediction method. It had been shown (Knafe and Landsea 1997) that nearly all forecast models—statistical and dynamic—also failed in this period. The rare occurrence of a prolonged ENSO warm event and the less coherent behavior between the SST and atmospheric conditions in the first half of the 1990s was shown to be responsible for the lower forecasting skill of every existing prediction model during this period (Ghil and Jiang 1998).

To examine if the T-EOF prediction method also suffered the spring barrier problem, Fig. 6 shows the monthly stratified correlation coefficients for the Niño-3 SSTA index in the validation period as a function of the starting months. The predominant feature of this figure is the quasi-annual oscillatory behavior of these monthly stratified correlation coefficients. The minima generally occurred when predictions were initiated from the months of March or April, while the maxima generally occurred when predictions were initiated from the months near the end of the year. Therefore, the T-EOF prediction method still encountered the spring barrier in the validation period. However, the minima of these correlation coefficients never dropped below 0.3 and the maxima showed only slightly interannual variation. It appears that the spring barrier is relatively weak in the 1990s.

Is the extended forecasting capability of this T-EOF prediction method dependent on the strength of the spring barrier? The same forecast experiments with a reversed time series for each Niño SSTA index were used as cross validations. Results (not shown) indicated that none of the T-EOF prediction models or the AR models could yield useful extended forecast skill. The autocorrelation of the tropical Pacific SSTA was more strongly phase-locked with the annual cycle in the 1950s and 1960s than in the 1980s and 1990s (see Fig. 12 of Latif et al. 1998), which suggested that these prediction models encountered a stronger spring barrier in the cross-validation period. Therefore, it seems that the predicting skill of the T-EOF based method does strongly depend on the strength of the spring barrier.

## 5. Summary and conclusions

In this study, a new T-EOF based statistical prediction method was introduced. Results from the above forecast experiments showed that this method had certain extended forecasting capability for ENSO in the period when the “spring barrier” was weak. Compared to the AR method, one can easily identify that the extended forecast capability of this T-EOF based prediction method comes mainly from the use of an optimal subset of T-EOFs to prevent forecasts from being contaminated by the noise component of the historical data. However, this built-in filter process also appears to filter out too many high-frequency signals from the data to deteriorate the correlation and rmse skills of this method at short leads. Thus, the T-EOF based method is probably more suitable to be used in predicting time series with predominant low-frequency variability. On the other hand, because this method is based solely on data manipulations and is explicitly independent of any physical rules, it can potentially be applied to any time series to produce valuable forecasts. Furthermore, the easy extension of SSA to multichannel SSA (Plaut and Vautard 1994; Zhang et al. 1998) suggests that it can be easily extended to a spatial–temporal empirical orthogonal function based prediction method. In summary, this study demonstrated that a wise use of historical data could greatly improve our predicting capability. Furthermore, it indicates that the potential predictability of ENSO may be greater than previously thought.

## Acknowledgments

I thank anonymous reviewers and the editor for their critical review of this paper. Their comments were essential in improving the final version of this paper. This study was supported by the National Science Council of Taiwan through Grants NSC 89-2111-M-008-005 and NSC89-2111-M-008-061.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* Dr. Yung-An Lee, Institute of Atmospheric Physics, National Central University, Chung-Li 32054, Taiwan. Email: yalee@atm.ncu.edu.tw