Cross wavelet analysis: significance testing and pitfalls

<jats:p>Abstract. In this paper, we present a detailed evaluation of cross wavelet analysis of bivariate time series. We develop a statistical test for zero wavelet coherency based on Monte Carlo simulations. If at least one of the two processes considered is Gaussian white noise, an approximative formula for the critical value can be utilized. In a second part, typical pitfalls of wavelet cross spectra and wavelet coherency are discussed. The wavelet cross spectrum appears to be not suitable for significance testing the interrelation between two processes. Instead, one should rather apply wavelet coherency. Furthermore we investigate problems due to multiple testing. Based on these results, we show that coherency between ENSO and NAO is an artefact for most of the time from 1900 to 1995. However, during a distinct period from around 1920 to 1940, significant coherency between the two phenomena occurs.                     </jats:p>

section_content: Time series analysis is a fundamental issue in climatology as in many other fields of empirical research (von Storch and Zwiers, 1999; Ghil et al., 2002). Considering climate records one almost always faces a composition of numerous scales ranging from days to decades or even longer periods. The classical method to investigate such data by frequency decomposition is Fourier analysis (Priestley, 1992; Brockwell and Davis, 1987). On the considered time scales climate processes are often non-stationary and time resolved methods become necessary. The straightforward time resolved extension of Fourier analysis is a windowed or sliding Fourier Transformation (Kaiser, 1994). A disadvantage of this method is, that the window width and thus the time resolution is constant for all investigated frequencies. Here, continuous wavelet analysis (Daubechies, 1992 ) has proven to be superior: The time resolution is intrinsically adjusted Correspondence to: D. Maraun (maraun@agnld. uni-potsdam. de) to the scales. Thus optimal time resolution for every scale is given. Another advantage of wavelet analysis is the flexible choice of the mother wavelet according to the characteristics of the investigated time series (in the following, we always omit the term "continuous"). For an introduction to discrete wavelet transformation and its applications see Kaiser (1994) and Percival and Walden (2000). 

During the last decade, wavelet analysis has been successfully applied in climate science, e. g. to El Niño/Southern Oscillation (ENSO) data (Gu and Philander, 1995; Kestin et al., 1998; Torrence and Webster, 1998), NAO proxy data (Appenzeller et al., 1998) or output of globally coupled ocean atmosphere models (Timmermann et al., 1999). A good introduction to the basic theory and application of wavelet analysis is given by Torrence and Compo (1998). The most important contribution of the latter paper is developing detailed significance tests for the wavelet power spectrum to address the critics, that wavelet analysis is rather producing colorful images than reliable results. 

An important field in time series analysis, especially in climate sciences, is multivariate analysis. When comparing two different variables like temperature or pressure, or when analyzing teleconnections, one needs the bivariate extension of wavelet analysis. This cross wavelet analysis, which was introduced by Hudgins et al. (1993) and others, is briefly discussed by Torrence and Compo (1998). It has already been applied to various problems in climatology like the ENSO-Monsoon system (Torrence and Webster, 1999), ENSO-North Atlantic Oscillation (NAO) teleconnections (Huang et al., 1998), a rainfall-runoff cross analysis (Labat et al., 2000) or the influence of NAO on European surface temperatures (Pozo-Vazquez et al., 2001). 

However, important questions concerning this technique remain open. We try to contribute to this discussion by investigating the wavelet cross spectrum (WCS) and the normalized wavelet coherency (WCO). For WCO, we develop a statistical test for zero coherency. Since test statistics cannot easily be calculated analytically, we apply Monte Carlo simulations. For the case that at least one of the processes can be modeled by Gaussian white noise, we present an approximative formula to easily calculate the critical value for significance on the 95% level. A second focus of the paper are typical pitfalls of cross wavelet analysis. We show, that not every structure in a wavelet plot has got a physical meaning but might be an artefact of WCS or multiple testing. To illustrate our discussion, we show that the coherency between ENSO and NAO for most of the moderate and strong El Niños between 1900 and 1995 as suggested by Huang et al. (1998) is an artefact of WCS. However, during a distinct period from around 1920 to 1940 significant coherency between the two phenomena occurs and invites to further investigation. The paper is organized as follows: In Sect. 2 we introduce the most important measures for wavelet analysis. The problem of normalizing wavelet spectra is briefly explained in Sect. 3. In Sect. 4, we interpret wavelet analysis as an inverse problem of parameter estimation and discuss the consequences arising. Sect. 5 provides a numerical significance test for time-and scale-dependent coherency. Important pitfalls related to cross wavelet analysis are studied in Sect. 6. We review the cross wavelet analysis of ENSO and NAO in Sect. 7.
section_name: Introduction

section_name: Definitions

section_content: The wavelet transformation W i (s) at time t i =i t on a scale s of a discrete time series x j =x(t j ) of length N with a sampling interval t can be interpreted as an extension of the discrete Fourier transformation F (ω)= j x j exp(iωt j ) (Kaiser, 1994). The latter one extracts contributions on frequency ω. Wavelet transformation replaces the periodic exponential exp(iωt j ) with a localized wavelet (t j -t i, s), which is located around the time t i and stretched according to the investigated scale s: Thus the time series can be decomposed scale-and time-dependent:

x j ((ji) t, s). 

(1)

If one considers arbitrary scales between the sampling interval and the length of the time series, one speaks of continuous wavelet transformation (CWT). In this case, the wavelets do not span an orthogonal set of basis functions. Hence neighboring scales and times contain redundant information and are correlated. 

The wavelet (t j -t i, s) is a stretched and translated version of a chosen mother wavelet, normalized with a factor c(s) (see Sect. 3):

(2)
section_name: Continuous wavelet transformation

section_content: In the following, we always consider the Morlet mother wavelet (see Fig. 1 )

where θ and ω 0 are unit-less. The Gaussian envelope exp(-θ 2 /2) localizes the wavelet in time. The time/scale resolution is adjusted by ω 0. For high values of ω 0, the scale resolution increases, whereas time resolution decreases and vice versa. Fourier frequency f and wavelet scale s are not direct reciprocals of each other. Instead, one has to rescale the result of wavelet analysis with a factor depending on the mother wavelet (Meyers et al., 1993; Torrence and Compo, 1998). For the Morlet wavelet, the conversion reads

For ω 0 =6, s•f is approximately one. For other mother wavelets refer to Kaiser (1994) and Torrence and Compo (1998).
section_name: Morlet wavelet

section_content: Transferring a concept of Fourier analysis, a wavelet power spectrum (WPS) can be defined as the wavelet transformation of the autocorrelation function. Following the Wiener-Khinchin theorem, this can be implemented as follows:

The brackets denote expectation values. WPS i (s) describes the power of the signal x(t) at a certain time t i on a scale s. 

The relevance of considering expectation values and how to estimate them is discussed in Sect. 4.
section_name: Wavelet power spectrum

section_content: As in Fourier analysis, the univariate WPS can be extended to compare two time series x(t i ) and y(t i ): One can define the wavelet cross spectrum WCS i (s) as the expectation value of the product of the corresponding W x i (s) and W y i (s):

In Sect. 6 we show, that WCS-peaks can even appear when x(t i ) and y(t i ) are independent, because WCS is not normalized. In contrast to the WPS, the WCS is complex valued (analogous to Fourier cross spectra) and can be decomposed into amplitude |WCS i (s)| and phase i (s):

The phase i (s) describes the delay between the two signals at time t i on a scale s. 

A normalized time and scale resolved measure for the relationship between two time series x(t i ) and y(t i ) is the wavelet coherency (WCO), which is defined as the amplitude of the WCS normalized to the two single WPS:

A value of 1 means a linear relationship between x(t i ) and y(t i ) around time t i on a scale s. A value of zero is obtained for vanishing correlation. Because of the normalization, spurious peaks can occur for areas of low wavelet power. Note the importance of estimating the expectation values of the single terms in Eq. ( 8 ): If they are not considered separately, nominator and denominator equal each other and the result will be identical to 1 no matter what the real WCO is. Thus computing WCO from observational data, one has to smooth nominator and denominator separately (see also Sect. 4. 2).
section_name: Cross wavelet analysis

section_content: In Fourier analysis, any normalization automatically provides the following two features:

-The Gaussian white noise spectrum is (by definition) flat. 

-Sines of the same amplitude have the same integrated power in the frequency domain. 

In wavelet analysis, we meet difficulties to obtain both. For the factor c(s) in Eq. ( 1 ), Torrence and Compo (1998) and Kaiser (1994) suggest different normalizations, which only preserve one of the mentioned features. Torrence suggests

which preserves a flat white noise spectrum, but sines of equal amplitude exhibit different integrated power proportional to their oscillation scale. This choice equals Kaiser's normalization with p = 1/2 (Eq. 3. 5 in Kaiser, 1994, p. 62). Wavelet power spectrum of Gaussian white noise and two sines of frequency 1/2 and 1/4 and equal amplitude. Solid line: Normalization according to Torrence and Compo (1998) does not preserve the integrated power, but a flat white noise spectrum. Dashed line: Normalization according to Kaiser (1994) preserves the integrated power, but the white noise spectrum is ∼1/scale. 

Using the normalization of Kaiser with p = 0,

sines of equal amplitude also exhibit equal integrated power (when scale is plotted logarithmically). On the other hand, the white noise spectrum is no longer flat but decreases proportional to 1/scale. Figure 2 illustrates this dilemma of different normalizations for time series of Gaussian white noise and a sine function. 

Fortunately, the normalization is only relevant for a first inspection by eye. A normalization following Eq. ( 9 ) emphasizes power on high scales and could lead to misinterpretations. However, when performing a significance test, the significance level already includes the chosen normalization.
section_name: Normalizations

section_name: Wavelet analysis of observational records

section_content: In time series analysis, one aims to reveal properties of an underlying process from a given observation. This is referred to as an inverse problem (Honerkamp, 1998), having two consequences: First, since most processes in geosciences are mixing (i. e. loosely spoken processes that lose the knowledge about their actual state for t→∞, see e. g. Karlin and Tylor, 1975) and additionally all measured time series are of finite size, one will never be able to obtain the exact values of the "true" properties. Instead one has to define a meaningful estimator, which is a function of the observation and (if possible) would yield the exact desired property if the time series was of infinite lengths. Secondly, in general it is not a priori clear if the inference of a process property from observational data is unambiguous (Moritz, 1993; Anger, 1993).
section_name: The inverse problem

section_content: The observed WPS of a measured time series is an estimator for the "true" WPS of the underlying process. A comparison with Fourier analysis illustrates the involved consequences:

For the Fourier spectrum, the first guess estimator is the periodogram, which is the squared absolute value of the Fourier transformation of a given time series (Priestley, 1992; Brockwell and Davis, 1987). If this time series is a realization of a mixing process, every value of the periodogram is χ 2 distributed with two degrees of freedom. The variance of this estimator will not converge to zero for an infinite time series, i. e. the periodogram is not a consistent estimator. To construct a consistent estimator for the true spectrum, one has to estimate the expectation value of the periodogram. In many practical cases, when only one realization is given, this has to be done by smoothing neighboring frequencies. 

This line of argumentation can be transferred to wavelet analysis: Commonly, as an estimator for the true WPS, the squared absolute value of the CWT is calculated (i. e. Eq. 5 without taking the expectation value). This measure is called wavelet periodogram (Nason et al., 2000) to emphasize that it is like the Fourier periodogram χ 2 -distributed with 2 degrees of freedom and thus not a consistent estimator. This does not contradict the finding of Percival and Walden (2000), who state that the global wavelet spectrum is a consistent estimator of the time independent spectrum. There the expectation value is estimated by averaging over all times. 

In wavelet analysis, the situation is quite subtle: Since neighboring points in time and scale are correlated (see Appendix A), the wavelet periodogram looks smooth. However, fluctuations around the true WPS are not smaller as in Fourier analysis, but vary jointly. Hence, the appearance of smooth patterns can suggest structures which are just coincidence. This has to be kept in mind, when interpreting wavelet plots. Ideally, the expectation value should be estimated. When, as usually in climatology, only one realization is measured, one either has to assume stationarity over a certain time interval and smooth in time direction, or one has to assume that neighboring scales exhibit similar power and thus smooth in scale direction (Sect. 5). For WCO the smoothing is essential, as shown in Sect. 2. 4.
section_name: Estimator for wavelet power

section_content: The wavelet transformation at a point in time t 0 always contains information of neighboring data points. The number of these points depends on the chosen wavelet and the scale considered. Thus, if the wavelet is centered close to the beginning or the end of the time series, edge effects occur. The area, where such effects are relevant, is called the cone of influence. Here, the results should be interpreted carefully. In the figures of this paper, the cone of influence is chosen as the e-folding time of the Morlet wavelet
section_name: Edge effects

section_content: 2s and marked as a shadow in the wavelet plot. For a discussion, refer to Torrence and Compo (1998).
section_name: √

section_content: The wavelet transform of the realization of a mixing process is a random number. Thus, even for independent processes, nonzero values of the coherency will result and the question arises, when a value has to be considered as being different from zero. The idea is to formulate a null hypothesis H 0 "the processes are not coherent" and to derive a significance test for H 0. Therefore, one has to estimate the probability distribution of the coherency under H 0. If the coherency lays outside a chosen α-quantile, the so-called critical value of this distribution, one rejects H 0 and the coherency is considered to be significantly different from zero on the (1 -α)level. For an introduction to significance testing refer to Honerkamp (1994), Lehmann (1986) or Cox and Hinkley (1994). 

In the case of Fourier analysis, a test for significant coherency can be calculated straightforwardly (Brockwell and Davis, 1987) : Irrespective of the two processes to be analyzed, the Fourier transform is asymptotically Gaussian distributed and thus the cross periodogram is χ 2 distributed. Since vicinal frequencies are asymptotically uncorrelated, the coherency thus follows a F distribution. 

However, this is not the case for different wavelet scales and times. Firstly, for small scales the asymptotic -process independent -distribution is not reached. Thus, considering small scales, one has to find suitable models for the two processes to be analyzed. Secondly, vicinal wavelet times and scales are not uncorrelated (Appendix A) and we show in Appendix B, that the WCS is not χ 2 distributed. Thus an analytical test statistic is highly non trivial if not impossible. 

To overcome the second problem, we performed Monte Carlo simulations to estimate the distribution under H 0 numerically. We first considered the most basic case: We simulated 10 000 realizations of two independent Gaussian white noise processes and estimated the 95% critical values for different smoothing windows in time and scale and for various time/scale resolutions defined by ω 0. Since WCO is a local and normalized measure, instationarity of the investigated process does not require instationary models to derive the test statistics. 

Given a choice of the smoothing window in time and scale, a scale independent critical value is desirable to ease the analysis (since white noise is stationary in time, the critical value will automatically be time independent). We found that this is obtained for the following smoothing procedures:

-When smoothing in time direction, one has to take into account the wavelet length of essential support, which is proportional to the scale s considered (see Eq. 3). Thus, a scale independent critical value is obtained, when the smoothing length is related to the actual scale s by a constant factor m. 

-Given a wavelet transform calculated for scales which are equidistant in logarithmic scale (i. e. more values for smaller scales), one has to smooth a constant number w of neighboring scales to obtain a scale independent critical value. 

(a) Figure 3 presents the 95% critical values WCO 95 obtained from the Monte Carlo simulations of two Gaussian white noise models. We use the Morlet wavelet and a fixed ω 0 =6 for different smoothing windows in (a) time direction and (b) scale direction respectively. The efficiency of smoothing is rather weak due to the strong correlation between neighboring times and scales (see Appendix A). Furthermore, smoothing in time direction appears to be much less effective than smoothing in scale direction (compare the values at the ordinates). This finding becomes intuitively clear, when we recall that wavelets already intrinsically smooth in time direction. As expected, we found the normalization discussed in Sect. 3 having no influence on the significance test, since different normalizations cancel down in WCO. The dependency of WCO 95 for the two Gaussian white noise models on the relevant parameters ω 0 and the smoothing length in scale direction w is shown in Fig. 4. The plotted values for WCO 95 are the mean values over all scales of the Empirical 95% critical values for 10 000 realizations of wavelet coherency of two Gaussian white noise time series for smoothing in scale direction with different smoothing window lengths w and different ω 0. The plotted quantiles WCO 95 are the mean values over all scales plotted in Fig. 3b, but for various ω 0. 

Monte Carlo simulations for a fixed pair of parameters ω 0 and w as depicted for ω 0 =6 in Fig. 3b. For practical purposes, we fitted a polynomial to these results and obtained the following approximation for significance testing WCO using the Morlet wavelet and smoothing in scale direction: WCO 95 = 8. 23 ×10 -5 ω 3 0 + 4. 24 ×10 -5 ω 2 0 w + 1. 13 ×10 -5 ω 0 w 2 + 1. 54 ×10 -5 w 3 -2. 30 ×10 -3 ω 2 0 -2. 19 ×10 -3 ω 0 w -7. 51 ×10 -4 w 2 + 2. 05 ×10 -2 ω 0 + 1. 27 ×10 -2 w + 0. 95, (11) where 2w+1 denotes the width of the smoothing window and ω 0 defines the time/scale resolution. This approximation is only valid in the interval shown in Fig. 4. Since smoothing in time direction appears to be too ineffective, we did not consider it in detail. 

The Gaussian white noise model is a simple and often unrealistic assumption. As mentioned above, for Fourier analysis critical coherency values are asymptotically independent of the model. This holds also for wavelet analysis in the limit of large scales s. However, since wavelet analysis is localized in time, for small scales s the asymptotic behavior is not reached and a proper model has to be chosen to estimate critical values. We simulated different stochastic models and found that the result from the Gaussian white noise model as given in Eq. ( 11 ) holds already when at least one of the two processes considered can be modeled by Gaussian white noise. However, when none of the processes can be modeled by Gaussian white noise, one has to find suitable models and perform Monte Carlo simulations. Further research is needed to understand the dependency of the critical value on the chosen model.
section_name: Significance testing of wavelet coherency

section_content: 6. 1 Peaks in wavelet cross spectra WCS describes the common power of two processes without a normalization to the single WPS. This can produce misleading results, because one essentially multiplies the CWTs of two time series. E. g. if one of the spectra is locally flat and the other exhibits strong peaks, this can produce peaks in the cross spectrum, which may have nothing to do with any relation of the two time series. Thus, WCS is not suitable for significance testing the relation between two time series. processes are independent by definition, a structure related to the peaks in the NINO3 WPS (i. e. El Niños) occurs (Fig. 5a ). WCO avoids this problem by normalizing to the single WPS (Fig. 5b ). However, if one has detected a significant interrelation using WCO, one can apply WCS to estimate the phase spectrum.
section_name: Pitfalls

section_content: When performing a test on the (1-α)-level, by definition one rejects H 0 with a probability of α% even when it is true. Hence when repeating a test for many independent realizations, around α% of the results should spuriously appear to be significant. This effect is referred to as multiple testing (Lehmann, 1986). Significance testing for WCO can also be seen as an example for multiple testing, since one investigates every point in the space/time domain separately. The important difference to a setting with independent realizations is the fact, that neighboring realizations are correlated according to Appendix A and due to smoothing. 

Figure 6 shows the effect of multiple testing: The 95% level contour line for wavelet coherency between two independent white noise time series is plotted for different smoothing windows. Even though the two processes are independent, peaks covering around 5% of the time/scale area appear. If it is not a priori known if two processes are coherent, the nontrivial issue is now to decide, which of these peaks can be interpreted as significant. A conservative solution is the Bonferroni correction (Lehmann, 1986) dividing the significance level α by the number of tested realizations (i. e. here the number of points in scale and time of the WCO analysis). However, this correction is too strict for correlated realizations. Also, a result based on the Bonferroni correction only states if any coherency exists at all between the two processes without specifying at what scale and time. Thus, deciding if a single peak is significant is not possible based on the Bonferroni correction. If a single peak covers more than 5% of the total area, the question is answered easily. For peaks with smaller area, a valuable hint is given in Fig. 6 : For different smoothing lengths, spurious peaks do not always overlap. Thus one should always compare different smoothing windows and sort out peaks that appear only in one of the corresponding plots. Additionally, it seems that the size of the considered patches plays an important role, as small structures appear much more often than larger ones. One could investigate the distribution of the peak areas under H 0 to draw more distinct conclusions about single peaks. This issue is beyond the scope of this paper.
section_name: Multiple testing

section_content: A prominent example of a bivariate time and scale resolved study is the WCS analysis of the ENSO-NAO teleconnections by Huang et al. (1998), stating that ENSO and NAO are coherent for most of the moderate and strong El Niños between 1900 and 1995. Keeping in mind the discussion in the previous section, we review the WCS results and compare them with those of WCO. We used the extended NINO3index updated from Kaplan et al. (1998) and the NAO-index defined in Jones et al. (1997). We considered N =1738 data points in the range from January 1856 to October 2000 in a monthly resolution. 

Figure 7a shows the results of the WCS analysis (ω 0 =6). As in the analysis of Huang et al. (1998), several peaks appear at a period of about six years around 1877, 1918 and 1940, and at a period between 2 to 4 years around 1888, 1969 to 1982. All these peaks are associated with El Niño Events. However, a significant peak in the WCS might simply arise due to a peak in the univariate WPS of one of the time series as shown in Sect. 6. 1. Thus we additionally calculated the normalized WCO (see Fig. 7b, ω 0 =6, 2w+1=31) and compared it to the WCS. For the calculation of the critical value, we chose Eq. ( 11 ). This critical value matches the result of additional Monte Carlo simulations, using an AR[1] process with an e-folding time of 12 months to model NINO3 and Gaussian white noise to model NAO. For WCO, most peaks from the WCS analysis vanish. However, this analysis shows significant peaks around 1920 at a scale of 4 to 8 years and around 1940 at a scale of about 4 to 14 years. The other small peaks cannot be distinguished from coincidence due to multiple testing. 

These results cast doubt on the finding that ENSO and NAO are coherent for most of the moderate and strong El Niños between 1900 and 1995. However, the period between around 1920 and 1940 seems to be distinct from other times, as a significance coherency occurs. The physical interpretation of this time-dependent interrelation requires further detailed investigation of possible instationary coupling mechanisms., 1888, 1918, 1940, 1969 and 1982. (b) Linear grey scale from black (WCO=0) to white (WCO=1). The black contour line denotes WCO 95. WCO analysis reveals that only a peak around 1920 at a scale of 4 to 8 years and one around 1940 at a scale of about 4 to 14 years are significant.
section_name: Cross wavelet analysis of ENSO/NAO teleconnections

section_content: In this paper we present a detailed evaluation of cross wavelet analysis. We investigate two wavelet measures for the linear scale and time dependent interrelation between two time series, the wavelet cross spectrum (WCS) and the normalized wavelet coherency (WCO). The latter one exhibits values between zero and one, characterizing a vanishing or perfect linear relationship respectively. Based on an inverse problem point of view, we recall that the wavelet periodogram is not a consistent estimator of the local wavelet power spectrum (WPS) and that the inference of peaks from wavelet analysis is not unambiguous. Particularly with regard to WCO, one has to estimate expectation values to get any reliable results. 

In practical cases, this has to be done by smoothing. 

Even for independent mixing processes, realizations will show nonzero WCO-values. Thus, a significance test is re-quired to decide, whether a WCO-value has to be considered as being significantly different from zero. In contrast to Fourier analysis, this test is not independent of the structure of the two processes considered and not easily calculated analytically. Instead, appropriate models have to be chosen and Monte Carlo simulations have to be performed to estimate critical values for a chosen significance level. We exemplified this discussion for the basic case of two processes which can be modeled by Gaussian white noise. The critical value WCO 95 is essentially a function of the smoothing parameters and the time/scale resolution defined by ω 0. An interesting result is that due to correlations of neighboring scales and times, smoothing has only a weak effect. Especially smoothing in time direction is not feasible. For practical purposes, we fitted a polynomial to our numerical results. The critical values derived for the two Gaussian white noise models and thus the approximation Eq. ( 11 ) are valid, if at least one of the processes considered can be modeled by Gaussian white noise. For other cases, appropriate models have to be chosen to perform Monte Carlo simulations. 

A second focus of the paper are pitfalls of cross wavelet analysis. We show that the WCS can show misleading peaks even for realizations of independent processes, just because the WPS of one of the time series exhibits strong peaks. Thus, WCS is not suitable for significance testing the interrelation between two processes and one should rather utilize WCO. Furthermore, due to multiple testing spuriously significant peaks have to appear by definition. A straightforward identification of a single peak is only possible, when this peak covers a percentage of the investigated time/scale area higher than the significance level. For the decision, if a small single peak is significant, no straightforward solution exists so far. However, a valuable hint is given by comparing the peak for different WCO parameters as smoothing length and time/scale resolution. 

On the basis of these results, we could show that the suggested coherency between ENSO and NAO for most of the moderate and strong El Niños between 1856 and 2000 is an artefact. The peaks visible in the WCS rather stem from high power in the single WPS of the NINO3 time series than from any coherency of the two time series. The WCO analysis shows significant peaks only around 1920 at a scale of 4 to 8 years and between 1920 and 1940 at a scale of about 4 to 14 years. 

Summarizing the discussion in this manuscript, we recommend the following steps for cross wavelet analysis:

-Choose suitable models for the significance test of WCO. If at least one of the models is Gaussian white noise, the approximation Eq. ( 11 ) can be used to calculate critical values for a chosen significance level. Otherwise, Monte Carlo simulations have to be performed. 

-Choose the time/scale resolution ω 0 (better compare different values). 

-Calculate WCO and plot a contour line at the critical value. 

-Check the result of different smoothing lengths. If peaks appear only for single smoothing lengths, they are probably artefacts due to multiple testing. 

However, further research is required. Understanding how the critical value for significant WCO depends on the stochastic models used to describe the investigated processes, could help to avoid time intensive Monte Carlo simulations. An important issue is the treatment of multiple testing effects. To our knowledge, there is no standard procedure for a test, that is not too conservative, but does not produce many false positive results. It seems that the area of the investigated peaks in the time/scale domain plays an important role, as small structures appear much more often than larger ones. Thus one could try to estimate the size distribution of spuriously significant patches to draw more distinct conclusions about single peaks.
section_name: Conclusions

section_content: Acknowledgements. We thank M. Kallache, U. Schwarz, M. Rosenblum and A. Pikovsky for inspiring discussions. Our special thanks are dedicated to R. Dahlhaus, who called our attention to the theorem of Ogasawara and Takahashi. D. M. was funded by Deutsche Forschungsgemeinschaft, SFB 555. J. K. was supported by SPP 1114. 

Edited by: M. Thiel Reviewed by: two referees

section_name: Appendix A Correlation of wavelet scales and times

section_content: In Fourier analysis, for a white noise time series neighboring frequencies are uncorrelated. Here we show that this is no longer the case for wavelet analysis. Given a continuous white noise time series x(t), the continuous formulation of Eq. ( 1 ) reads

where two different points in time are uncorrelated:

The correlation between the CWT of two different scales s 1 and s 2 at different times t 1 and t 2 is defined as

Inserting Eq. (A1) into Eq. (A3), we obtain

Since integrals and expectation values are linear operators, this can be simplified to

Using Eq. (A2) and some algebra lead to

Thus, in contrast to Fourier analysis, neighboring points in time and scale are not uncorrelated.
section_name: Appendix A Correlation of wavelet scales and times

section_content: A condition for a simple analytical test statistic for Fourier coherency is based on the fact, that the cross spectrum exhibits a χ 2 distribution (Brockwell and Davis, 1987). Since we have shown that wavelet times and scales are not uncorrelated, this condition is no longer trivially fulfilled. The theorem of Ogasawara and Takahashi (c. f. Rao (1965)) states, when a sum of squares of correlated Gaussian distributed random variables is χ 2 distributed. Given a vector Y of Gaussian random variables

with a covariance matrix, the product Y T AY is χ 2 distributed, when

holds. 

In the case of smoothing neighboring values of the WCS, we are interested in the sum of squares (in this example smoothing over three neighboring scales is considered. Generalization follows straight forward)
section_name: Appendix B Distribution of the cross wavelet spectrum