## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Predicting stochastic systems by noise sampling, and application to the El Niño-Southern Oscillation

Edited by* Alexandre J. Chorin, University of California, Berkeley, CA, and approved May 25, 2011 (received for review October 21, 2010)

## Abstract

Interannual and interdecadal prediction are major challenges of climate dynamics. In this article we develop a prediction method for climate processes that exhibit low-frequency variability (LFV). The method constructs a nonlinear stochastic model from past observations and estimates a path of the “weather” noise that drives this model over previous finite-time windows. The method has two steps: (*i*) select noise samples—or “snippets”—from the past noise, which have forced the system during short-time intervals that resemble the LFV phase just preceding the currently observed state; and (*ii*) use these snippets to drive the system from the current state into the future. The method is placed in the framework of pathwise linear-response theory and is then applied to an El Niño–Southern Oscillation (ENSO) model derived by the empirical model reduction (EMR) methodology; this nonlinear model has 40 coupled, slow, and fast variables. The domain of validity of this forecasting procedure depends on the nature of the system’s pathwise response; it is shown numerically that the ENSO model’s response is linear on interannual time scales. As a result, the method’s skill at a 6- to 16-month lead is highly competitive when compared with currently used dynamic and statistic prediction methods for the Niño-3 index and the global sea surface temperature field.

ENSO forecasting has a decade-long history and relies mainly on two classes of models: dynamical and statistical (1, 2). Still, a further distinction has to be made within the latter class: Some of the statistical models do not make any use of dynamical information, like Lorenz’s method of analogues (3) and its followers (4–6), while others do use a dynamical model—previously fitted to the observations from the past—to drive the statistics in the future (2, 7, 8). Empirical stochastic models belong to this hybrid category, and linear versions of such models have been used in ENSO forecasting for two decades; see ref. 9 for a survey. More recently, Kravtsov et al. (10) have extended this approach to nonlinear models by developing an EMR methodology that can include quadratic nonlinearities as well as state-dependent noise that parameterizes small-scale effects, without assuming a priori scale separation (11).

The purpose of this paper is to show that, under suitable circumstances, a better understanding of the role of the fast processes, weather or noise, can help predict the slow ones—namely, the climate. To achieve this purpose, we proceed in two steps: (*i*) develop a special prediction methodology, called past noise forecasting (PNF), using EMR models; and (*ii*) provide a theoretical framework for applying the PNF method—or any other forecasting method based on perturbations of the noise—to other empirical stochastic models. Of late, probabilistic forecasts in weather and climate prediction have become fairly widespread: They are grounded in an estimation of the probability density function (PDF: 12–14).

We take here a distinct, pathwise approach instead, and will show that this approach is particularly well adapted to empirical stochastic models and to phenomena in which a considerable part of the variability exhibits some form of repetitive regularity (15). We focus mainly on the EMR models of ENSO introduced and studied by Kondrashov et al. (16).

The paper is organized as follows. We present in the next section the concept of *pathwise sensitivity*, which describes the response of the stochastic system at hand with respect to perturbations of the noise realization. If this pathwise response is linear, then the system’s response to small perturbations of the noise may be small. In the subsequent section, we verify numerically that this linear-response property is fulfilled by the EMR model of ENSO used here. These two sections provide the theoretical framework for our PNF method, which is outlined in the fourth section.

The method’s key idea consists of selecting appropriate noise samplings—or snippets—from the past of the noise realization derived during the EMR model’s fitting procedure, and then using a subset of these snippets to drive the model’s dynamics into the future. An effective subset of snippets is selected by relying on ENSO’s well-known LFV modes (17) to key the selection to the analogous phases of these modes.

Essentially, the PNF method is a forecasting method based on perturbation techniques of the actual path of the future noise. The effects of changing the magnitude of these perturbations are analyzed by relying on the pathwise linear-response property. This property provides, in the fifth section, a clear theoretical explanation of the relatively good performance of the PNF method in the case at hand, as illustrated further by the statistical-significance tests in *SI Text*. The PNF approach is quite flexible, and further extensions are suggested in *Concluding Remarks*.

## Predictability and Pathwise Sensitivity

Random forcing in a climate model may increase the robustness of its long-term behavior (18), but it adds another source of uncertainty to short- or medium-term prediction. In stochastic linear models of the form *d***x** = **Ax***dt* + *dξ*_{t} (9)—where **x** is the state vector, *t* is time, and *ξ*_{t} is a vector-valued Wiener process—sensitivity to the initial state may be observed in the short run. In the long run, however, because stationarity of the stochastic process requires the linear operator **A** to be exponentially stable, there is no such sensitivity: To the contrary, trajectories started from different initial states **x**_{0} and driven by the same noise realization *dξ*_{t}(*ω*) are completely synchronized. Such *pathwise synchronization* characterizes all stochastic dynamical systems having only negative Lyapunov exponents (18, 19). We shall see that a weaker form of synchronization, namely so-called on-off synchronization, characterizes certain systems that do have a positive exponent (19).

Assuming the driving noise *dξ*_{t} is white, Ornstein–Uhlenbeck theory can be used to estimate **A** and the covariance of *dξ*_{t}. The Green’s function of the process can then be used for predicting the expected evolution of the state **x**(*t*) (9). Still, there is no information on the future of the sample path *t*↦*dξ*_{t}(*ω*) that has been fitted from the observations over a time window (0,*t*^{∗}); where *ω*∈Ω marks the sampled realization, and Ω is the appropriate probability space for the sample paths.

For a nonlinear stochastic model that exhibits a positive Lyapunov exponent when the noise is turned off, we are thus faced—once the noise is turned on and besides the sensitivity to the initial state—with the additional uncertainty in the realization of the noise. Hence the forecasting problem becomes, at first sight, even more difficult than in the purely deterministic case. In the remainder of this section, we show that for certain chaotic systems—which are of geophysical interest, as seen in the next section—the noise can help, rather than hinder, the prediction of the state.

We consider here a nonlinear stochastic model of the form, [1]the model is used to simulate a multivariate time series over a time interval (0,*t*^{∗}) and to predict its future evolution over the interval (*t*^{∗}, *t*^{∗} + *T*). In Eq. **1** and the sequel, the noise *dξ*_{t} is not necessarily white and hereafter will be denoted by *ξ*_{t}. One assumes that the deterministic functions **f** and **g** have been estimated, by EMR or in some other way. Unless one knows, though, how to approximate the sample path of the noise *t*↦*ξ*_{t}(*ω*) for *t*∈(*t*^{∗},*t*^{∗} + *T*)—with *ω* the unknown but fixed realization—one has to rely on predictions of the system’s PDF (12–14).

To exploit knowledge of past noise, as we propose to do here, requires first an estimate of the model’s response to a “perturbation” of *ω*. Consider an *observable* *ψ* of the system governed by Eq. **1**—i.e., a real-valued continuous function , where *X* is the system’s phase-space. For definiteness and clarity, let *ξ*_{t}(*ω*^{′}) be a perturbation of *ξ*_{t} over (0,*t*^{∗}), with *ω*^{′} a different realization and 0 < *ϵ* ≪ 1; the perturbed path of the noise is , although much more general perturbations can be considered. Denote by Φ(*t*,*s*,**x**_{0}; *ω*) the solution of Eq. **1** emanating from **x**_{0} at time *s* < *t*, when the system is driven by the path *t*↦*ξ*_{t}(*ω*), and by Φ_{ϵ}(*t*,*s*,**x**_{0}; *ω*^{′}) the one driven by and still emanating from **x**(*s*) = **x**_{0}.

Consider now a perturbation applied over the interval (*s*,*t*). We define the (local) deviation of *ψ* at *t*, when starting from **x**_{0} at time *s* and driven by the two noise paths and *ξ*_{t}(*ω*), respectively, from *s* to *t*: [2]*ω* here refers to the realization that drives the unperturbed system. The corresponding mean response at time *t* is provided by the expectation ; when no confusion is possible, we drop the indexing over *ω*^{′}∈Ω. We will be mainly interested in the expected response, averaged over (*s*,*t*) and given by , where .

At this point, we take a brief excursion beyond the scope of the present paper and note that, by taking an ensemble average over **x**_{0}∈*X* and letting *s* → -∞, one recovers, in our stochastic context, quantities that are analogous to those considered in Ruelle’s (20) response theory for smooth, time-dependent perturbations *F*(**x**,*t*) of autonomous systems with chaotic behavior. In that theory, the nature of the response—whether linear or nonlinear—is independent of the respective choices of the observable *ψ* and the perturbation in the (deterministic) forcing *F*.

There is no room to discuss here a rigorous framework in which to assess the response of a stochastic system governed by Eq. **1** to perturbations of the noise path *t*↦*ξ*_{t}(*ω*). Only numerical results for the short-range response can be given in the present paper; these will be justified rigorously elsewhere. Applications of linear-response theory to climate sensitivity have been recently investigated in a deterministic (21) and in a stochastic context (22), but only for the model’s PDF, while pathwise statistics are discussed in ref. 19.

To which extent does the local, short-range response depend on **x**_{0}, *ω*, *s* and *t*? It appears that the linear or nonlinear nature of , over a fixed interval (*s*,*t*), is unchanged for almost all **x**_{0} and almost all *ω*; as a consequence, we will drop the dependence in **x**_{0} and *ω*. This statement holds, for instance, in the case of the stochastic Lorenz system forced by *σ***x***dξ*_{t}, with *σ* > 0, as considered in ref. 19: For *t* - *s* fixed—and almost surely in **x**_{0} and *ω*— depends nonlinearly on the perturbation; this nonlinear response sets in above a certain threshold *ϵ*_{0} of the perturbation size, which is much smaller than the variance of the noise.

To the contrary, other nonlinear (and chaotic) stochastic systems can exhibit—over a wide range of perturbations and still for almost all **x**_{0} and *ω*—a linear response. For *t* - *s* fixed, such a response is visualized through a slope of vs. *ϵ* that is constant or changes only slightly with **x**_{0} and *s*. In certain cases, a loss of this linear dependence may occur above a relatively large perturbation size. In practice, the knowledge of , and the threshold *ϵ*_{0} past which the response becomes nonlinear, gives a quantitative assessment of the pathwise sensitivity: The smaller *ϵ*_{0} the more sensitive the system may be to a perturbation *ϵξ*_{t}(*ω*^{′}) of the path. It follows also that a “reasonable” value of *ϵ* ≤ *ϵ*_{0} could serve as a good indicator of an admissible noise perturbation. In the next section, we illustrate what is reasonable for an EMR-model of ENSO.

## Pathwise Linear Response of an EMR-ENSO Model

EMR models (10, 11, 16) are a subset of the stochastic systems described by Eq. **1**; they can be compactly written as [3]Here **x** represents the slowest and most energetic modes, **B**(**x**,**x**) is a quadratic nonlinearity, and **L** is a time-dependent linear operator obtained by integrating recursively from the *L*th level to the top level, *l* = 0, the linear stochastic equations that relate the auxiliary stochastic forcings and —i.e., ; **A** and **B** are estimated by a least-square procedure, and **b**_{l} are linear maps estimated recursively along the same lines. The procedure is stopped when , called the *L*th-level residual forcing, has a lag-1 vanishing autocorrelation. The stochastic forcings are ordered from the one with strongest memory, , to the most weakly autocorrelated one, *ξ*_{t}. Here and throughout this section, *ξ*_{t} in Eq. **3** is to be understood thus as an approximation over the sampling interval of *dξ*_{t} in Eq. **1**. This procedure allows one to parameterize the “fast” modes in terms of the “slow” ones; see *SI Text* for further details.

Kondrashov et al. (16) have shown that a two-level EMR model—i.e., *L* = 2 in Eq. **3**—can simulate key features of the global sea surface temperature (SST) field’s LFV and is quite competitive in predicting ENSO events on the seasonal-to-interannual scale. Let us assume that such a model has been fitted to the SST observations over some interval (0,*T*_{0}) and measure its pathwise sensitivity to changes in the model’s “weather”; i.e., in the EMR procedure’s residual noise .

Consider now the observable *ψ*(**x**)≔‖**x**‖, where ‖·‖ denotes the Euclidean norm of **x**. Our goal is to obtain an estimate of in terms of *ϵ*; to do so, we calculate [4]as a function of . Here (*t* - *s*) ≪ *T*_{0}, Φ(*t*) is the reference solution driven by the reference path of the noise *ξ*_{t}(*ω*), *σ*_{(0,T0)}(‖*ξ*‖) is the standard deviation on (0,*T*_{0}) of the norm of the residual noise *ξ*_{t}(*ω*) at the second level, while *σ*_{0,T0}(‖Φ‖) is the standard deviation of the norm ‖Φ‖ on (0,*T*_{0}).

The initial state at time *s* is chosen as a representative **x**_{0}, associated with a state found in the past of the SST field, and *t* is chosen so that *t* - *s* = *T* = 16 months, while the state of the system is known up to time *t*. The results are depicted as the blue points in Fig. 1, which clearly demonstrates a pathwise linear response, with a slope of about 0.8 for our two-level EMR model (blue curve in the figure) over a large range of perturbations of the noise.

This stochastic EMR model has a leading Lyapunov exponent that is positive, and it exhibits the intermittency phenomenon described by Chekroun et al. (19). This phenomenon is characterized by on-off synchronization of the trajectories for almost all fixed realizations; it is also associated with a sensitive dependence to initial states that is considerably weaker than observed in a stochastically perturbed Lorenz system (19).

In order to perform a more detailed analysis of the sensitivity to initial states of the ENSO model used here, we propose a simple methodology that is consistent with the forecasting method developed in this paper (see below). To be precise, we have considered 60 reference initial states spread over 60 different initial times *s*; to each of these 60 states, we have applied 120 perturbations corresponding each to a state occupied by the system over 10 yr prior to *s*. The system was driven by the same realization *ω* of the noise, and the divergences between the trajectories emanating from the reference initial state **x**_{0} and the 120 perturbed ones , have been plotted in Fig. 2 *A* and *B* at 1 mo (*t* - *s* = 1, blue points) and 16 mo (*t* - *s* = 16, red points), respectively.

The divergence at 1 mo exhibits a linear dependence with respect to the magnitude of the perturbation (*A*), but this dependence becomes nonlinear at 16 mo (*B*). Fig. 2*B*, however, clearly demonstrates that two distinct clusters of red points may be identified: Either the model responds actively to a perturbation of the initial state, and the scatter is large, or the response is weak; in the latter case, the divergence is less than 0.25—i.e., less than 25% of the SST field’s standard deviation. Fig. 2*B* shows, furthermore, that the “mass” of this compact cluster is larger than the mass of the former cluster.

This feature demonstrates that, at 16 mo, the trajectories have a stronger tendency to synchronize than to diverge, which indicates on-off synchronization with a stronger on mode than the off mode (19). It thus appears that a fairly strong tendency of trajectories to synchronize for given, fixed stochastic forcing at intermediate range is associated with a weaker dependence on perturbations of this forcing, on the one hand, and with smaller divergence of trajectories at short times, on the other.

The results in Figs. 1 and 2 are thus fairly encouraging in terms of developing a predictive methodology based on estimates of future noise, because errors in such estimates do not seem to produce unduly large errors in the forecast. We develop such a methodology in the next section.

## Past Noise Forecasting (PNF)

### Sampling Stationary Stochastic Processes.

A well-known fact from the theory of stochastic processes is that—given a probability space associated with a stochastic process {*t*↦*ξ*_{t}(*ω*)}_{ω∈Ω}, with the *σ*-algebra of events and the probability measure—a time-parameterization of Ω is obtained by introducing the shift operator *θ*_{t}*ξ*_{s}(*ω*)≔*ξ*_{t+s}(*ω*) and letting *ξ*_{s}(*ω*) = *ω*(*s*); see refs. 18 and 19 and references therein. If the stochastic process is stationary or with stationary increments, then the noise statistics are invariant under the transformation *θ*_{t} or its helix analogue (19)—i.e., .

It follows that shifting in time the path of such a stochastic process *t*↦*ξ*_{t}(*ω*), given for a fixed realization *ω*, leads to another path, *t*↦*ξ*_{t}(*ω*^{′}), of the same stochastic process but for another realization *ω*^{′}. An immediate consequence is that, given a path *t*↦*ξ*_{t}(*ω*) of length *t*^{∗} > 0, each continuous segment or snippet of length *T* of this path, with *T* < *t*^{∗}, provides another path *t*↦*ξ*_{t}(*ω*^{′}) of length *T*, for a realization *ω*^{′} of the same stochastic process. This property allows one, therefore, to select several sample paths of length *T*—associated with other realizations {*ω*^{′}}—from a single path of length *t*^{∗}≫*T*.

### Noise Sampling Refinement Based on LFV Analogs.

We consider now 0 < *T* ≪ *t*^{∗} and select *N* continuous, possibly overlapping segments of length *T* each, from a time series of the noise *ξ*_{t} of length *t*^{∗}, given for a fixed realization *ω*. Note that *N* ≤ *t*^{∗} depends on *t*^{∗}, and will be denoted by *N*_{t∗} when necessary. Starting the copying at time *t* = 0, this procedure leads to a set of *N* snippets of size *T*, denoted by , with *t*_{i}∈[0,*t*^{∗} - *T*], *t*_{0} = 0, and *t*_{i} < *t*_{i+1}. Each of these constitutes a path of length *T* of the noise for a realization *ω*_{i} that could potentially drive, from time *t*^{∗} to time *t*^{∗} + *T*, any stochastic system driven by the noise *ξ*_{t}; note that *ω*_{0} = *ω*.

First, we notice that—by running the two-level EMR–ENSO model of ref. 16 from an arbitrary *t*^{∗}, and driven by —the average over a large enough number *N*_{t∗} of simulations is essentially the same as the one obtained from a large ensemble of arbitrary realizations of the residual noise *ξ*_{t}, as given by Eq. **3**. This numerical observation points to two facts: (*i*) The ensemble contains sufficient statistics to recover the mean obtained from a larger ensemble; and (*ii*) no improvement is obtained with respect to forecasting by the mean, thus motivating the need to refine the set of snippets in some way. In other words, we ask if there exists a subset of that allows one to improve the forecast by simply taking the mean over this subset.

Because is derived merely by sampling the stochastic forcing from the past, we examine next the phase of the system that corresponds to a snippet , in order to select those snippets that correspond to the phase of the system at time *t*^{∗}. This selection is motivated by the fact that, for the EMR–ENSO model at hand, sensitivity with respect to the initial state is linear or weak on time scales of *T* = 1 - 16 months (cf. Fig. 2), and therefore the system’s phase is mainly determined by the stochastic forcing over these time scales. We thus select noise samples of length *T* from the past that have forced the system when it occupied a phase that resembles the one just preceding the currently observed state at *t*^{∗}, and use these to drive the system from *t*^{∗} on till *t*^{∗} + *T*. The aggregate forecast skill by relying on the individual forecasts so obtained will depend mainly upon the distribution of the selected snippets as perturbations of the actual path of the forcing into the future, provided the system exhibits a linear pathwise response to such perturbations. The latter point is discussed below, just before *Concluding Remarks*.

We still have to define the precise meaning we attribute to the system’s being in a given phase. In fact, several methods exist to examine the instantaneous phase of a stochastic, as well as a deterministic, system. In order to refine the set by considering the phase of our two-level EMR model of ENSO, we adopt here a heuristic approach and look at the phase of a smoothed version of the SST field’s time series. The smoothing is obtained by singular spectrum analysis (SSA) (23) and turns out to provide a sharp enough selection criterion for the snippets of interest in our PNF.

To simplify the discussion, we only consider here *x*_{1}—i.e., the first component of our EMR model—which gives the time evolution of the leading principal component (PC1) of the SST field (16) and captures about 80% of the variance (24). The idea of smoothing PC1 by SSA relies on the fact that the subset formed by the leading reconstructed components (RCs) of PC1 captures the LFV of this time series, while filtering out the spurious high frequencies. Ghil and Jiang (2) showed that much of ENSO prediction skill does lie in correctly forecasting its quasi-quadrennial (QQ) and quasi-biennial (QB) modes (25, 26).

Our approach consists of finding analogues for a continuous segment of a selected RC—of length 0 < Δ ≤ *T* months and ending at *t*^{∗}, the epoch from which we want to predict—by scanning the record of that RC over the available time series. The closeness of our analogues is determined by using both the rms (RMS) and the correlation coefficient (Corr) between two such segments. Consider the sum RC_{K} of the *k* = 1,…,*K* leading RCs of PC1 (cf. ref. 23 and *SI Text*). The selection procedure for analogues that correspond to prediction from an epoch *t*^{∗} is based on the following two steps:

Step S1. Let RC

_{K}(*t*_{j},*t*_{j}+ Δ) be the segment of RC_{K}that starts at*t*_{j}and has length Δ. We first select the times*t*_{j}<*t*^{∗}- Δ such that RC_{K}(*t*_{j},*t*_{j}+ Δ) is close enough to RC_{K}(*t*^{∗}- Δ,*t*^{∗}) in both RMS and Corr. More precisely, we compute the set: [5]with prescribed parameters 0 <*α*,*β*≪ 1. In theory, could be the empty set for very small α and β. In practice, α and β have to be tuned so to avoid this, while to ensure that is a strict subset of ; see step 2.Step S2. Define next the following subset of , [6]

The set is entirely based on the past observations and will serve us to drive the system into the future, as follows.

For any noise snippet , we denote by the solution of our two-level EMR model of ENSO at time *t*^{∗} + *t* starting from the current observation **x**_{t∗} at time *t*^{∗}, and forced by over the interval (*t*^{∗}, *t*^{∗} + *t*). Note that, by construction, 0 < *t* ≤ *T*.

Our predictor—from *t*^{∗} for *t* months ahead—will be: [7]where *μ*_{j}∈[0,1] and . In the sequel, we simply illustrate the method with for all *j*, and postpone the discussion of optimizing the choice of the *μ*_{j}s for future work. Note that, in Eq. **7**, can be replaced by a shifted *ζ*^{tj+s}, with 0 < *s* ≤ Δ, because we consider only an “average phase” of the system over Δ-months, rather than an instantaneous phase.

We call any method that consists of using Eq. **7** as a predictor, and steps **S1** and **S2** for the selection procedure of the noise snippets a PNF method; in other words, EMR fitting is but one way to determine the driving noise *ζ*^{t} from which the snippets are selected. A schematic illustration of the concept is given in *SI Text*.

In practice, the parameters *K* and Δ—needed to fully specify a PNF (*K*,Δ) method—have to be determined either by tuning, as discussed in *SI Text* or by a more systematic optimization method. PNF methods are being presented here in the context of ENSO and for a two-level EMR model, but they are quite general and can be used, in principle, for any empirical stochastic model fitted to some time-evolving dataset. Their forecasting skill depends, of course, on the nature of the model’s pathwise response introduced in the previous section; see also the last-but-one section.

## Numerical Results and Forecast Skill

We use here a two-level EMR model—i.e., *L* = 2 in Eq. **3**—of ENSO variability (16), embedded in the subspace of the 20 leading PCs of a time series of SST anomalies that is -month long (January 1950–September 1998), given on a 5° × 5° grid over the 30 S–60 N latitude band (27) and has the annual cycle removed. We predicted starting on successive months , with *δt*^{∗} = 1 month, and *q*∈{0,…,*Q* - 1}, with *Q* = 114, via (**7**), and the maximum number of months over which we do so is fixed at *T* = 16 months. Note that the total number of possible snippets here is the same for each —i.e., we do not update the EMR model from one value of to another; keeping it fixed only makes the test for the PNF method even harder.

For Δ = 5 months and using selection criterion **S1**—applied to RC_{K} of PC1, with *K* = 2, *α* = 0.5 and *β* = 0.05 in Eq. **5**—the PNF(2,5) method selects snippets that allow us to reach the goal set in the preceding section: improve the forecast obtained by taking the mean of the EMR model by refining the subsets . Fig. 3*A* illustrates this improvement in predicting the Niño-3 index by comparing the mean of the EMR model prediction (red curve) with the predictor given by (**7**) (blue curve). These predictions are issued each month *t*^{∗} using only data prior to that month, for *Q* = 114 values of , running from October 1998 to March 2008.

We note in Fig. 3*A* that the PNF method is able to capture episodes of strong anomalies in the evolution of the Niño-3 index—especially the large El Niño in 2003 and the large La Niña in 2008—much better than the mean EMR method does. While the PNF performance is not uniformly better, it is at no time substantially worse that the mean EMR either. The PNF improvement is most striking during the energetic phases of constructive interference between the LFV modes, QQ and QB (2, 26).

Fig. 3 *B* and *C* shows that the PNF method yields significantly better skill in Niño-3 prediction beyond 6 mo, compared with the standard EMR method of Kondrashov et al. (16). The optimal improvement occurs at relatively long lead times of 12–16 mo. This improvement is due, in part, to the fact that, in our EMR–ENSO model, the chaotic behavior is relatively weak and most trajectories, while starting from different initial states, synchronize for a fixed realization of the noise at longer lead times, as shown in Fig. 2*B*. The PNF improvement is also consistent with the characteristic decay time *τ* of the ENSO eigenmode associated with its most energetic LFV mode, the QQ mode; this decay time of *τ* ≈ 14 months (16) might ultimately determine the empirical limits of long-term ENSO predictability.

Even more strikingly, Fig. 4 shows that PNF skill in predicting the SST field itself at a 14-mo lead is uniformly better over the entire equatorial Tropical Pacific and Indian Ocean area, where ENSO is most active. This improvement is manifest in both forecast skill measures, Corr (*A* and *B*) as well as RMS (C and D).

In order to understand better the advantage of selecting appropriate noise snippets to drive the PNF forecasts, we compare the results in Fig. 3 *B* and *C* with those obtained by using a reshuffled version (green) of the residual noise *ξ*_{t}, obtained at the second level. The reshuffling modifies the autocorrelations of *ξ*_{t} but changes neither its mean nor its variance. This comparison demonstrates that the EMR forecasts with reshuffled noise are identical to those with the original EMR—i.e., the ensemble mean, over snippets, of the future paths is insensitive to random permutations of the residual noise.

To the contrary, the PNF(2,5) predictor in Eq. **7** is sensitive to being driven by snippets, and the skill of this reshuffled PNF drops dramatically, because it does not respect the relation between the driving noise and the LFV phases in the past. Statistical-significance tests in *SI Text* buttress the numerical results on PNF and EMR prediction skill presented here. Furthermore, the PNF methodology is illustrated in *SI Text* by introducing and analyzing a simple ecological “toy model” forced by white noise.

## PNF Method and Pathwise Linear Response

It remains to show that the relatively good skill obtained by the PNF method is in agreement with the pathwise sensitivity results obtained for the two-level EMR model of ENSO (cf. Fig. 1) and the fact that sensitivity with respect to initial data is weak at longer lead times (cf. Fig. 2*B*). Going back to Fig. 1, we have plotted in magenta and for the 114 epochs , *q* = 1,…,*Q* used in obtaining Fig. 3, the corresponding outputs, , with . The red dots are for the outputs driven by all the elements of , with no selection of a subset .

Clearly the PNF predictors (magenta) form a substantially smaller cluster than the EMR predictors (red cluster): The red outliers correspond to the strongest response and are associated with the worst forecasting skill. We observe furthermore in Fig. 1 that the magenta cluster is located within the linear-response regime, fairly close to a slope of about 0.8, but relatively far from the origin, for almost all initial states.

This result illustrates three facts: (*i*) Our LFV-conditioned PNF method selects snippets that correspond to small enough perturbations of the future noise realization {*ζ*_{t}: *t*^{∗} < *t* ≤ *t*^{∗} + *T*}, and are thus associated with a pathwise linear response of the system; (*ii*) still, this method does not provide paths that are close enough, in general, to *ζ*_{t} in order to allow one to recover the system’s actual future evolution; and (*iii*) certain extreme episodes are relatively well captured. Note, however, that we used here a relatively rough version of the PNF method, because the purpose was just to demonstrate the usefulness of considering past noise in forecasting future system evolution.

## Concluding Remarks

In this article, we have developed a forecasting method for dynamical systems that exhibit (*i*) oscillatory LFV modes and (*ii*) noise effects. The method relies on the use of an empirical stochastic model in which a path of the noise is estimated over a finite-time window. This method uses only information from past observations, and it was illustrated by using an EMR model of ENSO. We showed that the domain of validity of such a method depends strongly on the nature of the system’s pathwise response, and that this response is both weak-to-moderate and linear at finite times for the model at hand; see Figs. 1 and 2. The fact that even the “approximately right” noise can help, rather than hinder, prediction is surprising, to say the least.

The forecast skill of the method was evaluated in terms of rms error RMS and anomaly correlations Corr, for lead times of 6–16 months, against a currently operational method, the standard EMR prediction. The latter has been tested over the last 5 yr against 15 dynamical and another 7 statistical methods, and has proven quite competitive at lead times of 6–12 mo, in the context of ENSO prediction at the International Research Institute for climate and society (IRI), http://iri.columbia.edu/climate/ENSO/currentinfo/.

The PNF method seems to improve on the EMR method at 12 mo and beyond (Figs. 3 and 4). Further improvements seem possible and include use of multivariate, rather than univariate, SSA, and better criteria for choosing the most appropriate weather noise by LFV phase. These practical aspects—along with the very intriguing theoretical ones associated with the method’s basic concepts—could present interesting topics for future studies.

## Acknowledgments

We are grateful for discussions with and encouragement by S. Kravtsov, J. C. McWilliams, J. D. Neelin, A. W. Robertson, E. Simonnet, and I. Zaliapin. Two anonymous referees provided constructive and insightful suggestions. This study was supported by US Department of Energy grants DE-FG02-07ER64439 and DE-FG02-02ER63413 and by National Science Foundation grant DMS-1049253.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: mchekroun{at}atmos.ucla.edu.

Author contributions: M.D.C., D.K., and M.G. designed research; M.D.C. and D.K. performed research; M.D.C. and D.K. contributed new analytic tools; M.D.C., D.K., and M.G. analyzed data; and M.D.C., D.K., and M.G. wrote the paper.

The authors declare no conflict of interest.

*This Direct Submission article had a prearranged editor.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1015753108/-/DCSupplemental.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Van den Dool HM

- ↵
- ↵
- Keppenne CL,
- Ghil M

- ↵
- ↵
- ↵
- Palmer TN,
- Williams P

- Kravtsov S,
- Kondrashov D,
- Ghil M

- ↵
- ↵
- ↵
- Kleeman R,
- Majda A,
- Timofeyev I

- ↵
- Ghil M,
- Robertson AW

- ↵
- ↵
- ↵
- ↵
- Chekroun MD,
- Simonnet E,
- Ghil M

- ↵
- ↵
- ↵
- Majda AJ,
- Wang X

- ↵
- Ghil M,
- et al.

- ↵
- Preisendorfer RW

- ↵
- ↵
- Jiang N,
- Ghil M,
- Neelin D

- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Applied Mathematics