The theoretical predictability limit of El Niño–Southern Oscillation has been shown to be on the order of years, but long-lead predictions of El Niño (EN) and La Niña (LN) are still lacking. State-of-the-art forecasting schemes traditionally do not predict beyond the spring barrier. Recent efforts have been dedicated to the improvement of dynamical models, while statistical schemes still need to take full advantage of the availability of ocean subsurface variables, provided regularly for the last few decades as a result of the Tropical Ocean–Global Atmosphere Program (TOGA). Here we use a number of predictor variables, including temperature at different depths and regions of the equatorial ocean, in a flexible statistical dynamic components model to make skillful long-lead retrospective predictions (hindcasts) of the Niño-3.4 index in the period 1970–2016. The model hindcasts the major EN episodes up to 2.5 years in advance, including the recent extreme 2015/16 EN. The analysis demonstrates that events are predicted more accurately after the completion of the observational array in the tropical Pacific in 1994, as a result of the improved data quality and coverage achieved by TOGA. Therefore, there is potential to issue long-lead predictions of this climatic phenomenon at a low computational cost.
Skillful long-range forecasts of El Niño–Southern Oscillation (ENSO) are still in high demand. After decades of extensive efforts, dynamical models nowadays represent the best available tools to issue ENSO forecasts at lead times of up to two seasons, although they are still largely constrained by the lack of complete understanding of the physics of the phenomenon, by problems arising from the initialization of the components of the climate system, or by the need for accurate parameterization of important physical processes (Barnston et al. 2012). Statistical models, on the other hand, largely depend on the availability of ocean and atmosphere historical data, so that the longer the length of the data, the more robust is the predictor–predictand relationship identified by the model (Barnston et al. 2012). In addition to these factors, the low signal-to-noise ratio in boreal spring (Sarachik and Cane 2010), the influence of high-frequency atmospheric winds (Fedorov et al. 2003, 2015), and the natural irregularity of the climate system (Wittenberg 2009) all limit the long-term dynamical and statistical forecasting of the phenomenon. Some of the classical ENSO theories view the oscillation as self-sustained (Cane et al. 1990; Jin et al. 1994; Jin 1997) and support the claim that it is potentially predictable several years in advance (Cane et al. 1986; Goswami and Shukla 1991; Latif et al. 1998; Chen and Cane 2008; Wittenberg et al. 2014; Gonzalez and Goddard 2016; Luo et al. 2016; DiNezio et al. 2017; Astudillo et al. 2017), but only a handful of studies document such long-lead retrospective forecasts of past events (Latif et al. 1998; Chen et al. 2004; Luo et al. 2008; Izumo et al. 2010; Ludescher et al. 2013, 2014; Petrova et al. 2017; Gonzalez and Goddard 2016; Ramesh et al. 2017; Luo et al. 2017) and most of them use dynamical models. Statistical models are assumed to be less skillful at long lead times, and comparable in performance to dynamical schemes at shorter lead times of about half a year (Barnston 1994; Chen and Cane 2008). To some extent this is explained by the fact that a new generation of statistical models has not been added to the ENSO forecasting plume, and the majority of the old models have not been substantially revised in the recent years, and some since they were created in the 1980s and early 1990s (Barnston et al. 2012).
One of the strongest events on record—the 1982/83 El Niño (EN)—surprised the scientific community (Cane et al. 1986; McPhaden and Yu 1999) as it was neither predicted nor identified until very late in its development. This triggered a decade-long effort to put in place a monitoring system in the tropical Pacific with the aim of studying ENSO better and improving the predictive capacity of models (McPhaden and Yu 1999), which led to the inauguration of the TOGA research program in 1985 (McPhaden and Yu 1999). It deployed a three-dimensional array in the tropical Pacific that since then regularly samples the subsurface temperature down to 500-m depth. The number of monthly temperature profiles increased dramatically (see Fig. S1 in the online supplemental material). The system was completed in 1994, just in time to track the stronger-than-normal trade winds in 1995/96, which generated a buildup of warm waters in the western tropical Pacific more than one year before the peak of the record-breaking 1997/98 EN (McPhaden and Yu 1999). This was the first time when the scientific community and the public could see the benefits of TOGA. A number of studies now fully recognize the fundamental role that the intensification of the trade winds and the subsurface heat buildup in the western equatorial Pacific play in the onset of EN events (Wyrtki 1985; Cane et al. 1986; Jin 1997; Clarke and Van Gorder 2003; McPhaden 2003, 2004; McPhaden et al. 2006; Ramesh and Murtugudde 2013; Ballester et al. 2015; Petrova et al. 2017), and statistical models can benefit from available data to represent in more detail these processes that occur early on in the generation of the events.
In the present study we use an improved version of the flexible statistical dynamic components ENSO model described in Petrova et al. (2017). At long lead times it incorporates predictor variables designed to capture the three-dimensional shape of the warm pool subsurface heat buildup at different depths, as well as zonal wind stress anomalies in the central and western equatorial Pacific (see section 2). The aim is to capture the low-frequency deterministic and state-dependent portions of the variability and coupling between the ocean and atmosphere (Eisenman et al. 2005; Gebbie and Tziperman 2009; Hu and Fedorov 2016; Levine and Jin 2017), from which predictability can be derived (Latif et al. 1998; Chen et al. 2015). The model consists of several stochastic cycle components with frequencies corresponding to the main peaks in the spectrum of the Niño-3.4 index (see Petrova et al. 2017), as well as predictor regression variables such as sea surface and subsurface temperature and zonal wind stress. These variables enter the model equations in the form of lagged time series with respect to the monthly value of the Niño-3.4 index, and are selected to be consistent with the EN dynamical evolution. In this way, different covariates are used for predictions at different lead times, depending on the average temporal progression of EN events. In the present study we show hindcasts from the model, and as is normally the case, a somewhat poorer performance is expected in operational mode. In fact, the model has been operational since 2015, and while it correctly detected the 2015 EN and the mild La Niña (LN) in 2016, it failed to foresee the recent LN in 2017 and predicted neutral conditions instead (not shown). This paper is organized as follows: in section 2 we describe the data and methods used in the analysis; in section 3 we present the results and then discuss them in section 4, where we also provide concluding remarks.
2. Data and methods
The model used in this study is an advanced version of the statistical dynamic components model proposed by Petrova et al. (2017) and developed specifically for prediction of the average sea surface temperature in the Niño-3.4 region defined as the box 5°N–5°S, 170°–120°W. It is a statistical model that belongs to the class of dynamic components time series models. The distinctive feature of this type of models is that they decompose the time series of interest into dynamic components that represent linear stochastic processes with separate evolutions (Durbin and Koopman 2012). The addition of predictor variables, in this case derived from lead–lag climate composites, is done using regression. We refer to the appendix for complete and mathematically precise details.
The model first presented in Petrova et al. (2017) is built in terms of two main subgroups of elements. The first subgroup contains the so-called dynamic components, which include a trend (level), a seasonal, and three time-varying cyclical (quasi-periodic) components. The second subgroup contains a number of individually selected predictor variables, which enter the model equation in the form of regressed and lagged time series and will be described later. All of these separate components are added together in a linear fashion to form the final ENSO model given by
where yt represents the average monthly temperature in the Niño-3.4 region at time t, μt is the trend component, γt is the seasonal component with 12 seasonal effects (one fixed value for every month of the year), and ψ1t, ψ2t, and ψ3t are the stochastic cycle components. Also, Xtβ is a vector that contains the predictor variables, while εt is the noise term in the model.
Here we improve this first version of the model, by replacing the previously fixed seasonal component with two slowly varying annual and semiannual periodic components, and also by including one additional time-varying cycle component, so that the new model equation becomes
Thorough information about the different components and how they are modeled and estimated is provided in the appendix.
The previous model presented in Petrova et al. (2017) used three quasi-periodic cycle components that generally correspond to the near-annual (NA), quasi-biannual (QB), and quasi-quadrennial (QQ) modes of ENSO variability, while here we added one more stochastic cycle component associated with ENSO variability on decadal (D) time scales. In Petrova et al. (2017) we established that this low-frequency variability is important for the simulation of some EN events, and this feature was not explicitly resolved in the previous model version. We have also replaced the fixed seasonal component in the previous version of the model with two new cyclical components bound to annual (~12 month) and semiannual (~6 month) periodicities. They are allowed to vary slowly over time in order to address the finding in our previous study that the annual frequency of the seasonal component was not sufficiently well simulated, because the annual periodicity of the Niño-3.4 temperature is not strictly fixed at 12 months (Chen et al. 2016), and especially because during EN events the amplitude of the annual cycle is suppressed (Guilyardi 2006). As a result, we have a total of six stochastic cycle components in the new model version.
There are also different regression predictors included in the model at different lead times, all selected based on the general evolution of an average EN event. LN is assumed to be symmetrical, although we are aware that important asymmetries exist between the two and this problem will be addressed in future work. In the ocean we used both surface and subsurface temperatures at different depths (between 0 and 500 m) and regions for the extraction of the predictors. Regions are selected in the western and central equatorial Pacific where the ocean is typically warmed abnormally prior to EN and a heat buildup occurs during the growing and recharge phase (Ballester et al. 2016a; Petrova et al. 2017). Figures S2 and S3, as well as Table 1, show the selected regions and depths considered at different lead times. The selection is based on climate composites of EN events from the period 1978–2012 (also see Petrova et al. 2017). The sea surface temperature datasets used for the predictors and for the Niño-3.4 temperature time series are NOAA ERSST-V3 (https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.ersst.html) before 1982 and NOAA OISST V2 thereafter (https://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html). The subsurface temperature dataset used for the subsurface ocean predictors is the Subsurface Temperature and Salinity Analyses dataset by Ishii et al. (2005) archived at https://rda.ucar.edu/datasets/ds285.3/ before 2012 and the Hadley Centre EN4.0.2 analyses data thereafter (Good et al. 2013). In the atmosphere three different regions are used to extract zonal wind stress predictors for the model. The three regions are again located in the western and central equatorial Pacific (see Fig. S4 and Table 1) and the dataset is the NCEP–NCAR reanalysis (Kalnay et al. 1996).
During forecasting, the dynamic components (especially the stationary cycles) have larger weights for the mid- and short-term forecasts, while the impact of the predictors remains the same for short- and long-term forecasts. Hence, the predictors become relatively more important for long-term forecasting (also see the appendix for more information). Importantly, the predictor variables also affect the estimation of the cycle components parameters for each forecast. Parameter estimation relies on the Kalman filter methods (Kalman 1960; Harvey 1989) and on state space methods (Durbin and Koopman 2012).
Results for correlations and root-mean-square errors are obtained as follows: the Niño-3.4 predictions in the period 1972–93 are based on parameter estimates (calibration process) from data in the period 1952–70, while the predictions in the period 1994–2015 are based on parameter estimates from data in the period 1974–92. In this way, to avoid the heavy computations, we have produced the predictions on which correlations and root-mean-square errors are based using a prefixed period for calibration purposes. In comparison, all other predictions (including parameter estimation) are based on the observations available before each starting prediction point. Still, data for the prediction estimations were progressively excluded, in order to include only the more recent samples and discard earlier data of presumed lesser quality. Thus, predictions up to 1990 were made using the data from 1952 onward, predictions between 1991 and 1996 were made using the data from 1970 onward, and predictions thereafter were made using the data from 1982 onward.
A limitation of the study is that the performed predictions are not operational, as they are based on retrospective hindcasting experiments. Our system also strongly relies on the model variability skeleton, contributed by, among others, different cyclical components. However, all ENSO forecasting systems, including operational dynamical models, implicitly or explicitly rely on intrinsic ENSO variability generated at the cyclical low-frequency modes for prediction (Kirtman and Schopf 1998). As an example, we include the spectrum of Niño-3.4 from a long (500 yr) spinup simulation with the GFDL CM2.1, which is one of the operational models for ENSO prediction, in order to compare it with the spectrum of both the Niño-3.4 observations and their predictions with the model proposed here (Fig. S5). What can be clearly noticed is that the power density is distributed similarly in all cases, with main peaks corresponding to the NA, QB, QQ, and D modes of variability, respectively, also used as cyclical time-varying components in our model.
The observed and hindcast monthly Niño-3.4 anomalies at 6 and 24 months lead time are presented in Fig. 1. The 6-month lead hindcast predicts the timing and magnitude of all EN and LN events, and no false alarms are generated [root-mean-square error (RMSE) = 0.54; Fig. 1a]. Since an ENSO event is typically already underway half a year before its peak in December–February (DJF), the majority of the operational forecasting schemes are able to produce accurate predictions at this lead time (Barnston et al. 2012). The 24-month lead hindcast, and in general any lead time hindcast beyond the spring barrier (i.e., from 8 months onward; not shown), generally reproduces the crests and troughs in the time series (RMSE = 0.62; Fig. 1b). However, for the period before the prominent 1997/98 EN, we find that the predicted amplitudes of the larger events are notably smaller than the observed and sometimes an event is hardly or not detected. We highlight that this cannot be explained by a change in the interannual ENSO activity in the different time periods, as three sizeable EN (1972/73, 1982/83, 1986/87; also 1997/98, 2009/10, 2015/16) and LN (1973/74, 1975/76, 1987/88; also 1998–2000, 2007/08, 2010/11) episodes have occurred before and after 1994 (CPC 2016). In addition, it cannot be simply attributed to the design of the model and the predictor variables used, because EN events from both periods were considered for the composites on which the selection of predictor variables was based [see Petrova et al. (2017) for details].
To characterize better the difference between periods, Fig. 2 displays the regressions between the observations and hindcasts for two consecutive 22-yr subperiods (1972–93 in blue and 1994–2015 in red) at 6- and 24-month leads. No substantial difference is observed between the slopes of the regression lines for the two periods at the shorter lead time (regr1972–93 = 0.65, t = 23.88, regr1994–2015 = 0.74, t = 27.34, p < 0.001; Fig. 2a), indicating that the model performance is comparable. Conversely, the regression coefficients significantly increase for the long-range hindcasts made after 1994 (regr1972–93 = 0.35, t = 17.12, regr1994–2015 = 0.65, t = 30.93, p < 0.001; Fig. 2b), which represents a major improvement in the capacity of the model. The change in the overall similarity between the observations and the hindcasts at 24-month lead time is also assessed by the 16-yr moving RMSE shown in Fig. 1c. The RMSE decreases monotonically with time until the early 1990s and then stays relatively constant afterward. At the same time, data availability was constantly improving during TOGA, until the tropical Pacific network array of moorings was fully into place at the end of the program in 1994 (McPhaden et al. 1998).
To further explore the difference in the model performance over the two periods, Fig. 3 shows correlations and root-mean-square errors for the whole range of lead times up to 24 months. For lead times of about two seasons both the correlations and RMSE are similar among periods, while for lead times beyond 6 months they start to diverge. We also observe that correlations and RMSE stay relatively constant beyond this lead time. One possible reason for this stable behavior is that the stochastic quasi-periodic cycles are the main contributors to the skill of the predictions and their unknown parameters are estimated similarly by the Kalman filter at different lead times beyond two seasons. Generally, the skill is derived from information about subsurface heat anomalies of approximately the same intensity, and the different cycles capture similar oscillation phases. Previous studies (Chen et al. 1995, 2004; Chen and Cane 2008; Stockdale et al. 2011; Duan et al. 2016; Lee et al. 2018) have already concluded that the spring predictability barrier may not be an intrinsic barrier to the system itself, but it could rather depend on model skill, observational data availability (especially in the subsurface western tropical Pacific; Lee et al. 2018), and the precursors used. Warm water volume (WWV) in the tropical Pacific as a predictor (i.e., subsurface information) is not associated with a spring persistence barrier, and its correlation with the Niño-3.4 is above 0.7 for the February–April season when SST anomalies have the lowest correlation (McPhaden 2003). Here we add evidence to such claims, as we also find that the drop in forecast skill is slow and gradual for longer-lead predictions than a couple of seasons (Fig. 3).
The statistical model we use is linear, and while its stochastic cyclical components are mainly responsible for capturing the correct phase of the oscillation, the lagged predictor variables are expected to contribute to the correct forecasts of the amplitudes of the events, especially at longer lead times [see section 2 and the appendix herein, and Petrova et al. (2017) for details]. Below we analyze if the predictor variables add significantly to the EN hindcasts of the earlier period, which also coincides with a time when no regular subsurface temperature and wind stress data were being provided yet (McPhaden et al. 1998).
The hindcasts at several lead times of the strongest EN events in the study period (see CPC 2016) are displayed in Fig. 4. In all cases the model is capable of detecting a warming 29 months in advance (magenta curve), although there are evident errors in the amplitude and timing in some cases. A much better representation of the amplitudes in the long-lead hindcasts of the events in the second period (1997/98, 2009/10, and 2015/16), as compared to those occurring in the first period (1972/73, 1982/83, and 1986/87), is also clearly visible in the figure. The estimated coefficients and the corresponding t and p values for the predictor variables used in the 24-month lead predictions of all the warm events in the study period are listed in Table 2. Remarkably, none of the three predictor variables is found to be significant at the 90% level for the hindcasts of any of the events before 1994, while there is at least one significant variable for each hindcast of the episodes that occurred afterward. Similar results hold for the other long-lead predictions shown in Fig. 4 (Tables S1 and S2).
4. Discussion and conclusions
We demonstrated that the Tropical Pacific Observing System, and especially the provision of subsurface temperature data on a regular basis, has a vital contributing role (Newman et al. 2011) for the long-lead predictive capabilities of the model proposed here. With the end of TOGA in 1994 nearly the whole equatorial band between 10°N and 10°S was covered with moorings (McPhaden et al. 1998), and this is also the start of altimetry data (Stockdale et al. 2011). As can be seen from the Tropical Atmosphere Ocean–Triangle Trans-Ocean Buoy Network (TAO-TRITON) array development (NOAA 2018a), some subsurface data from the central Pacific were already streamed at the end of 1987, and at the end of 1991 data were also coming in from the western Pacific, which represents a key region for the forecast of the phenomenon at lead times beyond the spring barrier. Thus, almost three decades have passed since the three-dimensional observations began in the tropical Pacific, and the limited span of the data is now less of a problem for the robust definition of statistical predictive schemes (Barnston et al. 2012).
As seen in the previous section, there is a well-defined shift between the lack of significance of the predictor variables for the hindcasts of the warm events before the end of TOGA and their significance thereafter. Our results strongly support the view that the improved hindcasts are due to the availability of regular and higher-resolution subsurface data ensured by the implementation of the observational network array (NOAA 2018b). This is also confirmed by Fig. S6, which shows the same hindcasts as in Fig. 4, but made without the inclusion of the predictor variables in the model framework. The lack of predictors in the model results into a clear deterioration of the hindcasts of the EN events from the period after 1994, but in no substantial difference in the hindcasts of the events from the earlier period (also see Table S3).
The correct and relevant subsurface information also has implications for the forecasting of the magnitudes of the warm events (Ballester et al. 2016b, 2017). In the linear framework of the model that we use, at the longer lead times the predictor variables have more forecast weight than they do at the shorter lead times (see section 2 and the appendix). The predicted amplitudes of the three earlier events shown in Figs. 4a–c do not exceed 1.5°C at the long lead times of 21 and 29 months (green and magenta curves). At the same time, the predicted amplitudes of the three events that took place in the later period, when the predictor variables are shown to have an impact (Table 2), are consistent with the occurrence of a strong EN event (green and magenta curves in Figs. 4d–f). Some of the underestimation of the amplitudes of the events predicted at long lead times is also due to the stochastic noise component of the zonal wind (Penland 1996; Hu and Fedorov 2016; Levine and McPhaden 2016) as more extreme EN events have been found to result from more intense and frequent westerly wind bursts (Chen et al. 2015). Additionally, in the case of the 1982/83 EN anomalous solar radiation and suppressed convection may have played a more decisive role, setting this event apart from the others (Kim and An 2018) and making the predictors used here at long lead times less relevant (Kirtman and Zebiak 1997).
In essence, the same conclusion as the one reached here was also made by Stockdale et al. (2011), where a large reduction of the errors in Niño-3.4 SST forecasts made after 1994 is detected with the European Centre for Medium-Range Weather Forecasts (ECMWF) Seasonal Forecast System 3. The results are also in agreement with an earlier study with the same system (Balmaseda and Anderson 2009), in which the effect of Argo floats is removed from the observations, and it is established that improvements in the forecast are clearly explained by the improved observing system. Further, it was found that the information from the mooring array is the main contributor for the increased skill of the prediction system in the equatorial Pacific region. In addition, McPhaden et al. (2006) using an empirical ENSO model with two predictors—WWV in the equatorial Pacific and an index of the Madden–Julian oscillation—documents much better estimations of the Niño-3.4 after 1995, with lower than observed amplitudes before this year, just as in the results presented here. Similarly, the authors attribute the improvement to the better observations after the placement of the TAO array.
Conversely, a more recent study (Kumar et al. 2015) concluded that the increase of the number of observations after 1994 did not result in a clear improvement of the prediction skill of the National Centers for Environmental Prediction (NCEP) System 2. We note, however, that our results are not directly comparable, because the forecasts discussed therein are performed at up to 6 months lead time, when essentially an SST anomaly signature of a developing EN or LN event is already present in the eastern equatorial Pacific, and subsurface information is generally not as crucial as it is at the longer lead times discussed here. The authors themselves admit that the evolution of the ocean–atmosphere system at this short lead is affected much more by the surface wind and ocean circulation feedbacks. SST is in fact found to be a more useful predictor for forecasts started 6 months before the event than WWV (McPhaden 2003).
Some of the existing statistical systems already include measures of integrated equatorial heat content (Barnston et al. 2012). However, our model uses temperature data from a selection of dynamically relevant regions and depths to maximize its predictive power. These values may not always be well represented by spatially integrated measures of heat content, and our analysis suggests that the integration sometimes masks the intensity of the heat buildup in specific regions in the subsurface at long lead times, and more importantly, does not allow the systems to properly track the eastward propagation of heat along the equatorial thermocline (Ballester et al. 2015; Petrova et al. 2017). WWV anomalies along the whole equatorial Pacific present in late boreal winter and spring (February–May) are persistent until next boreal winter, but those in early boreal winter are not. Hence, as a predictor it could extend the lead time to about a year in advance, but not more (Izumo et al. 2019). Alternatively, WWV calculated only in the western equatorial Pacific is significantly correlated with Niño-3 SST anomalies for much longer lead times of more than 20 months (McPhaden 2003), and is a significantly better predictor beyond the spring barrier (Izumo et al. 2019). Sea surface height (SSH), on the other hand, is also not always representative of the heat accumulation in the warm pool, because sometimes positive and negative heat anomalies exist at different depths of the water column near the thermocline, and the net result is a lack of a prominent SSH anomaly (see Figs. S7 and S8). The combination of the memory effect represented by subsurface information, weakly varying seasonality and nonlinearity, on the other hand, could be sufficient for reproducing the overall ENSO variability (Chen et al. 2016), and our model design attempts to incorporate these particular effects.
Although there is a marked difference in the predictive capacity of the model during the earlier and later subperiods, it still exhibits high skill (i.e., correlations and RMSE in Fig. 3) in both periods. We conclude that statistical models should be improved in the direction of using the available subsurface information that is fundamental for ENSO in a more discrete and targeted way, so that they can provide early and useful information about EN and LN events to decision-makers around the world, which could prevent threats to human lives and reduce economic costs.
D.P. acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under grant agreements 727852 (project Blue-Action). J.B. gratefully acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under grant agreements 727852 (project Blue-Action), 730004 (project PUCS), and 737480 (Marie Sklodowska-Curie fellowship ACCLIM). X.R. gratefully acknowledges funding from the Daniel Bravo Foundation (project WINDBIOME), as well as from the PERIS PICAT project (SLT002/16/466) and the NEW INDIGO project (PCIN-2013-038).
The most basic version of the class of dynamic components time series models is the local level model for a univariate time series yt and is given by yt = μt + εt, where μt is a linear stochastic process, dynamically evolving over time, and εt is a noise term. We can consider μt to follow a random walk process that captures the long-term trend features in the time series and εt to be an independent and identically distributed (IID) variable that represents the short-term deflections from the trend or the noise in the time series. The trend signal μt is the key feature of interest in the local level model, also for generating long- or medium-term forecasts of yt. Its random walk process is given by μt+1 = μt + ηt where the noise term ηt is IID. Under the assumption that both noise terms εt and ηt are normally distributed with mean zero and variances and , respectively, the celebrated Kalman filter equations (Kalman 1960) compute the minimum mean square error (MMSE) estimates of μt given realizations for y1, y2, …, yt, in a recursive real-time fashion, for t = 1, 2, …, T, where T is the length of the time series under investigation. The estimate of μt can be expressed as the weighted average where the weights are normalized (they sum up to unity; w0 + w1 + w2 + … = 1), are exponentially decaying, and are a function of the signal-to-noise ratio . When q is relatively large ( is small relative to , implying that yt behaves close to a random walk process as yt ≈ μt), the weights are decaying fast to zero and we obtain a “noisy” estimate of μt. This estimate may be representative as it is close to the local level (small estimation bias), but given that only a few observations are used for the estimation, the precision is typically small (i.e., large estimation variance). When q is relatively small ( is small relative to , implying that μt is evolving slowly over time as μt+1 ≈ μt), the weights are decaying slowly to zero and we obtain a smooth estimate of μt. In the latter case, the estimation bias may be larger (less local targeting), but the estimation variance is smaller since more observations are used for estimation. The appropriate value for the signal-to-noise ratio q for a particular time series depends on the dynamic features of the time series. We estimate q by the method of maximum likelihood, which entails the numerical maximization of the likelihood function that is computed using the Kalman filter (for a specific value of q). The h-step ahead forecasting (that is the estimation of yT+h, given realizations for y1, y2, …, yT) is also computed by the forward-moving Kalman filter (from 1 to T). The estimation methodology provides the MMSE optimal weights for forecasting: the forecasting weights gradually decline when observations are increasingly remote from the forecast point as these become increasingly less relevant. The estimation of all μt values given the realizations y1, y2, …, yT (all data) is referred to as signal extraction and relies on Kalman smoothing, which is a backward-moving filter (from T to 1); see Durbin and Koopman (2012, ch. 2) with all methods for filtering, forecasting, signal extraction, and parameter estimation, and with related details for the local level model.
The local level model is a special case of the dynamic components model adopted in Petrova et al. (2017) where the observation equation yt = μt +εt is extended with regression effects (predictors) and more linear stochastic processes that represent key dynamic features of the Niño-3.4 temperature time series including seasonal and cyclical effects. The model then becomes , where Xt is the exogenous 1 × K vector of covariates (or predictor variables) measured at time t, β is the K × 1 vector of predictor coefficients, and ψit is the ith dynamic cycle component which is modeled as a stationary process for i = 1, …, M, where M is the number of cycles in the model (in our case M = 6). The model specification for the cycle component is given by ψi,t+1 = ρi cos(λi)ψit + ρi sin(λi)ψit + ωit with the auxiliary dynamic process given by , where ρi is the autoregressive coefficient (determines the persistence of the cycle process), λi is the frequency of the cycle measured in radians, and ωit and are two IID noise terms, which are independent of each other, and all other noise terms, but they have the same (common) variance , for i = 1, …, M. It can be shown that we can formulate ψit as a stationary autoregressive moving average (ARMA) process. As long as the coefficient pairs (ρi, λi) are sufficiently different for different i = 1, …, M, the M cycle components ψ1,t, …, ψM,t can be uniquely extracted from the observed time series yt. The parameter constraints for each cycle process are 0 < ρi < 1 (stationarity), 0 < λi < 2π (circularity), and , for i = 1, …, M. The signal-to-noise coefficient for the ith cycle process is given by , for i = 1, …, M. The complete model for yt (in our case for the Niño-3.4 temperature time series) can be represented as a linear Gaussian state space model such that the Kalman filter methods can be used in a similar way as for the local level model. The dynamic level and cycle (including the auxiliary cycle variables) components are placed in the state vector, denoted by αt, which is subject to a multivariate dynamic stochastic process. The predictor coefficients in vector β are treated as time-invariant, fixed parameters. Both αt and β are simultaneously estimated as part of the Kalman filter [see Harvey (1989) and Durbin and Koopman (2012, Part I therein) for its general treatment]. Also in this more general context of the state space model, the Kalman filter methods remain to provide the MMSE optimal weights to the observations for signal extraction and forecasting.
The statistical dynamic components model can be viewed as a linear time series model with time-varying parameters. The introduction of time-varying parameters as done with stochastically evolving level and cycle components can address and approximate nonlinear features in the time series via piecewise linearization. The number of nodes for the linearization (or the smoothness of the piecewise approximation) is implicitly determined via the signal-to-noise parameters of the time-varying components. We therefore may claim that the introduction of the dynamic components also makes the analysis more robust to nonlinear features in the time series.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-18-0877.s1.