Pascal Yiou

Didier Sornette

Michael Ghil

Data-Adaptive Wavelets and Multi-Scale SSA

Using multi-scale ideas from wavelet analysis, we extend singular-spectrum analysis (SSA) to the study of nonstationary time series of length N whose intermittency can give rise to the divergence of their variance. The wavelet transform is a kind of local Fourier transform within a finite moving window whose width W , proportional to the major period of interest, is varied to explore a broad range of such periods. SSA, on the other hand, relies on the construction of the lag-covariance matrix C on M lagged copies of the time series over a fixed window width W to detect the regular part of the variability in that window in terms of the minimal number of oscillatory components; here W = M ∆t, with ∆t the time step. The proposed multi-scale SSA is a local SSA analysis within a moving window of width M ≤ W ≤ N . Multi-scale SSA varies W , while keeping a fixed W/M ratio, and uses the eigenvectors of the corresponding lagcovariance matrix CM as a data-adaptive wavelets; successive eigenvectors of CM correspond approximately to successive derivatives of the first mother wavelet in standard wavelet analysis. Multi-scale SSA thus solves objectively the delicate problem of optimizing the analyzing wavelet in the time-frequency domain, by a suitable localization of the signal's covariance matrix. We present several examples of application to synthetic signals with fractal or power-law behavior which mimic selected features of certain climatic and geophysical time series. A real application is to the Southern Oscillation index (SOI) monthly values for 1933-1996. Our methodology highlights an abrupt periodicity shift in the SOI near 1960. This abrupt shift between 4 and 3 years supports the Devil's staircase scenario for the El Niño/Southern Oscillation phenomenon. 1

Introduction and motivation

Intermittent time behavior is common in a large variety of systems in the natural and socio-economic sciences. Multiple regimes of behavior and the transitions between them are ubiquitous in climate dynamics, from the alternation of glacials and interglacials to that of atmospheric blocking [START_REF] Ghil | Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory and Climate Dynamics[END_REF]. River discharge and rainfall similarly exhibit long periods of relative quiescence interspersed with large floods and bursts [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF]. The solid earth is characterized by even stronger intermittency that manifests itself by long tails of earthquake-size distribution with strong clustering of seismic activity in time and space [START_REF] Kagan | Observational evidence for earthquakes as a nonlinear dynamic process[END_REF]. Even mantle convection is now believed to be intensely intermittent, with internal boundary-layer instability leading to "avalanches" that come in a large range of sizes and follow a power-law distribution [START_REF] Peltier | Phase transition modulated mixing in the mantle of the earth[END_REF]. The number of contagious-disease incidences is also intermittent, and has recently been documented to be distributed according to power laws [START_REF] Rhodes | Power laws governing epidemics in isolated populations[END_REF]. The evolutionary path of species on earth with its "punctuated dynamics" provides a paradigm of intermittent behavior [START_REF] Gould | Full House[END_REF]. Economic and financial time series have long been recognized to have nonstationary variances, possibly described by Lévy processes [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF].

These are only a few examples among many that exhibit the ubiquity of intermittent behavior, and hence of extreme events. Additional examples range from volcanic eruptions, hurricanes and tornadoes, landslides, avalanches, lightning strikes, and meteorite and asteroid impacts, to the failure of engineering structures, crashes in the stock market, social unrest leading to large-scale strikes and upheaval, economic drawdowns on national and global scales, regional power blackouts, and traffic gridlock.

Many of these problems involve nonstationary time series, while standard methods of time-series analysis assume explicitly or implicitly stationarity. Stationarity in the wide sense is defined by (i) the existence of first (mean) and second (variance) moments of the time series and (ii) their invariance with respect to time translations. The failure of either assumption can lead to spurious spectral-analysis results.

Wide-sense stationarity implies that the covariance operator of a time series exists and is a function of the separation between epochs only, i.e., of "lags." The spectral decomposition of the resulting lag-covariance operator is connected to the power spectrum of the time series through the Wiener-Khinchin theorem [START_REF] Hannan | Time Series Analysis[END_REF][START_REF] Percival | Spectral Analysis for Physical Applications[END_REF]. For time series in continuous time and of infinite length, the corresponding eigenfunctions are always complex exponentials or sine-and-cosine pairs, and one time series differs from another solely through the associated eigenvalues.

For a time series of given, finite length -in discrete time -singular-spectrum analysis (SSA) [START_REF] Broomhead | Extracting qualitative dynamics from experimental data[END_REF][START_REF] Vautard | Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series[END_REF] relies on the Karhunen-Loève decomposition of an estimate of the covariance matrix that is based on M lagged copies of the time series. The eigenvectors of this matrix are often called empirical orthogonal functions, or EOFs and differ, in general, from sines and cosines [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF]. For the given time series, the EOFs form an optimal basis that is orthonormal at zero lag; they permit one to decompose the signal into its, possibly anharmonic, oscillatory components and aperiodic ones. The data-adaptive character of the SSA basis functions presents an advantage over the standard Fourier basis of the classical Blackman-Tukey [START_REF] Hannan | Time Series Analysis[END_REF][START_REF] Percival | Spectral Analysis for Physical Applications[END_REF] method: it allows one to capture a nonlinear oscillation by a single pair of EOFs -that maximizes the variance associated with the given frequency -rather than by a number of sine-and-cosine pairs associated with that basic frequency and its harmonics.

SSA seems to be less sensitive than classical Fourier analysis to deviations of a given time series from wide-sense stationarity assumptions. Trends are captured quite successfully, most often by the leading one or two EOFs of the given time series [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF][START_REF] Ghil | Interdecadal oscillations and the warming trend in global temperature time series[END_REF] and simple discontinuities in the local value of an oscillation's frequency can also be detected by SSA [START_REF] Yiou | Spectral analysis of climate data[END_REF][START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF].

Wavelet analysis, on the other hand, has become a tool of reference for intermittent, complex and self-similar signals, because it works as a mathematical microscope that can focus on a specific part of the signal to extract local structures and singularities [START_REF] Meyer | Ondelettes et Opérateurs I: Ondelettes[END_REF][START_REF] Daubechies | Ten Lectures on Wavelets[END_REF][START_REF] Meyer | Wavelets: Algorithms and Applications[END_REF]. The first step in the definition of a wavelet transform is the choice of the analyzing function, often called "mother wavelet." Considerable work has been devoted to finding mother wavelets that provide an orthogonal expansion that is complete but not (highly) redundant. It is also often desirable to concentrate the spectral energy in an optimal band for the problem at hand, while keeping the localization property.

A large number of analyzing wavelet functions have been introduced to satisfy these often conflicting requirements, with their relative advantages and drawbacks. To provide an optimal multi-scale ing wavelet to the signal rather than search through the extensive "libraries" of mother wavelets (http://www.mathsoft.com/wavelets.html) and (B) modify this shape in time and scale, especially if the data set's nonstationarity implies different structures as one slides along time or scale. The intrinsic constraint of a unique mother wavelet does not allow for this flexibility.

We introduce here a simple and powerful resolution to the limitations of both types of methodology. The standard forms of the wavelet transform and SSA are recalled in Secs. 2.1 and 2.2 respectively. The SSA method -modified to use varying windows with width W proportional to the order M of the lag-covariance matrix C -provides data-adaptive wavelet transforms with analyzing functions given by the leading EOFs of the matrix C; it is described in Sec. 2.3. A set of synthetic time series that exhibit intermittent and self-similar properties, to wit, a Cantor set, a log-periodic process, and a multiplicative noise, are introduced in Sec. 3.1; the details of their construction appear in Appendices 7.1-7.4. We also introduce a climatic time series, as a "real-world" test set. The multi-scale SSA (MS-SSA) of Sec. 2.3 is applied to these four test sets in Secs. 4.1 through 4.4, respectively.

MS-SSA provides a new approach for testing self-similarity and its breakdown by comparing the data-adaptive wavelets at different scales, while it retains the usual scaling of wavelet coefficients which is obtained from any wavelet transform [START_REF] Arneodo | Beyond classical multifractal analysis using wavelets: Uncovering a multiplicative process hidden in the geometrical complexity of diffusion limited aggregates[END_REF]. A summary of our results and comparison with complementary approaches to address the requirements (A) and (B) above appear in Sec. 5.

Methodology

Wavelet transform

A wavelet transform requires the choice of an analyzing function ψ, with general admissibility properties [START_REF] Meyer | Ondelettes et Opérateurs I: Ondelettes[END_REF][START_REF] Daubechies | Ten Lectures on Wavelets[END_REF][START_REF] Meyer | Wavelets: Algorithms and Applications[END_REF], and with the more specific property of time and frequency localization, i.e., ψ and its Fourier transform ψ must decay rapidly outside a given interval. Functions based on a Gaussian, f (x) = exp(-x 2 ), possess the localization property even if they do not verify the admissibility condition that the integral of the mother wavelet vanish. This zero-mean condition ensures that knowledge of all position-and scale-dependent wavelet coefficients allows one to retrieve the initial function through an inversion formula (similar to a two-parameter inverse Fourier transform) [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF].

In the sequel, we shall follow [START_REF] Arneodo | Beyond classical multifractal analysis using wavelets: Uncovering a multiplicative process hidden in the geometrical complexity of diffusion limited aggregates[END_REF] and use a Gaussian wavelet and its first derivative for the sake of simplicity. The Gaussian itself does not satisfy the admissibility condition, while its first derivative does. The overall statistical characterization of complex structures depends only weakly on the choice of the mother wavelet [START_REF] Arneodo | Beyond classical multifractal analysis using wavelets: Uncovering a multiplicative process hidden in the geometrical complexity of diffusion limited aggregates[END_REF] and we choose this simple example to illustrate its multi-scale low-pass filter properties and contrast them to the proposed MS-SSA technique.

A ψ-wavelet transform W ψ in continuous time and frequency is simply a projection of a signal

X(t), -∞ < t < ∞, onto b-translated and a-dilated versions of ψ: W ψ (a, b) = 1 √ a ∞ -∞ X(t)ψ t -b a dt (1) = 1 √ a ∞ -∞ X(t + b)ψ t a dt. (2)
If most of ψ is concentrated between, [-1, 1] say (up to a rescaling), then Eq. ( 1) is clearly an analysis of X in the interval [ba, b + a], where the integral is nonvanishing a priori.

Using the successive derivatives ψ (n) of a given wavelet ψ in Eq. ( 1) is equivalent to a ψ-analysis (up to a normalization factor) of the successive derivatives of the time series X; this is easy to see through an integration by parts. It is also easy to show, for a wide class of smooth analyzing functions ψ, that the number of oscillations in ψ (n) increases with n; this is true, in particular for our Gaussian wavelet of choice.

The original signal (or a filtered version of it) can be reconstructed from the family of wavelet transforms. Hence, for scale values a in an interval I, a reconstructed version X I of the signal X(t) is:

X I (t) = A ψ a∈I ∞ b=-∞ W ψ (a, b)ψ t -b a dadb a 2 ; ( 3
)
A ψ is a normalization factor which only depends on the mother wavelet ψ. This formulation is essentially a bandpass filter of X through I. If I is the positive real line, I = R + , then X I = X. Formula (3) applies when the mother wavelet is the derivative of the Gaussian, which we use below. As mentioned previously, the Gaussian ψ(x) = exp(-x 2 /2) itself cannot be used in this reconstruction formula, because it does not satisfy ψ(x)dx = 0. However, the forward transform of Eq. ( 2) is well-defined and provides nevertheless a useful multi-scale decomposition of the signal. The counterpart of Eqs. ( 1) and ( 3) for discrete time series can be obtained in a straightforward way. Given a discrete signal X i = X(t i ), t i = i∆t, 0 ≤ i ≤ N , there are numerous ways of estimating Eq. ( 1) numerically that avoid the problems of high-frequency noise and aliasing, while minimizing computational cost [START_REF] Meyer | Wavelets: Algorithms and Applications[END_REF][START_REF] Holschneider | Wavelets: An Analysis Tool[END_REF]. Computationally efficient algorithms based on the fast Fourier transform can be used [START_REF] Meyer | Wavelets: Algorithms and Applications[END_REF], but we will only need here the simplest numerical quadrature formulae to process Eq. (1).

Singular spectrum analysis (SSA)

SSA has been described theoretically, for time-continuous signals X(t), by Vautard and Ghil [START_REF] Vautard | Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series[END_REF]. Since it is still less familiar to a general readership than the wavelet transform, we summarize here its form for time-discrete signals {X i = X(i∆t), 0 ≤ i ≤ N }, which is generally used in practice, as it will be in the present study. Methodological and review papers include [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF][START_REF] Yiou | Spectral analysis of climate data[END_REF][START_REF] Dettinger | Software expedites singularspectrum analysis of noisy time series[END_REF][START_REF] Ghil | Spectral methods: What they can and cannot do for climatic time series[END_REF]. SSA proceeds by finding the eigenvectors, called EOFs, of the lag-covariance matrix C of the time series {X i } as above. The true lag-covariance matrix C of the process that generated {X i : 1 ≤ i ≤ N } is estimated by a lag-covariance matrix C M ; assuming stationarity of X yields a Toeplitz structure

C M,i = 1 N -i N -i j=1 X j X j+i , (4)
where i (0 ≤ i ≤ M -1) is the diagonal index of the matrix, and the time series X is normalized to have zero mean. The eigenvectors of this symmetric matrix are orthonormal and provide the Karhunen-Loève basis for expanding the time series X i with respect to its M lagged copies. These EOFs ρ k of C M are sorted by decreasing order of the associated eigenvalue

λ k , 0 ≤ λ M ≤ • • • ≤ λ 2 ≤ λ 1 . The principal component (PC) a k = {a ki : 0 ≤ i ≤ N -M } is computed by projecting the time series onto ρ k = {ρ kj : 0 ≤ j ≤ M }: a ki = M j=1 X i+j-1 ρ kj . (5)
For a given k, this equation is formally similar to Eq. ( 2), in which the integration interval is essentially restricted to the domain over which the wavelet function is numerically nonnegligible. Furthermore, the EOFs ρ k have an oscillation property which resembles that of the derivatives ψ (n) : the higher the order k, the more sign changes ρ k has [START_REF] Vautard | Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series[END_REF][START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF]. This oscillation property is due to the formal analogy between the Toeplitz structure of the covariance matrix C M and that of finitedifference discretizations of Sturm-Liouville problems on an interval W = M ∆t [START_REF] Courant | Methods of Mathematical Physics[END_REF]. In our case, ∆t is the sampling interval of the time series; in the discretized Sturm-Liouville problem, it is the mesh size. In both cases, the true operator being approximated is positive definite (or, more precisely, non-negative semi-definite), i.e., its eigenvalues are all positive (possibly one or more, but not all, being zero).

As in wavelet and Fourier analyses, time-dependent reconstructions of the entire signal {X i } or of a suitable filtered version thereof are possible from its PCs and EOFs. Hence, for a given set of eigenelements K, the corresponding reconstructed component (RC) is [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF] :

r K,i = 1 M i k∈K M j=1 a k,i-j ρ kj , (6)
where M i is a scaling factor equal to M when i is sufficiently far from the end points, and contains

Multi-scale SSA (MS-SSA)

Systematic comparisons between SSA, the wavelet transform and other spectral-analysis methods have been carried out in Ref. [START_REF] Yiou | Spectral analysis of climate data[END_REF][START_REF] Ghil | Advanced spectral analysis methods[END_REF] (see, in particular, Table 1 in [START_REF] Ghil | Advanced spectral analysis methods[END_REF]). Further analogies between certain mathematical features of SSA and wavelet analysis were mentioned in [START_REF] Yiou | Dynamique du Paléoclimat: Des Données et des Modèles[END_REF]. Table 1 here summarizes the most useful mathematical parallels between the two time-series analysis methods; these parallels provide the basis for extending global SSA analysis to a local one. In SSA, the largest scale at which the signal X is analyzed in Eq. ( 5) is approximately N (the length of the time series), and the largest period is M . As a consequence, the EOFs ρ k contains information from the whole time series, as in the Fourier transform.

Method SSA Wavelet transform Analyzing EOF ρ k Mother wavelet ψ function Basic facts ρ k eigenvectors of C M ψ chosen a priori Decomposition M t ′ =1 X(t + t ′ )ρ k (t ′ ) X(t)ψ t -b a dt Scale W = αM a Epoch t b Average ρ 1 ψ (0) & trend Derivative ρ 2 ψ (1)
Table 1: Analogy between SSA and the wavelet transform.

In order to define a local SSA, we propose to extend the SSA methodology by using a timefrequency analysis within a running time window whose size W is proportional to the order M of the covariance matrix. Varying M , and thus W in proportion, we obtain a multi-scale representation of the data. We perform local SSA on a time series by sliding windows of length W ≤ N , centered on times b = W/2, . . . , N -W/2. When using this method, we assume that considerable information content resides in the local variance structure and that the time series is locally the sum of a trend, statistically significant variability, and noise.

A priori, we can vary the two scales W and M independently, as long as W is larger than M , W/M ≥ α > 1, and α is large enough [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF]. In the wavelet transform, however, for instance the Mexican hat which corresponds to the second derivative of the Gaussian function, the number of oscillations of the mother wavelet is fixed and independent of the scale (width) of the analyzing wavelet. In this spirit, we fix the ratio W/M = α and rely therewith on the oscillation property of the EOFs described in Sec. 2.2 above to provide a fixed number of zeroes for the data-adaptive "wavelet" ρ k of each local SSA analysis; α = 3 in the calculations presented below and, as we shall see, it will in fact suffice to use k = 1 or k = 2 on each W -interval. This provides an analysis at a fixed scale W . Sampling a set of W values that follow a geometrical sequence, for instance in powers of 2 or 3, provides a multi-scale analysis very similar to the wavelet transform.

For a given position b and fixed W , we obtain local EOFs, for which PCs or RCs can be computed as described in Sec 2.2. The EOFs are the direct analogs of analyzing wavelet functions, as summarized in Table 1. Similarly to successive analyzing wavelets (for instance the derivatives of the Gaussian function), the number of EOF oscillations increases roughly with order and the zeroes of ρ k+1 separate those of ρ k . This is consistent with the orthogonality of distinct EOFs, with their being even or odd about the mid-point of the window, and with their tending in certain cases to pairs of sines and cosines [START_REF] Vautard | Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series[END_REF] or of Legendre polynomials [START_REF] Gibson | An analytic approach to practical space reconstruction[END_REF]. The first EOF thus corresponds approximately to an analyzing wavelet function with a single extremum and no zero, for instance the Gaussian. The second EOF is reminiscent of the first derivative of the Gaussian and so on.

Then for each b and each EOF ρ k , it is possible to obtain local PCs a k and RCs r k . The kth

PC at time b is a b ki = M j=1 X i+j-1 ρ b kj , (7)
and the corresponding RC is

r b ki = 1 M i M j=1 a b k,i-j ρ b kj , (8)
with b -W/2 ≤ i ≤ b + W/2.
The crucial difference between this local version and global SSA is that the RCs are obtained here from local lag-covariance matrices. As b varies from W/2 to N -W/2, this implies that the RCs will be truncated near the edges of the time series. We thus see that the local SSA method provides simultaneous "wavelet transforms" of the data by a set of analyzing wavelet functions, corresponding to the M different EOFs of the lag-covariance matrix. When W (and thus M ) is small, local SSA provides a small-scale analysis of the signal with a few distinct analyzing functions, namely a (sub)set of EOFs indexed by k. This is reasonable as there are not many possible structures at scales that approach the sampling time scale. On the other hand, at large scales, local SSA can also provide the simultaneous analysis by many different analyzing mother wavelet functions, {ρ k : 1 ≤ k ≤ M }, and thus reflect the large possible complexity of the structures that can develop over the entire time series.

The most important property of this local SSA analysis is that the analyzing functions (analogous to the wavelets) are data adaptive: they are just the EOFs of the lag-covariance matrix with time lags up to M , within a time window of size W = αM . In other words, the shape of these analyzing functions is not imposed a priori, like in a wavelet analysis, but explicitly depends on the time series itself. This property is particularly suitable in analyzing time series that exhibit dominant structures which differ along the signal. For instance, an oscillatory behavior could be followed by white or colored noise and then by deterministically intermittent behavior; this could indicate regime transitions that the system which generates the signal underwent while under observation. If so, one would like to have an analyzing wavelet which is adapted to each section of the signal. Our data-adaptive scheme will definitely help follow such regime transitions in time.

In addition, there is information to be gained from the shape of the EOFs that are provided by the data. Wavelet transforms are sometimes used to test for self-similarity: one tests for the existence of a power-law dependence of a wavelet coefficient at a given time as a function of scale. If the power law is observed, one concludes on the existence of a local singularity associated with such a law. This is useful in analyzing multifractal structures characterized by a set of distinct singularities [START_REF] Arneodo | Beyond classical multifractal analysis using wavelets: Uncovering a multiplicative process hidden in the geometrical complexity of diffusion limited aggregates[END_REF].

The approach presented here provides an additional test. Indeed, if a signal is locally self-similar, the shape of the EOFs must also be the same, as we change scale (i.e., M and W ). We can thus carry out an analysis at different scales, obtain the data-adaptive analyzing functions and test for their self-similarity. This possibility arises naturally out of our multi-scale SSA (MS-SSA) method and does not exist in any technique that assumes a priori the shape of the analyzing function, however carefully selected. We present a few tests of MS-SSA below, and compare the results with those of a standard wavelet analysis.

The interesting properties of MS-SSA come however at a significant computational cost. Indeed, the numerical cost for matrix diagonalization far exceeds that of fast Fourier transforms or simple algorithms for wavelet analysis, even when several wavelet functions are used. In order to circumvent this problem and reduce the computational cost, we carry out the local SSA analysis only for a subset of time steps, b = β∆t. This sampling interval for the analysis must be substantially smaller than the scale parameter W , i.e. β ≪ W . Then we use Eq. ( 8) for b = β∆t to obtain the local RCs where the MS-SSA algorithm is not applied. We have checked that this interpolation procedure is very accurate when compared to exhaustive local SSA analysis at each time step, i.e. computing r b kb in Eq. ( 8) for each b and hence discarding its values at epochs surrounding b. This procedure is based on the relative robustness of the EOFs, which is also used in SSA-based time-series prediction [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF][START_REF] Keppenne | Adaptive filtering and prediction of the Southern Oscillation index[END_REF][START_REF] Ghil | Recent forecast skill for the El Nio/Southern Oscillation[END_REF].

3 Data sets

Synthetic data sets

We would like to investigate the behavior of MS-SSA on time series that are self-similar or fractal, as they pose a priori the biggest challenge to techniques based on well-behaved finite covariances.

As already mentioned, such features are ubiquitous in the geosciences [START_REF] Dubrulle | Scale Invariance and Beyond[END_REF].

The first time series P1 that we shall investigate was obtained by using an Iterated Function System (IFS) [START_REF] Barnsley | Fractals Everywhere[END_REF] to approximate a simple triadic set. The procedure that we have used is described in Appendix 7.1 and the resulting time-series is shown in Fig. 1. A closely related time series P 1 is obtained from the characteristics of the exact triadic Cantor set. Its construction is described in Appendix 7.2. In the "perfect" construction scheme of P 1 , we anticipate large "resonances" in the lag-covariance matrix that should be smoothed out by the noise present in the construction of P1 . Finally, a multiple-rule Cantor set is generated in Appendix 7.2 to yield the time-series P1 . The second process is a log-periodic signal with a square-root singularity. In the context of critical phenomena, the log-periodicity corresponds to the existence of an imaginary part in the critical exponent, and is associated with a discrete scale invariance [START_REF] Sornette | Discrete scale invariance and complex dimensions[END_REF], i.e. to the invariance of the system or of its properties only under magnifications that are integer powers of a fundamental ratio. Initially, complex exponents were proposed as formal solutions of renormalization-group equations in the seventies [START_REF] Jona-Lasinio | The renormalization group: a probabilistic view[END_REF][START_REF] Nauenberg | Scaling representations for critical phenomena[END_REF][START_REF] Th | Phase Transitions and Critical Phenomena[END_REF]. In the eighties, they have been shown to emerge in various physical problems that arise in discrete hierarchical systems [START_REF] Derrida | Fractal structure of zeros in hierarchical models[END_REF][START_REF] Bessis | Mellin transforms of correlation integrals and generalized dimension of strange sets[END_REF].

Recently, it has been realized that discrete scale invariance and its associated complex exponents and log-periodicity may appear "spontaneously," without the need for a pre-existing hierarchy, as a result of dynamical processes. Such behavior has been found in models and experimental data from irreversible growth processes, rupture in heterogeneous systems, earthquakes, percolation and even in financial crashes [START_REF] Sornette | Large financial crashes[END_REF] (see [START_REF] Sornette | Discrete scale invariance and complex dimensions[END_REF] for a review, and references therein). Another source of time series with log-periodic behavior is provided by Boolean delay equations, whose solutions are Booleanvalued processes in continuous time [START_REF] Dee | Boolean difference equations, I: Formulation and dynamic behavior[END_REF]. These equations have been used to model paleoclimatic variability [START_REF] Darby | A Boolean delay equation model of an interdecadal arctic climate cycle[END_REF][START_REF] Ghil | Deep water formation and Quaternary glaciations[END_REF], interdecadal [START_REF] Wright | A note on Quaternary climate modelling using Boolean delay equations[END_REF] and most recently, seasonal-to-interannual [START_REF] Saunders | A Boolean Delay Equation Model of ENSO Variability[END_REF] climate change, although not all of these applications exhibit the behavior in question. In the case of Boolean delay equations, the log-periodicity exhibits a period and amplitude that increase as t → t c = +∞ [START_REF] Dee | Boolean difference equations, I: Formulation and dynamic behavior[END_REF] and their treatment by MS-SSA is left for subsequent study.

The log-periodicity we concern ourselves with here is reflected in oscillations of decreasing period and amplitude that are superimposed on a power-law behavior and culminate at a singularity t c in finite time. At the critical point t c , the instantaneous frequency of the oscillation becomes infinite, but its amplitude vanishes. Again, a global spectral analysis method would fail to represent this signal in a satisfactory manner since its nonstationarity is characterized by a continuously varying period and a nonlinear trend. The log-periodic time series P 2 is constructed according to the procedure given in Appendix 7.3 and is shown in Fig. 2. While the first two processes (P 1 , P1 ) and P 2 are purely deterministic, the third process we have studied still presents self-similar properties, albeit with a large stochastic component. Its self-similarity stems from its multiplicative-noise nature. The process P 3 has the property of being statistically stable under affine transformations [START_REF] Sornette | Convergent multiplicative processes repelled from zero: power laws and truncated power laws[END_REF][START_REF] Sornette | Linear stochastic dynamics with nonlinear fractal properties[END_REF]. A 512-point realization of P 3i , for 1 ≤ i ≤ 512, is shown in Fig. 3 and a histogram is provided to illustrate the power-law behavior and its long tail.

Figure 3: Multiplicative-noise time series P 3 .

Climatic data set

In addition to the synthetic data sets that possess well-understood properties, we chose to analyze a "real-world" climatic time series, the Southern Oscillation Index (SOI). This time series illustrates well another type of behavior often encountered in the geosciences, viz., the superposition of quasiperiodic and aperiodic behavior [START_REF] Ghil | Recent forecast skill for the El Nio/Southern Oscillation[END_REF][START_REF] Ghil | Nonlinear dynamics and predictability in the atmospheric sciences[END_REF]. SOI is a climatic index connected with the recurring El Niño/Southern Oscillation (ENSO) conditions in the tropical Pacific; it is essentially the monthly mean difference in sea-level pressure between Darwin, Australia, and Tahiti (Fig. 4). An anomalously negative value of this index indicates a warm ENSO event (El Niño), while a highly positive value is associated with a cold event (La Niña). SOI has been used in numerous studies for understanding the dynamics of ENSO [START_REF] Jiang | Quasi-quadrennial and quasi-biennial variability in the equatorial Pacific[END_REF] and improving its prediction [START_REF] Keppenne | Adaptive filtering and prediction of the Southern Oscillation index[END_REF][START_REF] Ghil | Recent forecast skill for the El Nio/Southern Oscillation[END_REF].

In the data set we use here, the annual cycle was removed and the time series was normalized by its variance. The time interval considered goes from January 1933 to December 1996, during which very few observations are missing at either station. ENSO models have shown that, in the Tropical Pacific, phase locking of the coupled oceanatmosphere system's self-sustained oscillation to the annual cycle leads to a Devil's staircase [START_REF] Chang | Interactions between the seasonal cycle and the Southern Oscillation: Frequency entrainment and chaos in an intermediate coupled ocean-atmosphere model[END_REF]. This behavior has been confirmed by studies of sea-surface temperature data sets using both classical (global) SSA and (fixed-analyzing function) wavelet methods. Such studies have been performed on the sea-level pressure record at Darwin only [START_REF] Wang | Temporal structure of the Southern Oscillation as revealed by waveform and wavelet analysis[END_REF], which spans a longer time interval, using wavelet and waveform analysis, and Tropical Pacific sea-surface temperatures with global SSA [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF][START_REF] Rasmusson | The biennial component of ENSO variability[END_REF][START_REF] Jiang | Quasi-quadrennial and quasi-biennial variability in the equatorial Pacific[END_REF]. Two of these analyses have shown a fairly sudden shift in the characteristic period of ENSO's lowfrequency component that occurred in the 1960s, with a longer period before and a shorter one prevailing since [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF][START_REF] Wang | Temporal structure of the Southern Oscillation as revealed by waveform and wavelet analysis[END_REF].

Numerical results

Cantor-set experiments

SSA analysis

Global methods of time-series analysis will clearly miss the multi-scale feature of the IFS construct P1 (t), and essentially see one or more nearly periodic components that arise from a smoothed version of the exact Cantor set P 1 (t). A Monte-Carlo SSA [START_REF] Ghil | Spectral methods: What they can and cannot do for climatic time series[END_REF][START_REF] Allen | Monte Carlo SSA: Detecting irregular oscillations in the presence of coloured noise[END_REF] of P1 (t) with a lag of M = 40 shows two frequencies that appear as pairs of eigenvalues (0.12 and 0.25 cycles per time unit) emerging high above a red-noise-like background; two additional eigenvalues stand slightly above the red-noise background at frequencies 0.006 and 0.06 cycles per time unit (Fig. 5). The pair at frequency 0.25 (period 4) represents the small-scale oscillations in P1 ; the smallest possible scale at this step of the IFS is 2, but it is not systematically reached throughout the interval of interest, so that the smallest scale which is always present is 4 (see Fig. 1b). The pair at frequency 0.12 cycles per unit (0.12 cpu; period 8.33) accounts for the grouping of two smaller-scale oscillations. Finally, the pair near 0.09 cpu (period 11.1; not shown) appears to be closely associated with the 0.12 cpu pair; hence the latter two pairs should be viewed as a quadruplet. The first eigenvalue has an associated period (156) larger than the lag window M = 40, so that it has no statistical significance; visually, it corresponds to roughly half the length of the largest-scale alternation between quiescence and oscillations (1/6 of the entire interval). SSA enhances one characteristic scale (which depends on M ) and smears out the faster ones. On the other hand, the RCs of the low-order EOFs exhibit mixtures of fast and slow oscillations. This is so because all the oscillations have nearly the same variance by construction (Fig. 6): the so the variance only depends on the distribution of the oscillations over the interval. The fast oscillations can thus be viewed as spurious because they do not provide information on the way the data set was constructed. Other global methods (Maximum Entropy: [START_REF] Yiou | Spectral analysis of climate data[END_REF][START_REF] Ghil | Spectral methods: What they can and cannot do for climatic time series[END_REF][START_REF] Ghil | Advanced spectral analysis methods[END_REF][START_REF]Modern Spectrum Analysis[END_REF]; Multi-Taper: [START_REF] Yiou | Spectral analysis of climate data[END_REF][START_REF] Ghil | Spectral methods: What they can and cannot do for climatic time series[END_REF][START_REF] Ghil | Advanced spectral analysis methods[END_REF][START_REF] Thomson | Spectrum estimation and harmonic analysis[END_REF]) do not provide more relevant information on P1 (t) either (not shown).

Wavelet transforms

We next performed wavelet transforms of P1 with the Gaussian analyzing function ψ(x) = exp(-x 2 /2) and its derivative ψ ′ . Geometric scale increments a = 3 n , n = 0, . . . , 6, were used. The choice of an increment of 3 is suggested by the triadic structure of the Cantor set; other choices (e.g., a = 2 k ) clearly cannot provide the same optimal scale decomposition (see below). More generally, when the data set contains a discrete scale invariance with respect to dilations, the use of the preferred scaling ratio -such as the ratio 3 in the Cantor set studied here -is optimal [START_REF] Arneodo | Beyond classical multifractal analysis using wavelets: Uncovering a multiplicative process hidden in the geometrical complexity of diffusion limited aggregates[END_REF][START_REF] Dubrulle | Scale Invariance and Beyond[END_REF]. The wavelet analyses of P1 are shown in Fig. 7. The triadic and self-similar structure of the Cantor set is clarified, as each scale picks up one step of the construction process. A wavelet analysis with ψ ′ as the mother wavelet (right column of Fig. 7), rather than ψ (left column), reveals further aspects of the process's singularities. The triadic construction of P1 is also apparent in this decomposition, albeit not as clearly as with ψ. The wavelet transforms of P1 in Fig. 7 do not differ significantly from those of its "exact" counterpart P 1 (not shown).

MS-SSA analysis

• RC analysis. We performed MS-SSA analyses of P1 with scales W varying with geometric increments, W = 2 • 3 n , and a ratio W/M = 3. We kept the first two EOFs for each scale. The first RCs of P1 reveal very well the triadic structure of the Cantor set (left column of Fig. 8), as in the Gaussian-wavelet analysis. Unlike the wavelet analysis, MS-SSA restricted to EOF-1 reproduces almost exactly the steps taken to generate P1 , i.e., the divisions do not overlap, even though parts of the reconstruction are lost due to side effects. Thus, in the case of P1 , the triadic structure seems to be more faithfully captured by this method. A slight drawback, on the other hand, is that SSA exhibits a mild "Gibbs effect," i.e., an overshoot at certain discontinuities; it does not guarantee, therewith, the positiveness of the signal reconstruction which is ensured by a wavelet analysis with a positive function, like the Gaussian (but not its derivative: see right column of Fig. 7).

The RCs corresponding to EOF-2 (right column of Fig. 8) provide an analysis of the Cantor set's "derivative," as with the ψ ′ -wavelet transform. Hence, EOF-2 is analogous to ψ ′ (see Table 1 and ref. [START_REF] Gibson | An analytic approach to practical space reconstruction[END_REF]), and the resulting RCs yield an analysis of the derivative of P1 (not shown).

• EOF analysis. As anticipated in Section 2, we expect the EOFs ρ (W ) k to be similar for fixed order k and on different scales W , due to the choice of the scaling ratio 3 being equal to that used in constructing the time series. Figure 9 shows EOFs 1 through 4 of P1 , for scales W in a fixed ratio of 3 between 9 and 243 (a scale of 3 implies M = 1, for which the analysis is trivial), and at a location b = 333, situated at one third of the time series. The EOFs are scaled by 3 n/2 , so that their amplitudes match. We notice a near-convergence in the shape of EOFs 1-3 as the scale increases; this indicates that the time series is self-similar with scale increments of 3. EOF 4, on the other In contrast, we show in Fig. 10 the same EOFs 1-4 of P1 , for scales W in a fixed ratio of 2, between 16 and 256, at the same position b = 333 along the time series. There is practically no convergence with W for fixed k, reflecting the mismatch between the scale factor 2 of the multiscale analysis and the preferred scale factor 3 of the time series. This suggests a novel method for detecting the existence of a preferred scaling factor λ in a time series by optimizing the scale ratio l in MS-SSA: only when l ≈ λ m , where m is an integer, will the convergence of the EOFs ρ This phenomenon reflects a "scale-locking" between the analyzing (multi)scale factor and the intrinsic scale ratio(s) present in the time series. This provides a scaling test that conventional wavelet analysis cannot offer, due to the a priori choice of ψ. In MS-SSA, the functional-shape constraint is relaxed and new information on the signal can thus be obtained.

(W ) k → ρ (∞) k be good as W → ∞.
• Resonance. An interesting phenomenon appears when applying MS-SSA to the "perfect" Cantor set P 1 (t). The local lag-covariance matrices seem to exhibit resonances between scales so that a given "slow" scale contains information on faster scales and hence exhibits oscillations that blur the triadic decomposition (not shown). This was not seen in the MS-SSA of the "approximate" time series P1 (t), which contains noise in both time and scale that smoothes out such resonances and leads to a clear reconstruction of the triadic structure. The presence of a small amount of noise often increases the robustness, and hence ease of identification, of a phenomenon's main features. Wellknown instances are the randomization method in general probability theory [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF] and the addition of process or observation noise in sequential estimation theory [START_REF] Ghil | Data assimilation in meteorology and oceanography[END_REF].

• Change of geometrical structure along the scale axis. MS-SSA analysis of the multiple-rule Cantor set P1 (Fig. 11) shows that the method can capture the difference in structure beyond a threshold scale, even though the transition may be smoothed out by the two construction procedures. The RCs do not give as clearcut a triadic reconstruction of this Cantor set as for the single-rule set P1 in Fig. 8. But, in spite of the relatively small data set we used (only 729 points), we can still detect in Fig. 11 the difference between the use of rule [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF] for the larger scales and rule [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF] for the smaller ones at the passage between the scales W = 9 and W = 27. • Reconstruction process by summing over scales for a single EOF. By construction, there is a finite number M of EOFs for each scale W in MS-SSA. The EOFs form an orthonormal basis and hence the local decomposition of the time series occurs on a finite number M of modes, at each scale. Thus, for a given scale W , the sum of the corresponding RCs reproduces the original time series in each window. In wavelet language, this corresponds to reconstructing a signal by summing, at a fixed scale a (see Table 1), over a finite number of wavelet transforms; each of these transforms uses a different mother wavelet -for instance, the successive derivatives ψ (n) , 0 ≤ n ≤ M -1, of the same mother wavelet ψ = ψ (0) -such that the finite set of mother wavelets form a (not necessarily orthogonal, but still nondegenerate) basis. This, however, does not correspond to the usual reconstruction formula (3).

In continuous wavelet analysis, the reconstruction is usually performed by summing over the scales for a fixed mother wavelet, according to Eq. ( 3), often called a "resolution of the identity (operator)" on admissible functions [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]. For discrete wavelet transforms, there does not exist, in discrete wavelet basis is not orthogonal in general and special iterative methods must be developed to derive an inversion formula (see Chapter 3 in [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]).

Likewise, for our MS-SSA method, the expansion over discrete scales of a function with a fixed EOF order is not orthogonal in general. The MS-SSA approach is thus more closely related to the theory of "frames" [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]. It is easy to show from the definition of a frame that the family of EOFs ρ W k of given order k for discrete scales W n = a n with n integer constitutes such a frame. This yields an algorithm to reconstruct the signal from the discrete set of wavelet transforms, in our case the RCs of same order k at different scales, for special choices of the scale W and translation parameter b [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF] :

X i = m M j=1 ρWm k,j a Wm k,i-j ; (9)
here the sum over m corresponds to the different scales, and a Wm k,i-j is given by

a Wm k,i-j = M l=1 ρ Wm k,l X i-j+l-1 , (10)
while ρWm k,j is a vector related to the EOF ρ Wm k,j through an iterative algorithm [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]. The inversion can thus be formulated in principle, but in practice it is quite laborious. Emphasizing the practical aspect, we propose an approximate inversion formula similar to (9) but with ρWm k,j replaced by ρ Wm k,j :

X i = m M j=1 1 M ρ Wm k,j a Wm k,i-j . ( 11
)
This formula is tested in Fig. 12. The top panel shows the time series P 1 . The second panel shows the RC using EOF-2 (k = 2) for the single scale W = 9. This corresponds to formula [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF] with a single value of k = 2, while a = 3 and thus W 2 = 9. The third, fourth and fifth panels show the sum [START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF] with two, three and four terms; these terms correspond to W 2 = 9, W 3 = 27, W 4 = 81 and W 5 = 243, while the order of the EOF ρ Wm k is always k = 2.

Figure 12: Reconstruction of P 1 using Eq. ( 11) for scale increments of 3.

A log-periodic process

In this subsection, we focus on the detection and quantification of self-similar structures in a logperiodic process. In the case of Eq. ( 17), we expect discrete self-similarity near t = t c , with a preferred scaling of λ = 2. Indeed, in Fig. 2 the distance between successive maxima of the oscillations decreases geometrically with a scaling of about 2, and thus converges to zero as t approaches t c .

In the case of P 2 , a straightforward application of MS-SSA (with EOFs 1 and 2) and of the wavelet transform (with ψ and ψ ′ ) give extremely similar results; they both act as low-pass filters with geometrically increasing width (not shown) and hence do not provide a satisfactory reconstruction of the signal. We analyzed therefore P 2 (t) by MS-SSA with scaling increments of 2 but placing ourselves near t c = 100 in the sense that the local windows always end at t c , and their midpoint moves to the left. The first four EOFs, shown in Fig. 13, clearly converge rapidly as the size of the windowand thus the scale analyzed -increases. If we choose another ratio, say 3, to scale the windows no such convergence is observed (not shown), as was the case in comparing Fig. 9 with Fig. 10 for the Cantor set. This sensitivity of MS-SSA's EOFs to the correct scaling factor suggests a new way to characterize log-periodicity.

Figure 13: EOFs 1-4 of P 2 near t c , with 2 n scaling.

EOFs 1 and 3 capture the local mean of the signal and scale the period of its oscillation, respectively, having no zero and two zeroes (see Fig. 13). RC-2 describes best the local fast oscillations near the critical point t = t c , while RC-1 at different scales captures well the overall behavior of the time series in Fig. 2 (see Fig. 14).

Figure 14: RCs 1 and 2 of P 2 for scale increments of 2.

We look next at what happens near t = 0, i.e., we use windows starting at the beginning of the time series. The convergence of EOFs 2-4 is still excellent as W increases (not shown). Only EOF-1, which is sensitive to the overall behavior of the time series, exhibits slower convergence with W , due to the lack of exact self-similarity of P 2 (t) near the origin. This robust, albeit imperfect, convergence of its EOFs for log-periodic time series makes the MS-SSA approach potentially useful for the analysis of real log-periodic data, in the presence of noise and truncated singularities. The application of this technique to real data on earthquakes [START_REF] Sornette | Complex critical exponents from renormalization group theory of earthquakes : Implications for earthquake predictions[END_REF][START_REF] Johansen | Discrete scaling in earthquake precursory phenomena : evidence in the Kobe earthquake, Japan[END_REF] and financial [START_REF] Sornette | Large financial crashes[END_REF] time series will be investigated separately.

Multiplicative noise

As expected, classical global methods which assume the existence of finite variance do not give a particularly good insight into the process P 3 (t). Monte-Carlo SSA yields a periodicity that is statistically significant above 95%; the actual period, though, lies between 3 and 4 cycles and fluctuates from realization to realization (not shown). The time series cannot be otherwise distinguished from red noise; this is obviously wrong and can be explained by the lack of stability in the global variance estimate.

The global RCs (Fig. 15) provide a sensible scale-dependent reconstruction for the slower scales: the characteristic scale halves from RC-1 to RCs 2-3, on to RC-4 and, again to RCs 5-8, but the reconstructions become noisier as the order increases further (not shown). The total reconstruction of P 3 (t) from the eight leading RCs (Fig. 15e) is an excellent coarse-scale representation of P 3 (t), as depicted in Fig. 3. Next, wavelet and MS-SSA analyses were both performed with geometric scale increments of 2, while the latter used a ratio W/M = 3 as in Secs. 4.1 and 4.2. A wavelet analysis using the Gaussian mother wavelet (Fig. 16, left column) progressively smoothes the details of the time series as the scale parameter a increases. The analysis reveals the natural self-similarity of P 3 , as the analysis at scale a = 32 can be normalized to match the first half of the analysis at scale a = 8 (Fig. 16). In turn, the analysis at scale a = 8 can be normalized to fit the analysis at scale a = 2 in the same way. Thus, we can infer that scale increments of about 4 are sufficient for the analysis of P 3 .

Figure 16: Wavelet analysis of P 3 . MS-SSA using EOF-1 (Fig. 17, left column) behaves in a similar way, with a progressive smoothing of the time series and the same scaling result. RC-2 captures the oscillating variability (right column of Fig. 17) around EOF-1, as does ψ ′ with respect to ψ (right column of Fig. 16). In this case, the parallel between the behavior of wavelet analysis and MS-SSA is striking. As we saw in Fig. 15e, the reconstruction of P 3 can also be completed with a finite sum of EOFs, for a fixed scale of analysis, whether global of not.

In contrast to the results found in Fig. 9 for the Cantor set P1 , while the EOFs of the same order have a similar shape across scales, they cannot be exactly superimposed (not shown). This simply reflects the fact that the time series P 3 is not exactly self-similar with a preferred scaling ratio of 2. In fact, no choice of a single scaling ratio provides a good convergence of the EOFs for fixed k: P 3 is not self-similar but rather exhibits a more complex scaling structure, that is approximately wavelets. The data-adaptive virtues of the EOFs are illustrated in Fig. 18, where the shape of the local EOF-1 follows closely the variations of P 3 .

Figure 18: The evolution in time of the local EOF-1 for P 3 .

On this data set, the two methods under consideration perform in a rather similar way. The advantage of MS-SSA is that it fully uses the information contained in the time series in order to compute analyzing functions, and that the reconstruction is intrinsically complete, given a small set of EOFs, which is not the case for wavelets.

Southern Oscillation Index

We performed MS-SSA on the monthly SOI data for the years 1933-1996 (see Section 3.2). The parameters were, as before, W/M = 3 and geometric scale increments of 2. In addition to the previous analyses, we also computed an "instantaneous" frequency for each local EOF. This was simply done by least-square fitting a sine wave to each local EOF of interest. The instantaneous frequency can also be obtained from a complex wavelet transform [START_REF] Farge | Wavelet transforms and their applications to turbulence[END_REF], by using information on the phase of the transform.

Our analysis does not reveal any evidence of self-similarity or fractality, as obtained for the previous synthetic data sets. Instead, we find a preferred scale of variability between 3 and 5 years (not shown), which corresponds to ENSO's low-frequency mode [START_REF] Rasmusson | The biennial component of ENSO variability[END_REF][START_REF] Jiang | Quasi-quadrennial and quasi-biennial variability in the equatorial Pacific[END_REF]? The first two local EOFs are consistently paired and in phase quadrature, which shows that the nonlinear oscillation associated with this mode is robust and persists throughout the sixty odd years being examined.

The computation of the instantaneous frequency allows us to detect an abrupt frequency shift of the ENSO mode near 1960 (Fig. 19). The characteristic periodicity goes from 57 months (between 1943 and 1961) to 39 months (between 1963 and 1980). A comparable decrease in period was observed in the early 1960s by Moron et al. [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF] in Tropical Pacific sea-surface temperatures, by using multi-channel global SSA, and by Wang and Wang [START_REF] Wang | Temporal structure of the Southern Oscillation as revealed by waveform and wavelet analysis[END_REF] in a sea-level pressure record at Darwin, using a wavelet transform. Our MS-SSA approach improves on these previous results in two respects. First, compared to both the global SSA analysis [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF] and the wavelet and waveform (matching pursuit) approach [START_REF] Wang | Temporal structure of the Southern Oscillation as revealed by waveform and wavelet analysis[END_REF], the combination of our local SSA in a sliding window with the local sine-wave fit gives a much sharper signature of the increase in instantaneous frequency. Second, we can check this regime change at different scales by performing the same analysis separately at each scale. In contrast, the wavelet and waveform approaches need all the scales together in order to be able to detect the frequency change. We thus have at our disposal an additional degree of freedom to test for the reliability and robustness of the signal.

We find that, for W = 64 and 32, the transition is not as sharp as with W = 128 (the value used in Fig. 19), but is still identifiable. This degradation is expected since the maximum period M = W/3 that can in principle be retrieved at these values of W is less that the ENSO period of interest. This means in practice that only scales of the order of W = 128 or larger will capture the ENSO period in question and the change in it. At scales that are even larger, however, the analysis becomes limited in its reliability by the finite length of the time series. Hence the value W = 128 provides a good compromise between the two conflicting constraints of resolution and statistical significance [START_REF] Vautard | Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series[END_REF][START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF][START_REF] Ghil | Spectral methods: What they can and cannot do for climatic time series[END_REF][START_REF] Ghil | Advanced spectral analysis methods[END_REF].

We have presented a multi-scale extension of the orthogonal component decomposition of a time series known as singular spectrum analysis (SSA). This multi-scale extension, which we called MS-SSA, has been performed by varying the width W of the analyzing window, by analogy with wavelet analysis. The successive analyses used a fixed ratio -equal to 2 or 3 in the numerical examplesbetween the width W and the order M of the local covariance matrix.

Our multi-scale SSA (MS-SSA), or data-adaptive wavelet analysis, extends various other approaches in the same spirit. Mallat et al. [START_REF] Mallat | Adaptive covariance estimation of locally stationary processes[END_REF] describe a hierarchical binary tree of basis functions for locally stationary time series. Lilly and Park [START_REF] Lilly | Multiwavelet spectral and polarization analyses of seismic records[END_REF] select their Slepian wavelets [START_REF] Slepian | Prolate spheroidal wave functions, Fourier analysis and uncertainty-V: The discrete case[END_REF] at a given scale by optimizing the spectral energy in a given frequency band [START_REF] Thomson | Spectrum estimation and harmonic analysis[END_REF] and vary the scale. Mann and Park [START_REF] Park | Interannual temperature events and shifts in global temperature: A "multiwavelet" correlation approach[END_REF] have applied this approach to a multivariate analysis of interannual temperature variations. Wavelet packets [START_REF] Wickerhauser | Lectures on Wavelet Packet Algorithms[END_REF][START_REF] Learned | A wavelet packet approach to transient signal classification[END_REF] can be defined as smooth versions of the Walsh system used in binary signal processing. While they extend wavelet analysis and take a step in the direction of satisfying our requirements (A) and (B) above, the wavelet-packet approach does not allow the same degree of flexibility and data-adaptivity as provided by our method.

Coifman and Saito [START_REF] Coifman | The local Karhunen-Loève bases[END_REF] use the Karhunen-Loève expansion that underlies SSA to construct a dictionary of orthonormal bases that has a tree structure. Their starting point is thus very similar to ours; Coifman and Saito, however, use additional criteria to choose the orthonormal bases from the overdetermined dictionary in order to perform a multi-scale analysis. Our Karhunen-Loève decomposition is performed in moving windows of size W that maintain a fixed ratio to the order M of the decomposition. As a consequence, our orthogonal basis is unique, minimal and offers a very intuitive interpretation of the decomposition.

The empirical orthogonal functions (EOFs) play the same role in MS-SSA as the "mother wavelets." In contrast to standard wavelet transforms, our basis functions are data-adaptive, i.e. they are determined from a diagonalization of the covariance matrix C M of the data within the observing window. The data-adaptive character of the MS-SSA basis allows us to test further for self-similar properties of a time series by comparing the EOFs at different scales.

In MS-SSA's orthogonal decomposition, the sum of a finite set of coefficients weighted by the corresponding "wavelets," i.e., the reconstructed components (RCs) at the same scale, reconstructs exactly the initial signal. This is in contrast to the standard wavelet decomposition in which one needs to sum over all scales or sum over an infinite set of wavelets (for instance corresponding to the successive derivatives of a Gaussian) at a given scale. We also propose an approximate "resolution of identity" inversion formula to reconstruct the signal from a nonorthogonal sum over scales of RCs that correspond to a fixed EOF order. This formula is the counterpart of the usual wavelet inversion formula and arises from the theory of "frames" [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF].

We tested MS-SSA on a few irregular time series for which conventional single-scale spectral analysis techniques would fail to reveal their properties: Cantor sets and intermittently amplifying multiplicative noise. In many respects, wavelets and MS-SSA behave in similar ways : MS-SSA decomposition onto the first and second EOFs is very similar to a wavelet analysis that uses the Gaussian mother wavelet and its first derivative. But in the case of the approximate Cantor set P1 , the MS-SSA decomposition gives a clearer insight on the way the set was generated. In addition, the shape of the EOFs also allowed us to check precisely for the self-similarity of that set.

The latter property is probably the most interesting one to investigate in the future, as the change in shape of the EOFs as a function of scales provide a new method to test for discrete scale invariance. Indeed, the detection of the preferred scaling ratio in a discretely scale-invariant system may be obtained by finding the ratio of the multi-scale analysis that minimizes the difference between the data-adaptive EOF "wavelets" obtained at different scales. This test appears to be more robust than previously used methods [START_REF] Sornette | Discrete scale invariance and complex dimensions[END_REF] and has been demonstrated successfully on a time series that presents the scaling structure of the triadic Cantor set.

Finally, we applied the MS-SSA methodology to a real climatic time series, given by the monthly Southern Oscillation index (SOI) for 1933-1996. The SOI is a broadly used indicator for the El Niño/Southern Oscillation (ENSO). The sharp transition in ENSO's low-frequency mode from a period of roughly 5 years to roughly 3 years in the early 1960s refines earlier results obtained using global SSA [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF] and standard wavelet [START_REF] Wang | Temporal structure of the Southern Oscillation as revealed by waveform and wavelet analysis[END_REF] analysis. It supports the Devil's staircase theory of ENSO that an interdecadal change in mean thermocline depth in the Tropical Pacific [START_REF] Ghil | Recent forecast skill for the El Nio/Southern Oscillation[END_REF][START_REF] Knutson | Simulated ENSO in a global coupled ocean-atmosphere model: Multidecadal amplitude modulation and CO 2 sensitivity[END_REF] causes the coupled ocean-atmosphere system there to "jump" from one broad step of the staircase to another. We hope that similar physical inferences will be facilitated in other areas of the geosciences to which MS-SSA might be applied. The time series P1 was obtained by using an Iterated Function System (IFS) [START_REF] Barnsley | Fractals Everywhere[END_REF] to approximate a simple Cantor set. An IFS is a contracting stochastic map constructed to converge towards a self-similar attractor of zero measure. The example we study consists in a time series whose characteristics mimic the triadic Cantor set. To generate this time series, we use the following one-dimensional IFS:

Acknowledgments

x n+1 = 1 3 x n with probability 1 2 , 1 3 x n + 2 3
with probability 1 2 .

(

) 12
Starting from an arbitrary point x 0 ∈ [0, 1], it is straightforward to verify that this stochastic IFS map converges to a triadic Cantor set on [0, 1], characterized by the fractal dimension of log 2/ log 3 [START_REF] Barnsley | Fractals Everywhere[END_REF]. We iterated Eq. ( 12) 10000 times and discarded the first 1000 values of the set so generated in order to obtain a set {x i } that has converged close to the Cantor set target.

We sorted these values, x i ∈ [0, 1], into increasing order, 0 ≤ x 0 ≤ ... ≤ x i ≤ x i+1 ≤ ... ≤ x N ≤ 1. A time series P1 with, say N = 1000 equidistant points, is generated by interpolating the set {x i } every ∆t = 1/N and scaling it to N :

P1j = 1 if ∃x i ∈ [j/N, (j + 1)/N ], 0 otherwise, (13)
with 0 ≤ j ≤ N . As a result of the way this time series is generated with the stochastic IFS, it contains residual noise with respect to a perfect Cantor set both in time and scale, because of the interpolation process (Fig. 1). This means that the smallest scales are not always of unit length, due to incomplete interval division.

The exact triadic Cantor set

The second time series P 1 was obtained recursively by iterating the binary map:

1 → 101, ( 14
) 0 → 000 , (15)
on the initial list s 0 = (1). At the L-th iteration, the list s L contains 3 L bits which have an exact Cantor set structure. We generated the time series P 1 with L = 6 iterations, so that P 1j = s L (j), 1 ≤ j ≤ 3 L . A "noisy" Cantor time series could be obtained by assigning a probability of occurence between rule ( 14) and 1 → 111. This procedure emulates the incomplete division scheme towards smaller scales that we observe in P1 .

A large variety of Cantor sets can be obtained by modifying the rules [START_REF] Moron | Trends, interdecadal and interannual oscillations in global sea-surface temperatures[END_REF][START_REF] Meyer | Ondelettes et Opérateurs I: Ondelettes[END_REF] in a straightforward manner. For instance, we generated a multi-rule Cantor set by iterating rules (14, 15) three times, then used the rules 1 → 110, [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF] 0 → 000.

three more times, in order to obtain a set P1 of 729 (= 3 6 ) points.

The fractal dimension of a set generated after iterating infinitely many arbitrary permutations of rules ( 16) and ( 14) is still log 2/ log 3, but the geometry of the resulting set obviously depends on the sequence in which the rules were applied. This is a clear instance where the fractal dimension is not enough to characterize the geometric structure [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF][START_REF] Blumenfeld | Lévy dusts, Mittag-Leffler statistics, mass fractal lacunarity, and perceived dimension[END_REF][START_REF] Solis | Lacunarity of random fractals[END_REF]. Other rules and combinations of rules can be used to obtain other fractal dimensions and, for a given dimension, distinct geometric structures of the resulting set.

Log-periodic time series

The functional form of the log-periodic process P 2 is simply:

P 2 (t) = A + B(t -t c ) α {1 + C cos [2π log(t -t c )/ log(λ)]} . (17)
Here, we take A = 2, B = -1, C = 0.2, α = 1/2 and λ = 2 (see Fig. 2). The critical time is t c = 100; we regularly sample 1000 time steps between 1 and t c with a time interval ∆t = 0.1, and define the P 2 time series by P 2i = P 2 (i∆t), 1 ≤ i ≤ 1000. Note that in our case α = 1/2 and thus lim t→tc P 2 (t) = A; hence P 2 is indeed continuous on the [0, t c ] interval, although it is not differentiable at t = t c .

7.4

The statistically self-similar process P 3

The process P 3 is defined by

P 3,i+1 = a i P 3i + b i , (18)
where a i is a uniformly distributed random variable that can take values larger than 1 -here, a i ∈ [0.48, 1.48] -and b i is uniformly distributed over [0, 1]. In order to ensure that P 3 does not grow to infinity but is "globally" stationary, we enforce the technical condition that the average growth rate ln a i be negative. For the parameters taken above, it follows that ln a i = -0.06747. The apparently innocuous self-affine map (18) leads to time series with surprisingly rich properties [START_REF] Sornette | Convergent multiplicative processes repelled from zero: power laws and truncated power laws[END_REF][START_REF] Sornette | Linear stochastic dynamics with nonlinear fractal properties[END_REF]. The probability density function of P has an asymptotic tail for large P of the form of a power law p(P ) ∼ P -(1+µ) , where the exponent µ is the smallest real solution of [a i ] µ = 1. Complex solutions µ of this probabilistic equation also exist that have equal or larger real parts; they correspond to higher-order corrections to the leading power-law behavior of p(P ) with logperiodic modulations. In practice, however, these modulations are not always observable due to their averaging by the large fluctuations in the a i and b i values [START_REF] Jögi | Fine structure and complex exponents in power law distributions from random maps[END_REF].

In the present numerical example, the smallest real exponent is µ ≈ 1.5. Since µ > 1, the mean P 3i is well-defined and finite. However, since µ < 2, the variance [P 3i -P 3i ]

2 is infinite. In practice, this means that, for a given realization and a fixed time length T , the estimation of the variance by integration over [t 0 , t 0 + T ] is very unstable as a function of the position t 0 of the time window and its length T . Furthermore, the variance diverges approximately as T 2-µ with the length T of the time series. Figure 16 : Wavelet analysis of P 3 with a Gaussian wavelet and its derivative. The scales vary from a = 2 to 32 in powers of 2; same layout as for Fig. 7.

Figure 17 : MS-SSA of P 3 with RC-1 and RC-2. The scales W vary from 4 to 64, in powers of 2; same layout as in Fig. 8.

Figure 1 :

 1 Figure 1: A 1000-point realization of the P1 process (top) and its blow-up (bottom).

Figure 2 :

 2 Figure 2: Log-periodic time series P 2 .

Figure 4 :

 4 Figure 4: SOI variations between 1933 and 1996.

Figure 5 :

 5 Figure 5: Monte-Carlo SSA of P1 .

Figure 6 :

 6 Figure 6: Global RCs of P1 .

Figure 7 :

 7 Figure 7: Wavelet transforms of P1 with scale increments of 3.

Figure 8 :

 8 Figure 8: RCs 1 and 2 of P1 for scale increments of 3.

Figure 9 :

 9 Figure 9: EOFs 1 -4 of P1 at b = 333, with a scaling of 3 n/2 .

Figure 10 :

 10 Figure 10: EOFs 1 -4 of P1 at b = 333, with a scaling of 2 n/2 .

Figure 11 :

 11 Figure 11: RCs 1 and 2 of P1 for scale increments of 3.

Figure 15 :

 15 Figure 15: Global RCs of P 3 (t).

Figure 17 :

 17 Figure 17: RCs 1 and 2 of P 3 with scale increments of 2.

Figure 19 :

 19 Figure 19: Instantaneous SOI frequencies, based on EOF-1 of MS-SSA at scale W = 128.

  It is a pleasure to thank A. Fournier, M. Mann, J.-F. Muzy and N. Saito for stimulating conversations and references. This work was supported by NSF grant ATM95-23787 (MG and PY), a NATO travel grant (PY) and research funds of UCLA's Institute of Geophysics and Planetary Physics (IGPP) during PY's leave from LSCE-CEA, and NSF grants EAR96-15357 and EAR97-06488 (DS). This is contribution No. XXXX of LSCE and contribution No. YYYY of IGPP at UCLA.

7

  Appendix: Construction of the synthetic time series 7.1 Approximate triadic Cantor set

Figure 1 :

 1 Figure 1 : Iterated funtion system (IFS) simulation of a classical triadic Cantor set. (a) 1000-point interpolation to generate the P1 time series; (b) blow-up of P1 between 1 and 333 in order to visualize the (incomplete) triadic structure of the time series.

Figure 2 :

 2 Figure 2 : Log-periodic time series P 2 , with critical point t c = 100, regularly sampled with 1000 points; see text and Eq. (17) for parameter values.

Figure 3 :

 3 Figure 3 : Multiplicative noise process P 3 : (a) individual realization over 512 equidistant points; (b) histogram of the variations of P 3 , with 100 bins -the ordinate represents the range of values (same as in panel a) and the abscissa is a frequency count.

Figure 4 :

 4 Figure 4 : Variations of the Southern Oscillation Index (SOI) between January 1933 and December 1996 (from the Climate Research Unit Data at the University of East Anglia, U.K.). Time on the abscissa in years and SOI on the ordinate centered to have mean zero and normalized by its standard deviation.

Figure 5 :

 5 Figure 5 : "Global" Monte-Carlo SSA of P1 , with window width (i.e. number of lags) M = 40.The vertical bars are 95% confidence intervals for 100 red-noise realizations with the same variance and exponential-decorrelation time as P1 (t). The diamonds represent the singular values of the signal; those that stand out above the error bars are statistically significant. The frequencies associated with the latter are obtained by least-square fitting the corresponding EOF to a sine function.

Figure 6 :

 6 Figure 6 : Reconstructed components (RCs) of P1 with global SSA. (a-c) Reconstructions using groups of eigenelements, as indicated in the legend of each panel; (d) the sum of RCs 1-6. The SSA window is M = 40, as in Fig. 5.

Figure 7 :

 7 Figure 7 : Wavelet transforms of P1 , with scale increments of 3. The left panels use the Gaussian analyzing function ψ(x) (shown as "Psi" in the legend) and the right ones use its derivative ψ ′ (x) (shown as "dPsi"); the scales a = 3, 9, 27 and 81 (for n = 1-4, see text) are also indicated from top to bottom on the left.

Figure 8 :

 8 Figure8: Multi-scale SSA (MS-SSA) of P1 , with scale increments of 3. The left panels show the reconstruction of P1 using the local RCs associated with EOF-1, the right panels that using the RCs obtained by projection onto the local EOF-2; the window widths W = 9, 27, 81 and 243 are shown from top to bottom on the left (like the scale a in Fig.7).

Figure 9 :

 9 Figure 9 : Normalized EOFs 1-4 (panels a-d) of P1 at b = 333. The EOFs at scale 3 n were multiplied by a factor of 3 (n-1)/2 ; the abscissa in each panel is linearly normalized to [0, 1]. The scales are W = 9, 27, 81 and 243, respectively in thick continuous, thin continuous, dotted and dashed lines (see also legend in panel a).

Figure 10 :

 10 Figure 10 : Normalized EOFs 1-4 of P1 at b = 333. The EOFs at scale 2 n are multiplied by a factor of 2 (n-1)/2 ; abscissa normalized as in Fig. 9. The scales are W = 16, 32, 64, 128 and 256, respectively in thick continuous, thin continuous, dotted, dashed and dot-dashed lines (see also legend in panel b).

Figure 11 :

 11 Figure11: MS-SSA of P1 , with scale increments of 3. The left panels contain the RCs associated with EOF-1, the right panels those associated with EOF-2; the same window widths W as in Fig.8(shown on the left).

Figure 12 :

 12 Figure12: Reconstruction of P 1 using formula[START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF] for scales that are powers of 3. The top panelshows the time series P 1 derived from the exact triadic Cantor set. The remaining panels show the reconstructions of P 1 using the local RCs based on EOF-2. The second panel uses RC-2 for the single scale N = 9; this corresponds to Eq. (11) with a single value of m = 2, a = 3 and thus W 2 = 9. The third panel shows the sum[START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF] with two terms, corresponding to W 2 = 9 and W 3 = 27; the fourth panel shows the sum[START_REF] Vautard | Singular spectrum analysis: A toolkit for short noisy chaotic signals[END_REF] with three terms, corresponding

Figure 13 :

 13 Figure 13 : Normalized EOFs 1-4 (panels a-d) of P 2 at t = t c -W/2, so that each window ends at t c . The EOFs at scale 2 n are multiplied by a factor of 2 (n-1)/2 , and the abscissa is normalized to [0, 1]. The scales are W = 16, 32, 64, 128 and 256, respectively in thick continuous, thin continuous, dotted, dashed and dash-dotted lines (see panel a).

Figure 14 :

 14 Figure 14 : MS-SSA of P 2 with RC-1 (left column) and RC-2 (right column). The scales W vary from 16 to 256, going from the top to the bottom (indicated along the left edge of the figure).

Figure 15 :

 15 Figure 15 : RCs of P 3 using global SSA. Panels (a-d) represent reconstructions that use groups of RCs as indicated to the left of each panel; panel e contains the sum of RCs 1-8.

Figure 18 :

 18 Figure 18 : Evolution in time of the local EOF-1's for P 3 , at W = 32 (M = W/3 = 10). The axes are the time and the lag.

Figure 19 :

 19 Figure 19 : Instantaneous frequencies of the local EOFs 1 (solid) and 2 (dotted) for the SOI time series in Fig. 4; the MS-SSA scale is W = 128 and the number of lags is M = W/3 = 42.