# Contents

## 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 . River discharge and rainfall similarly exhibit long periods of relative quiescence interspersed with large floods and bursts . 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 . 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 . The number of contagious-disease incidences is also intermittent, and has recently been documented to be distributed according to power laws . The evolutionary path of species on earth with its “punctuated dynamics” provides a paradigm of intermittent behavior . Economic and financial time series have long been recognized to have nonstationary variances, possibly described by Lévy processes .

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 . 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)  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 . 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  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  and simple discontinuities in the local value of an oscillation’s frequency can also be detected by SSA .

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 . 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 decomposition of a signal, it may be desirable to (A) automatically adjust the shape of the analyzing 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 $`𝐂`$ – provides data-adaptive wavelet transforms with analyzing functions given by the leading EOFs of the matrix $`𝐂`$; 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 . A summary of our results and comparison with complementary approaches to address the requirements (A) and (B) above appear in Sec. 5.

## 2 Methodology

### 2.1 Wavelet transform

A wavelet transform requires the choice of an analyzing function $`\psi `$, with general admissibility properties , and with the more specific property of time and frequency localization, *i.e.*, $`\psi `$ and its Fourier transform $`\widehat{\psi }`$ must decay rapidly outside a given interval. Functions based on a Gaussian, $`f(x)=\mathrm{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) .

In the sequel, we shall follow  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  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 $`\psi `$-wavelet transform $`W_\psi `$ in continuous time and frequency is simply a projection of a signal $`X(t)`$, $`\mathrm{}<t<\mathrm{}`$, onto $`b`$-translated and $`a`$-dilated versions of $`\psi `$:

$`W_\psi (a,b)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{a}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}X(t)\psi \left({\displaystyle \frac{tb}{a}}\right)𝑑t`$ (1)   
$`=`$ $`{\displaystyle \frac{1}{\sqrt{a}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}X(t+b)\psi \left({\displaystyle \frac{t}{a}}\right)𝑑t.`$ (2)   
If most of $`\psi `$ 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 $`\psi ^{(n)}`$ of a given wavelet $`\psi `$ in Eq. (1) is equivalent to a $`\psi `$-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 $`\psi `$, that the number of oscillations in $`\psi ^{(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_\psi _{aI}_{b=\mathrm{}}^{\mathrm{}}W_\psi (a,b)\psi \left(\frac{tb}{a}\right)\frac{dadb}{a^2};$$
(3)   
$`A_\psi `$ is a normalization factor which only depends on the mother wavelet $`\psi `$. This formulation is essentially a bandpass filter of $`X`$ through $`I`$. If $`I`$ is the positive real line, $`I=𝐑^+`$, 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 $`\psi (x)=\mathrm{exp}(x^2/2)`$ itself cannot be used in this reconstruction formula, because it does not satisfy $`\psi (x)𝑑x=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\mathrm{\Delta }t`$, $`0iN`$, there are numerous ways of estimating Eq. (1) numerically that avoid the problems of high-frequency noise and aliasing, while minimizing computational cost . Computationally efficient algorithms based on the fast Fourier transform can be used , but we will only need here the simplest numerical quadrature formulae to process Eq. (1).

### 2.2 Singular spectrum analysis (SSA)

SSA has been described theoretically, for time-continuous signals $`X(t)`$, by Vautard and Ghil . 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\mathrm{\Delta }t),0iN\}`$, which is generally used in practice, as it will be in the present study. Methodological and review papers include . SSA proceeds by finding the eigenvectors, called EOFs, of the lag-covariance matrix $`𝐂`$ of the time series $`\{X_i\}`$ as above. The true lag-covariance matrix $`\stackrel{~}{𝐂}`$ of the process that generated $`\{X_i:1iN\}`$ is estimated by a lag-covariance matrix $`𝐂_M`$; assuming stationarity of $`X`$ yields a Toeplitz structure

$$𝐂_{M,i}=\frac{1}{Ni}\underset{j=1}{\overset{Ni}{}}X_jX_{j+i},$$
(4)   
where $`i`$ ($`0iM1`$) 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 $`\rho _k`$ of $`𝐂_M`$ are sorted by decreasing order of the associated eigenvalue $`\lambda _k`$, $`0\lambda _M\mathrm{}\lambda _2\lambda _1`$. The principal component (PC) $`𝐚_k=\{a_{ki}:0iNM\}`$ is computed by projecting the time series onto $`\rho _k=\{\rho _{kj}:0jM\}`$:

$$a_{ki}=\underset{j=1}{\overset{M}{}}X_{i+j1}\rho _{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 $`\rho _k`$ have an oscillation property which resembles that of the derivatives $`\psi ^{(n)}`$: the higher the order $`k`$, the more sign changes $`\rho _k`$ has . This oscillation property is due to the formal analogy between the Toeplitz structure of the covariance matrix $`𝐂_M`$ and that of finite-difference discretizations of Sturm-Liouville problems on an interval $`W=M\mathrm{\Delta }t`$ . In our case, $`\mathrm{\Delta }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 $`𝒦`$, the corresponding reconstructed component (RC) is  :

$$r_{𝒦,i}=\frac{1}{M_i}\underset{k𝒦}{}\underset{j=1}{\overset{M}{}}a_{k,ij}\rho _{kj},$$
(6)   
where $`M_i`$ is a scaling factor equal to $`M`$ when $`i`$ is sufficiently far from the end points, and contains corrections when $`i`$ approaches the interval’s endpoints to within a distance $`M\mathrm{\Delta }t`$.

### 2.3 Multi-scale SSA (MS-SSA)

Systematic comparisons between SSA, the wavelet transform and other spectral-analysis methods have been carried out in Ref.  (see, in particular, Table 1 in ). Further analogies between certain mathematical features of SSA and wavelet analysis were mentioned in . 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 $`\rho _k`$ contains information from the whole time series, as in the Fourier transform.

In order to define a local SSA, we propose to extend the SSA methodology by using a time-frequency 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 $`WN`$, centered on times $`b=W/2,\mathrm{},NW/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\alpha >1`$, and $`\alpha `$ is large enough . 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=\alpha `$ 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” $`\rho _k`$ of each local SSA analysis; $`\alpha =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 $`\rho _{k+1}`$ separate those of $`\rho _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  or of Legendre polynomials . 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 $`\rho _k`$, it is possible to obtain local PCs $`a_k`$ and RCs $`r_k`$. The $`k`$th PC at time $`b`$ is

$$a_{ki}^b=\underset{j=1}{\overset{M}{}}X_{i+j1}\rho _{kj}^b,$$
(7)   
and the corresponding RC is

$$r_{ki}^b=\frac{1}{M_i}\underset{j=1}{\overset{M}{}}a_{k,ij}^b\rho _{kj}^b,$$
(8)   
with $`bW/2ib+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 $`NW/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, $`\{\rho _k:1kM\}`$, 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=\alpha 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 .

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=\beta \mathrm{\Delta }t`$. This sampling interval for the analysis must be substantially smaller than the scale parameter $`W`$, i.e. $`\beta W`$. Then we use Eq. (8) for $`b\beta \mathrm{\Delta }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_{kb}^b`$ 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 .

## 3 Data sets

### 3.1 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 .

The first time series $`\stackrel{~}{P}_1`$ that we shall investigate was obtained by using an Iterated Function System (IFS)  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 $`\stackrel{~}{P}_1`$. Finally, a multiple-rule Cantor set is generated in Appendix 7.2 to yield the time-series $`\widehat{P}_1`$.

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 , 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 . In the eighties, they have been shown to emerge in various physical problems that arise in discrete hierarchical systems .

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  (see  for a review, and references therein). Another source of time series with log-periodic behavior is provided by Boolean delay equations, whose solutions are Boolean-valued processes in continuous time . These equations have been used to model paleoclimatic variability , interdecadal  and most recently, seasonal-to-interannual  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 $`tt_c=+\mathrm{}`$  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,\stackrel{~}{P}_1)`$ 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 . A 512-point realization of $`P_{3i}`$, for $`1i512`$, is shown in Fig. 3 and a histogram is provided to illustrate the power-law behavior and its long tail.

### 3.2 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 quasi-periodic and aperiodic behavior . 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  and improving its prediction .

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 ocean-atmosphere system’s self-sustained oscillation to the annual cycle leads to a Devil’s staircase . 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 , which spans a longer time interval, using wavelet and waveform analysis, and Tropical Pacific sea-surface temperatures with global SSA . Two of these analyses have shown a fairly sudden shift in the characteristic period of ENSO’s low-frequency component that occurred in the 1960s, with a longer period before and a shorter one prevailing since .

## 4 Numerical results

### 4.1 Cantor-set experiments

#### 4.1.1 SSA analysis

Global methods of time-series analysis will clearly miss the multi-scale feature of the IFS construct $`\stackrel{~}{P}_1(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  of $`\stackrel{~}{P}_1(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 $`\stackrel{~}{P}_1`$; 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 amplitude of small-scale oscillations is exactly the same as the one of large-scale oscillations, and 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: ; Multi-Taper: ) do not provide more relevant information on $`\stackrel{~}{P}_1(t)`$ either (not shown).

#### 4.1.2 Wavelet transforms

We next performed wavelet transforms of $`\stackrel{~}{P}_1`$ with the Gaussian analyzing function $`\psi (x)=\mathrm{exp}(x^2/2)`$ and its derivative $`\psi ^{}`$. Geometric scale increments $`a=3^n`$, $`n=0,\mathrm{},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 . The wavelet analyses of $`\stackrel{~}{P}_1`$ 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 $`\psi ^{}`$ as the mother wavelet (right column of Fig. 7), rather than $`\psi `$ (left column), reveals further aspects of the process’s singularities. The triadic construction of $`\stackrel{~}{P}_1`$ is also apparent in this decomposition, albeit not as clearly as with $`\psi `$. The wavelet transforms of $`\stackrel{~}{P}_1`$ in Fig. 7 do not differ significantly from those of its “exact” counterpart $`P_1`$ (not shown).

#### 4.1.3 MS-SSA analysis

$``$ RC analysis. We performed MS-SSA analyses of $`\stackrel{~}{P}_1`$ with scales $`W`$ varying with geometric increments, $`W=23^n`$, and a ratio $`W/M=3`$. We kept the first two EOFs for each scale. The first RCs of $`\stackrel{~}{P}_1`$ 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 $`\stackrel{~}{P}_1`$, i.e., the divisions do not overlap, even though parts of the reconstruction are lost due to side effects. Thus, in the case of $`\stackrel{~}{P}_1`$, 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 $`\psi ^{}`$-wavelet transform. Hence, EOF-2 is analogous to $`\psi ^{}`$ (see Table 1 and ref.), and the resulting RCs yield an analysis of the derivative of $`\stackrel{~}{P}_1`$ (not shown).

$``$ EOF analysis. As anticipated in Section 2, we expect the EOFs $`\rho _k^{(W)}`$ 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 $`\stackrel{~}{P}_1`$, 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 hand, needs larger scales (and more data) in order to converge.

In contrast, we show in Fig. 10 the same EOFs 1-4 of $`\stackrel{~}{P}_1`$, 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 $`\lambda `$ in a time series by optimizing the scale ratio $`l`$ in MS-SSA: only when $`l\lambda ^m`$, where $`m`$ is an integer, will the convergence of the EOFs $`\rho _k^{(W)}\rho _k^{(\mathrm{})}`$ be good as $`W\mathrm{}`$.

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 $`\psi `$. In MS-SSA, the functional-shape constraint is relaxed and new information on the signal can thus be obtained.

$``$ 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 $`\stackrel{~}{P}_1(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. Well-known instances are the randomization method in general probability theory  and the addition of process or observation noise in sequential estimation theory .

$``$ Change of geometrical structure along the scale axis. MS-SSA analysis of the multiple-rule Cantor set $`\widehat{P}_1`$ (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 $`\stackrel{~}{P}_1`$ 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 (14) for the larger scales and rule (16) 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 $`\psi ^{(n)},0nM1`$, of the same mother wavelet $`\psi =\psi ^{(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 . For discrete wavelet transforms, there does not exist, in general, a “resolution of identity” formula analogous to (3): the expansion of a function over a discrete wavelet basis is not orthogonal in general and special iterative methods must be developed to derive an inversion formula (see Chapter 3 in ).

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” . It is easy to show from the definition of a frame that the family of EOFs $`\rho _k^W`$ 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`$  :

$$X_i=\underset{m}{}\underset{j=1}{\overset{M}{}}\widehat{\rho }_{k,j}^{W_m}a_{k,ij}^{W_m};$$
(9)   
here the sum over $`m`$ corresponds to the different scales, and $`a_{k,ij}^{W_m}`$ is given by

$$a_{k,ij}^{W_m}=\underset{l=1}{\overset{M}{}}\rho _{k,l}^{W_m}X_{ij+l1},$$
(10)   
while $`\widehat{\rho }_{k,j}^{W_m}`$ is a vector related to the EOF $`\rho _{k,j}^{W_m}`$ through an iterative algorithm .

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 $`\widehat{\rho }_{k,j}^{W_m}`$ replaced by $`\rho _{k,j}^{W_m}`$ :

$$X_i=\underset{m}{}\underset{j=1}{\overset{M}{}}\frac{1}{M}\rho _{k,j}^{W_m}a_{k,ij}^{W_m}.$$
(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 (11) with a single value of $`k=2`$, while $`a=3`$ and thus $`W_2=9`$. The third, fourth and fifth panels show the sum (11) 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 $`\rho _k^{W_m}`$ is always $`k=2`$.

### 4.2 A log-periodic process

In this subsection, we focus on the detection and quantification of self-similar structures in a log-periodic process. In the case of Eq. (17), we expect discrete self-similarity near $`t=t_c`$, with a preferred scaling of $`\lambda =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 $`\psi `$ and $`\psi ^{}`$) 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 window — and 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.

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).

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  and financial  time series will be investigated separately.

### 4.3 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`$.

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 $`\psi ^{}`$ with respect to $`\psi `$ (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 $`\stackrel{~}{P}_1`$, 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 multifractal . The shape of the EOFs (not shown) is in fact quite similar to that of the Gaussian 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`$.

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.

### 4.4 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 , 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 ? 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.  in Tropical Pacific sea-surface temperatures, by using multi-channel global SSA, and by Wang and Wang  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  and the wavelet and waveform (matching pursuit) approach , 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 .

## 5 Concluding remarks

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 examples – between 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.  describe a hierarchical binary tree of basis functions for locally stationary time series. Lilly and Park  select their Slepian wavelets  at a given scale by optimizing the spectral energy in a given frequency band  and vary the scale. Mann and Park  have applied this approach to a multivariate analysis of interannual temperature variations. Wavelet packets  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  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” .

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 $`\stackrel{~}{P}_1`$, 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  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  and standard wavelet  analysis. It supports the Devil’s staircase theory of ENSO  in the following sense : a sudden shift in frequency can be reconciled most easily with the idea that an interdecadal change in mean thermocline depth in the Tropical Pacific  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.

## 6 Acknowledgments

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

The time series $`\stackrel{~}{P}_1`$ was obtained by using an Iterated Function System (IFS)  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:

$$x_{n+1}=\{\begin{array}{cc}\frac{1}{3}x_n\hfill & \text{with probability }\frac{1}{2},\hfill \\ \frac{1}{3}x_n+\frac{2}{3}\hfill & \text{with probability }\frac{1}{2}.\hfill \end{array}$$
(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 $`\mathrm{log}2/\mathrm{log}3`$ . 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, $`0x_0\mathrm{}x_ix_{i+1}\mathrm{}x_N1`$. A time series $`\stackrel{~}{P}_1`$ with, say $`N=1000`$ equidistant points, is generated by interpolating the set $`\{x_i\}`$ every $`\mathrm{\Delta }t=1/N`$ and scaling it to $`N`$:

$$\stackrel{~}{P}_{1j}=\{\begin{array}{cc}1\hfill & \text{if }x_i[j/N,(j+1)/N],\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}$$
(13)   
with $`0jN`$. 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.

### 7.2 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)`$, $`1j3^L`$. A “noisy” Cantor time series could be obtained by assigning a probability of occurence between rule (14) and $`1111`$. This procedure emulates the incomplete division scheme towards smaller scales that we observe in $`\stackrel{~}{P}_1`$.

A large variety of Cantor sets can be obtained by modifying the rules (14, 15) 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,`$ (16)   
$`0`$ $``$ $`000.`$   
three more times, in order to obtain a set $`\widehat{P}_1`$ 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 $`\mathrm{log}2/\mathrm{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 . 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.

### 7.3 Log-periodic time series

The functional form of the log-periodic process $`P_2`$ is simply:

$$P_2(t)=A+B(tt_c)^\alpha \left\{1+C\mathrm{cos}\left[2\pi \mathrm{log}(tt_c)/\mathrm{log}(\lambda )\right]\right\}.$$
(17)   
Here, we take $`A=2`$, $`B=1`$, $`C=0.2`$, $`\alpha =1/2`$ and $`\lambda =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 $`\mathrm{\Delta }t=0.1`$, and define the $`P_2`$ time series by $`P_{2i}=P_2(i\mathrm{\Delta }t)`$, $`1i1000`$. Note that in our case $`\alpha =1/2`$ and thus $`lim_{tt_c}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_iP_{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 $`\mathrm{ln}a_i`$ be negative. For the parameters taken above, it follows that $`\mathrm{ln}a_i=0.06747`$.

The apparently innocuous self-affine map (18) leads to time series with surprisingly rich properties . The probability density function of $`P`$ has an asymptotic tail for large $`P`$ of the form of a power law $`p(P)P^{(1+\mu )}`$, where the exponent $`\mu `$ is the smallest real solution of $`[a_i]^\mu =1`$. Complex solutions $`\mu `$ 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 log-periodic 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 .

In the present numerical example, the smallest real exponent is $`\mu 1.5`$. Since $`\mu >1`$, the mean $`P_{3i}`$ is well-defined and finite. However, since $`\mu <2`$, the variance $`\left[P_{3i}P_{3i}\right]^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\mu }`$ with the length $`T`$ of the time series.

## 8 Figure captions

: Iterated funtion system (IFS) simulation of a classical triadic Cantor set. (a) 1000-point interpolation to generate the $`\stackrel{~}{P}_1`$ time series; (b) blow-up of $`\stackrel{~}{P}_1`$ between 1 and 333 in order to visualize the (incomplete) triadic structure of the time series.

: 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.

: 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.

: 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.

: “Global” Monte-Carlo SSA of $`\stackrel{~}{P}_1`$, 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 $`\stackrel{~}{P}_1(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.

: Reconstructed components (RCs) of $`\stackrel{~}{P}_1`$ 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.

: Wavelet transforms of $`\stackrel{~}{P}_1`$, with scale increments of 3. The left panels use the Gaussian analyzing function $`\psi (x)`$ (shown as “Psi” in the legend) and the right ones use its derivative $`\psi ^{}(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.

: Multi-scale SSA (MS-SSA) of $`\stackrel{~}{P}_1`$, with scale increments of 3. The left panels show the reconstruction of $`\stackrel{~}{P}_1`$ 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).

: Normalized EOFs 1–4 (panels a–d) of $`\stackrel{~}{P}_1`$ at $`b=333`$. The EOFs at scale $`3^n`$ were multiplied by a factor of $`3^{(n1)/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).

: Normalized EOFs 1–4 of $`\stackrel{~}{P}_1`$ at $`b=333`$. The EOFs at scale $`2^n`$ are multiplied by a factor of $`2^{(n1)/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).

: MS-SSA of $`\widehat{P}_1`$, 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).

: Reconstruction of $`P_1`$ using formula (11) for scales that are powers of 3. The top panel shows 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 (11) with two terms, corresponding to $`W_2=9`$ and $`W_3=27`$; the fourth panel shows the sum (11) with three terms, corresponding to $`W_2=9`$, $`W_3=27`$ and $`W_4=81`$; the bottom panel shows the sum (11) with four terms, corresponding to $`W_2=9`$, $`W_3=27`$, $`W_4=81`$ and $`W_5=243`$ and still the same order $`k=2`$ of the EOF.

: Normalized EOFs 1–4 (panels a-d) of $`P_2`$ at $`t=t_cW/2`$, so that each window ends at $`t_c`$. The EOFs at scale $`2^n`$ are multiplied by a factor of $`2^{(n1)/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).

: 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).

: 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.

: 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.

: 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.

: 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.

: 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`$.