Wavelet analysis for non-stationary , nonlinear time series

Methods for detecting and quantifying nonlinearities in nonstationary time series are introduced and developed. In particular, higher-order wavelet analysis was applied to an ideal time series and the quasi-biennial oscillation (QBO) time series. Multiple-testing problems inherent in wavelet analysis were addressed by controlling the false discovery rate. A new local autobicoherence spectrum facilitated the detection of local nonlinearities and the quantification of cycle geometry. The local autobicoherence spectrum of the QBO time series showed that the QBO time series contained a mode with a period of 28 months that was phase coupled to a harmonic with a period of 14 months. An additional nonlinearly interacting triad was found among modes with periods of 10, 16 and 26 months. Local biphase spectra determined that the nonlinear interactions were not quadratic and that the effect of the nonlinearities was to produce nonsmoothly varying oscillations. The oscillations were found to be skewed so that negative QBO regimes were preferred, and also asymmetric in the sense that phase transitions between the easterly and westerly phases occurred more rapidly than those from westerly to easterly regimes.


Introduction
Spectral analysis is a tool for extracting embedded structures in a time series.In particular, Fourier analysis has been used extensively by researchers for extracting deterministic structures from time series but is incapable of detecting nonstationary features often present in geophysical time series.Wavelet analysis can extract transient features embedded in time series, with a wavelet power spectrum representing variance (power) of a time series as a function of time and period.Since the seminal work of Torrence and Compo (1998), wavelet analysis has been applied extensively to geophysi-cal time series such as the indices for the North Atlantic Oscillation (Olsen et al., 2012), Arctic Oscillation (Jevrejeva et al., 2003), Pacific Decadal Oscillation (Macdonald and Case, 2005;Newmann et al., 2003), El Niño-Southern Oscillation (ENSO; Torrence and Webster, 1999), Pacific-North American pattern, and west Pacific pattern (Gan et al., 2007).The application of wavelet coherence and cross-wavelet analyses (Grinsted et al., 2004), moreover, has proven useful in relating geophysical time series to other time series (Jevrejeva et al., 2003;Gan et al., 2007;Labat, 2010;Lee and Lwiza, 2008).
Many statistical methods, including power and crossspectral analyses, rely on the assumption that the variable in question is Gaussian distributed (King, 1996).For a linear system in which the output is proportional to the input, the first-and second-order moments (the mean and variance) can fully describe the distribution of a process.In the frequency domain, by analogy, the variable can be fully described by the power spectrum, the decomposition of variance as a function of frequency.Suppose, however, that the distribution is non-Gaussian so that higher-order moments such as skewness and kurtosis exist.In this case, the mean and variance, while useful, are unable to fully describe the distribution in question.In a time series context, non-Gaussian distributions can arise from nonlinear systems, systems for which the output is no longer simply proportional to the input.For a nonlinear system, if the input is the sum of two sinusoids with different frequency components the output will contain additional frequency components representing the sum and difference of the input frequencies (King, 1996).In such cases, it is necessary to examine the decomposition of higher-order moments in frequency space.
The frequency decomposition of the third-order moment, for example, results in a bispectrum or skewness function that measure deviations from Gaussianity (Nikias and Raghuveer, 1987;King, 1996).In fact, Hinich (1985) developed a bispectral test to determine if a time series is non-Gaussian and nonlinear.In some situations, higher-order nonlinearities such as cubic nonlinearities may exist, in which case the trispectrum or other polyspectra would have to be used (Collis et al., 1998).
Another advantage of higher-order spectral analysis is that the cycle geometry of oscillations, such as asymmetry with respect to a horizontal axis (skewed oscillation) or with respect to a vertical axis (asymmetric oscillation) can be quantified using the biphase.A pure sine wave, for example, is neither skewed nor asymmetric, whereas a time series resembling a sawtooth is asymmetric.Skewed and asymmetric cycle geometry can identify, for example, abrupt climatic shifts, sudden shifts in the climate system that exceed the magnitude of the background variability (King, 1996).Abrupt climate shifts have occurred numerous times in the past and have dire impacts on ecological and economic systems (Alley et al., 2003).An understanding of past abrupt climate shifts is essential to understanding future climate change and so there is a need to quantify nonlinearities present in climatic oscillations.
The quasi-biennial oscillation (QBO), as another example, has been shown to behave nonlinearly, transitioning from easterly phases to westerly phases more rapidly than from westerly to easterly phases (Lu et al., 2009).Another source of asymmetry in the QBO time series arises from the westerly shear zone descending more regularly than the easterly shear zone.Asymmetries in the QBO time series are not well captured by linear methods such as linear principal component and singular spectrum analyses (Lu et al., 2009) but are better captured using, for example, nonlinear principal component analysis (Hamilton and Hsieh, 2002).Another example of a nonlinear time series is the sunspot cycle.Solar activity undergoes an 11-year oscillation characterized by asymmetric cycle geometry, with solar maxima generally rising faster than they fall, indicating the presence of nonlinearities (Moussas et al., 2005;Rusu, 2007).ENSO, a climate phenomenon with regional-to global-scale impacts, has also been shown to exhibit nonlinearities (Timmermann, 2003).The presence of nonlinearities and possible nonstationarities in the QBO, ENSO and sunspot time series makes traditional Fourier and wavelet analysis inadequate for feature extraction, underscoring the need to develop methods for quantifying nonlinearities in a nonstationary geophysical setting.
The application of higher-order wavelet analysis has been rather limited compared to traditional wavelet analysis (van Millagan et al., 1995;Elsayed, 2006).One geophysical application of higher-order wavelet analysis is to oceanic waves (Elsayed, 2006), which was found to be capable of identifying nonlinearities in wind-wave interactions.However, the study lacked rigorous statistical significance testing, which is problematic because even a Gaussian process of finite length can produce nonzero bicoherence.Therefore, the first aspect of this paper is to apply significance testing methods for higher-order wavelet analysis to aid physical interpretation of results.
The number of bicoherence estimates to which the statistical test is applied will be large and multiple artifacts will result.The multiple-testing problem was already identified for traditional wavelet analysis (Maraun et al., 2007;Schulte et al., 2015;Schulte, 2016).The first objective of this paper will be therefore to apply statistical methods controlling false positive detection.It is also noted that the bicoherence spectra calculated are only sample estimates of the true bicoherence spectra.The second objective of this paper will be to develop a procedure for calculating confidence intervals corresponding to the sample estimates, which represent a range of plausible values for the sample estimates.
Another problem with the application of higher-order wavelet analysis is selection of a time interval on which to calculate the high-order wavelet quantities.Such an approach is subjective and the result of the analysis may depend on the time interval chosen.Objective three of this paper will address the time interval selection problem.Such an approach has already been adopted in wavelet coherence analysis (Grinsted et al., 2004).
Additionally, properties of the biphase have only been examined for Fourier-based bispectral analysis (Elgar and Sebert, 1989;Maccarone, 2013) and its usefulness in higherorder wavelet analysis has yet to be examined.For nonstationary time series, the biphase and cycle geometry corresponding to the time series may change with time and thus objective four of this paper will be to introduce a local wavelet-based biphase spectrum.
In this paper, higher-order wavelet analysis is put in a statistical framework and applied to the QBO time series to demonstrate the insights afforded by the methods.Before describing higher-wavelet analysis, a brief overview of wavelet analysis is first presented in Sect. 2. Higher-order wavelet analysis is described in Sect. 3 and a new local autobicoherence spectrum is introduced, eliminating the selection of a time interval on which to calculate nonlinear properties of time series.The new and existing methods are applied to an ideal time series and the QBO index.In Sect.4, a new procedure for estimating confidence intervals of global autobicoherence quantities is developed to estimate uncertainties in the sample autobicoherence spectra.The application of the new procedure to the sample autobicoherence spectrum of the QBO time series is then used to further assess confidence in results.

Wavelet analysis
The idea behind wavelet analysis is to convolve a time series with a function satisfying certain conditions.Such functions are called wavelets, of which the most widely used is the Morlet wavelet, a sinusoid damped by a Gaussian envelope: where ψ 0 is the Morlet wavelet, ω 0 is the dimensionless frequency and η is the dimensionless time (Torrence and Compo, 1998;Grinsted et al., 2004).In practical applications, the convolution of the wavelet function with a time series X = (x n ; n = 1, . .., N) is calculated discretely using where δt is a uniform time step, s is scale, η = s × t and W X n (s) is the wavelet transform.The wavelet power is given by W X n (s) 2 (Torrence and Compo, 1998;Grinsted et al., 2004).For the Morlet wavelet with ω 0 = 6, the wavelet scale and the Fourier period λ are approximately equal (λ = 1.03s).A more detailed discussion of wavelet analysis can be found in Torrence and Compo (1998).Shown in Fig. 1a is the time series of the QBO index and shown in Fig. 1b is the corresponding wavelet power spectrum.The QBO data from 1950 to 2013 were obtained from the National Oceanic Atmospheric Administration Earth System Research Laboratory (available at: http:// www.esrl.noaa.gov/psd/data/climateindices/list/).The QBO index is defined as the zonal average of the 30 hPa zonal wind at the Equator.As such, a positive index indicates westerly winds and a negative index indicates easterly winds.The most salient feature of the time series is the rather regular periodicity of approximately 28 months.Also note the asymmetry between the negative and positive phase, with the negative phases generally being stronger.The periodic behavior of the QBO was corroborated by examining the wavelet power spectrum.A well-defined 28-month periodicity is evident, with the associated wavelet power changing little throughout the study period.
There are also secondary features located at a period of approximately 14 months, primarily from 1985 to 2013.The appearance of significant power at a period of 14 months also coincides with most of the largest negative phases of the QBO.Such a correspondence may not have been a coincidence; the 14-month mode and the 28-month mode may have interacted constructively to generate large negative events but interacted destructively to create smaller positive events.However, additional tools are needed to confirm if the periodicities are interacting and to understand how the interactions were related to the behavior of the QBO.
3 Higher-order wavelet analysis

Wavelet-based autobicoherence
Higher-order spectral analysis provides the opportunity to quantify nonlinearities and allows for the detection of inter-  (Torrence and Compo, 1998).Light shading represents the cone of influence; i.e., the region in which edge effects cannot be ignored.acting oscillatory modes within a time series.More specifically, nonlinearities are quantified using bicoherence, a tool for measuring quadratic nonlinearities, where quadratic nonlinearities imply that for frequencies f 1 , f 2 and f 3 and corresponding phases φ 1 , φ 2 and φ 3 the sum rules and are satisfied; Eq. ( 3) implies frequency coupling and Eq. ( 4) implies phase coupling.To see from where Eqs. ( 3) and ( 4) originate let be the input into a system, whose output is related to the input by The multiplicative factor ε is used to represent the contribution of the nonlinear component of the signal and w(t) is Gaussian white noise.Note that if ε = 0, then the system is linear because the output contains the same frequency components as the input.The substitution of Eq. ( 5) into Eq.( 6) results in and thus the output has sinusoids with additional frequency components 2f 1 , 2f 2 , f 2 − f 1 and f 2 + f 1 , which arise from the second term in right-hand side of Eq. ( 6).Unlike the power spectrum, which is the Fourier transform of the second-order moment of a time series, the bispectrum is defined as the double Fourier transform of the third-order moment, or, more generally, the third-order cumulant, i.e., where C is the third-order cumulant, defined as and the t i are lags.If X(t) is zero-mean, then in Eq. ( 9), M 1 = E [X(t)] = 0 denotes the first-order moment (mean), M 2 = E [X(t)X(t + t 1 )] denotes the second-order moment (autocorrelation) and denotes the third-order moment (Nidal and Malik, 2013).Also note that for a zero-mean process, the third-order cumulant reduces to the third-order moment (Collis et al., 1998).A more useful quantity is the normalized version of the bispectrum, the autobicoherence spectrum (Collis et al., 1998), which can be computed using the following: where b 2 (f 1 , f 2 ) is bounded by 0 and 1 by the Schwarz inequality and X f denotes the Fourier transform of X. b 2 (f 1 , f 2 ) can be interpreted as the fraction of power at f 1 + f 2 due to quadratic phase coupling among f 1 , f 2 and gar and Chandran, 1993).For a more in-depth discussion of higher-order spectral analysis the reader is referred to Nikias and Raghuveer (1987).Phase information and cycle geometry can be obtained from the biphase, which is given by It was noted by Maccarone (2013), however, that the biphase should be defined on the full 2π interval and thus in this paper the four-quadrant inverse tangent is computed and not the inverse tangent as shown above.By doing so, statistically significant autobicoherence detected together with the biphase can be used to quantify cycle geometry.A biphase of 0 • indicates positive skewness and a biphase of 180 • indicates negative skewness (Maccarone, 2013).An example of a skewed oscillation time series with biphase close to 0 • is shown in Fig. 2a.Mathematically, the time series is written as where a = 0 (Maccarone, 2013).The time series is skewed because the positive spikes are not accompanied by negative spikes of equivalent magnitude and therefore the distribution of the time series would be positively skewed, with the right tail being larger than the left tail.
For asymmetric waveforms, a biphase of 90 • indicates that the time series is linearly rising but rapidly falling as shown in Fig. 3, whereas a biphase of −90 • indicates that the time series rises rapidly and falls linearly.A purely asymmetric time series will have a biphase of 90 • or −90 • , as shown in Fig. 3, where the sawtoothed time series obtained by setting a = π/2 in Eq. ( 12) rises more slowly than it falls.In a physical setting, asymmetric cycle geometry implies that phase transitions occur at different rates, as observed in the QBO time series.
According to Elsayed (2006), the wavelet-based autobicoherence is defined as where T is a time interval, W x (s, t) is the wavelet transform of a time series X at scale s and time t, and W * x (s, t) denotes the complex conjugate of W x (s, t).The wavelet-based autobicoherence measures the degree of quadratic phase coupling, where a peak at (s 1 , s 2 ) indicates a statistical dependence among the scale components s 1 , s 2 and s.
In practice, the autobicoherence is computed discretely so that Eq. ( 16) can be written as where n 1 ≥ 1, and n 2 ≤ N. Note that if n 1 = 1 and n 2 = N, then Eq. ( 16) represents the global autobicoherence spectrum.
The Monte Carlo approach to pointwise significance testing is adopted in this paper and is similar to that used in wavelet coherence (Grinsted et al., 2004).To estimate the significance of wavelet-based autobicoherence at each point (s 1 , s 2 ), Monte Carlo methods are used to (1) generate a large ensemble of red-noise processes with the same lengths and lag-1 autocorrelation coefficients as the input time series and (2) compute for each randomly generated red-noise process the autobicoherence spectrum.From the ensemble of autobicoherence spectra, the p = 100(1 − α p ) percentile of the autobicoherence estimates is computed for every point (s 1 , s 2 ), where p corresponds to the critical level of the test and α p is the pointwise significance level of the test.Given the symmetry of the autobicoherence spectrum, the critical level of the test can be computed using only half of the autobicoherence estimates, reducing computational costs.

Multiple testing
Let α p be the significance level of the pointwise significance test as described above and let K denote the number of autobicoherence estimates being tested, then there will be on average α p K false positive results.A similar problem occurs in traditional wavelet analysis (Maraun et al., 2007;Schulte et al., 2015;Schulte, 2016).In the case of simultaneously testing multiple hypotheses, the number of false positive results can be reduced by applying, for example, the Bonferonni correction (Lehmann, 1986).However, this simple correction often results in many true positives being rejected and is especially permissive in the case of autocorrelated data (Maraun et al., 2004).Other procedures also exist, including the Walker p value adjustment procedure, which has more statistical power than the Bonferonni correction.An even more powerful method is the Benjamini and Hochberg (1995) procedure, which controls the false discovery rate (FDR), where the FDR is the expected proportion of the false rejections that are actually true.An advantage of this method, in addition to its statistical power, is that it takes into account the confidence with which local hypotheses are rejected and is robust even in the case of autocorrelated data (Wilks, 1997).Benjamini and Yekutieli ( 2001) developed a modified version of the Benjmini and Hochberg (1995) procedure that works for any dependency structure among the local test statistics and thus this procedure will be used in this paper to control the FDR.
The procedure can be described as follows: suppose that K local hypotheses were tested.Let p (i) denote the smallest of the K local p values, then, under the assumption that the K local tests are independent, the FDR can be controlled at the q level by rejecting those local tests for which p (i) is no greater than so that the FDR level is equivalent to the global test level α global .For a local p value to be deemed significant using this procedure, it must be less than or equal to the largest p value for which Eq. ( 19) is satisfied.If no such local p values exist, then none are deemed insignificant, and, therefore, the global test hypothesis cannot be rejected.If the test statistics have an unknown dependency structure, q can be replaced with q/ K i=1 1 i , though this substitution makes the procedure less powerful (Groppe et al., 2011).This modified method will be applied to autobicoherence spectra at the 0.05 level throughout this paper.

Wavelet-based autobicoherence of an idealized time series
To demonstrate the features of a time series that can be extracted using higher-order wavelet analysis, an idealized non- stationary time series will first be considered.Consider the quadratically nonlinear time series where f is frequency, w(t) is Gaussian white noise and γ (t) is a time-dependent nonlinear coefficient given by Note that Eqs. ( 3) and ( 4) are satisfied because f 1 + f 2 = 2f 1 = 2f 2 and similarly for φ.The sinusoid with frequency 2f 1 is said to be the harmonic of the primary frequency component with frequency f 2 , where the amplitude of the harmonic depends on γ (t), the strength of the quadratic nonlinearity.X(t) and the corresponding wavelet power spectrum for the case when f 1 = 0.03 is shown in Fig. 4. The signal-tonoise ratio of the Gaussian white noise was set to 1 decibel.The primary frequency component results in a large region of 5 % pointwise significance at λ = 30, whereas its harmonic only results in a few small significance regions located from t = 700 to t = 1000.It also noted that the appearance of the significance power at λ = 15 from t = 700 to t = 1000 is accompanied by large positive spikes in the time series that result in the time series favoring positive values.Prior to the emergence of the significant power at λ = 15, the time series varied smoothly in the sense that negative phases were accompanied by positive phases of similar amplitude.
To determine if the oscillations are quadratically interacting, the autobicoherence of X (t) was computed (Fig. 5).The significant peak centered at (30,30) indicates that an oscillation with period 30 is phase coupled to an oscillation with λ = 15.The result implies that the variability at λ = 15 is partially related to the statistical dependence between the two modes.The fraction of variability is determined by the autobicoherence value corresponding to the significant peak.In the present case, W b (s 1 , s 2 ) = 0.5 so about half of the variability at λ = 15 is due to the nonlinear interaction.Note that no other peaks were found to be significant.

Wavelet-based autobicoherence of geophysical time series
Shown in Fig. 6 is the wavelet-based autobicoherence spectrum for the QBO time series.A large region of significance was identified, which contained the local maximum at (28, 28) months.The peak represents the phase coupling of the primary frequency component with its harmonic with a period of 14 months.The power at λ = 14 months therefore is partially related to the statistical dependence between its primary frequency component and its harmonic.The significance and magnitude of the autobicoherence in the QBO spectrum is consistent with how the QBO does not vary smoothly, shifting to the easterly phase more quickly than to the westerly phase and with the westerly phase tending to be stronger than the easterly phase.The asymmetry in both phase transition and magnitude are suggestive of nonlinearities.

Local wavelet autobicoherence
It may also be desirable to see how autobicoherence along slices of the full autobicoherence spectrum changes with time.To compute local autobicoherence, apply a smoothing operator S(W ) = S scale S time W X n (s) (Grinsted et al., 2004) to each term in Eq. ( 13) instead of summing in time, i.e." The smoothing operator for the Morlet wavelet is given by and where c 1 and c 2 are normalization constants determined numerically and is the rectangular function.
It is important to mention that the numerator of Eq. ( 22) contains a term with wavelet coefficients at two different scales so that the choice of smoothing is not as straightforward as for wavelet coherence.Smoothing autobicoherence estimates with respect to s min = min(s 1 , s 2 ) was found to result in larger autobicoherence estimates, whereas smoothing the autobicoherence with respect to s max = max(s 1 , s 2 ) resulted in smaller autobicoherence estimates.Given that the autobicoherence estimates are influenced by the choice of smoothing, it is inevitable that the significance of the autobicoherence estimates is also impacted.In particular, smoothing the autobicoherence spectrum with respect to s max allowed extrema to be smoothed out, eliminating spuriously large autobicoherence.For this reason, all local autobicoherence spectra in this paper will be computed by smoothing with respect to s max .
The advantage of using Eq. ( 22) is that transient quadratic nonlinearities can now be detected and the need for choosing an integration time interval has been eliminated.If s 1 = s 2 , then (t, s 1 , s 1 ) = (t, s 2 , s 2 ) = (t, s) and thus, in the case of this diagonal slice, the local wavelet-based bicoherence spectrum is a two-dimensional representation of the degree of local quadratic nonlinearity.The vertical axis corresponds to the primary frequency and the horizontal axis corresponds to time.As a concrete example, a peak at (64, 64) would indicate that at time index t = 50 the oscillation with a fundamental period λ = 1.03s ≈ 64 is locally coupled to an oscillation with period λ ≈ 32.
One can also compute a local biphase from the smoothed bispectrum by taking the four-quadrant inverse tangent of the smoothed imaginary part divided by the smoothed real part.The local biphase, for example, was computed for the skewed time series shown in Fig. 2a.As expected, the biphase fluctuates regularly around 0 • and the mean is 2 • .The local biphase for the sawtoothed time series is shown in Fig. 3b.The biphase fluctuates about 90 • and the mean biphase is 90 • as expected.
The procedure for the estimation of the statistical significance of local autobicoherence is the following: generate red-noise time series with the same lag-1 autocorrelation coefficients as the input time series and use the local autobicoherence estimates outside the COI to generate a null distribution of b w n (s 1 , s 2 ).Note that the calculation only needs to be performed at a fixed time outside of the COI because red-noise is a stationary process, which produces a stationary background spectrum.

Local wavelet-based autobicoherence of an idealized time series
The local autobicoherence spectrum of X (t) for (30, 30) is shown in Fig. 6b.Initially, there is no local autobicoherence that exceeds the 5 % significance level.At t = 250 and t = 500, on the other hand, small regions of 5 % significant autobicoherence emerge, indicating a transient nonlinear interaction.At t = 500 the nonlinearity is strong and results in a large region of significant local autobicoherence extending from t = 500 to the edge of the wavelet domain In order to determine if the peaks in autobicoherence are associated with a quadratic nonlinearity, it is important to compute the biphase, which is shown in Fig. 7b.From t = 0 to t = 400 there is an unstable phase relationships between the phase of the primary frequency component and its harmonic.Such a lack of phase coherence indicates a weak nonlinear interaction, which is consistent with how the autobicoherence is lower before t = 400.In contrast, after t = 400, the biphase becomes stable, changing little with time, indicating a consistent phase relationship between the primary frequency mode and its harmonic.It also noted that the biphase during this time fluctuates near 0 • , which implies that the phase relationships arise from a quadratic nonlinearity.The near-zero biphase is consistent with how X(t) was constructed from the sum of two cosines with zero phase and also suggests that the interaction results in skewed cycle geometry, where positive values of the time series are preferred.Indeed, by inspection of Fig. 4a the oscillations initially appear to be sinusoidal, varying smoothly, whereas after t = 400 spikes begin to appear and X (t) favors positive values.

Local wavelet-based autobicoherence of the QBO time series
The local autobicoherence spectrum of the QBO index at the point (28, 28) in the full autobicoherence spectrum is shown in Fig. 8. From 1950 to 1970 the magnitude of the autobicoherence fluctuated and consisted of one local significant peak at 1965.Significant autobicoherence was also found from 1975 to 1998, contrasting with the autobicoherence after 1998, which was not found to be significant until 2010.
To determine if the peaks indicated in the autobicoherence spectrum are associated with a quadratic nonlinearity, the local biphase was computed.Figure 8a shows the local biphase for the autobicoherence peak at (28,28).For most of the study period, the biphase was found to vary considerably, particularly during the 1950-1970 and 1995-2013 periods.On the other hand, the biphase varied smoothly from 1970 to 1995, consistent with how the autobicoherence during that period was large and stable (Fig. 8a).Also, during that period the biphase was nonzero; in fact, the mean biphase during the period was −100 • , suggesting that the phase coupling is not the result of a quadratic interaction.A biphase of −100 • indicated asymmetric geometry, which physically represents how phase transitions of the QBO occurred at different rates.Recall that it has already been discussed in the introduction that the QBO transitions from easterly phases to westerly phases more rapidly than from westerly to easterly phases (Lu et al., 2009).Another interesting feature is the general increase in the biphase from 1970 to 1995.In the beginning of the time period, the biphase was −180 • and after 1980 the biphase switched to −90 • .
The local autobicoherence and biphase corresponding to the peak (16, 26) was also computed (Fig. 9).The mean of the  absolute value of the biphase for the period 1950-2013 was 130 • , indicating a statistical dependency among the modes with periods of 10, 16 and 26 months resulted in skewed waveforms.In fact, because the biphases were close to 180 • the waveforms should have been skewed to negative values (Maccarone, 2013) and such skewness is evident by inspecting Fig. 1.Also note that some of the largest negative phases of the QBO occurred from 1995 to 2010, which coincided with the period of most significant autobicoherence as shown in Fig. 9a.

Block bootstrapping autobicoherence
Bootstrapping is a widely used technique to estimate the variance or uncertainty of a sample estimate.For independent data, one samples with individual replacement data points (Efron, 1979); for dependent data one must sample with replacement blocks of data to preserve the autocorrelation structure of the data (Kunsch, 1989).The latter technique is called block bootstrapping and should be used for variance estimation of global wavelet quantities, as wavelet coefficients are known to be autocorrelated in both time and scale.The use of traditional bootstrapping techniques would result in confidence intervals that are too narrow.It is expected, however, that the choice of the bootstrapping technique is more critical at larger scales, as the decorrelation length of the mother wavelet increases with scale.
A brief overview of the procedure is provided below but a more detailed discussion can be found in Schulte et al. (2016).To find the approximate 100(1 − β) % confidence interval of an autobicoherence estimate, divide the set of wavelet coefficients at each scale into overlapping blocks.The lengths of the blocks at each scale should be the same and the randomly resampled blocks chosen should be the same at each scale to avoid randomizing the data.The concatenation of the blocks then results in a synthetic set of wavelet coefficients at each scale.The synthetic set of wavelet coefficients can then be used to calculate a bootstrap replicate of the autobicoherence.The iteration of the procedure 1000 times results in a distribution of bootstrap replicates from which a 95 % confidence interval can be obtained.
As noted by Schulte et al. (2016), the appropriate block length to use can be determined by Monte Carlo methods.In that study, it was determined from a Monte Carlo experiment that a block length of N 0.6 produced accurate confidence bounds for wavelet coherence while also producing the widest confidence intervals at all scales.The Monte Carlo experiment was repeated for 95 % confidence in this study because bicoherence estimation requires the use of wavelet coefficients at three wavelet scales, with the wavelet coefficients at each scale having a different correlation structure.For wavelet coherence, the block length selection procedure is simpler because a single wavelet scale is used so that correlation structure of wavelet coefficients is similar.The Monte Carlo experiment was performed by generating rednoise processes of length 1000 with different lag-1 autocorrelation coefficients and computing 95 % confidence inter-vals around the estimated autobicoherence.Remarkably, the Monte Carlo experiment found that a block length of N 0.6 is also optimal for bicoherence confidence interval estimation.For block lengths exceeding N 0.6 , confidence intervals were found to be too narrow, with in some instances the estimated bicoherence falling outside the 95 % confidence interval.It is also noted that the results were insensitive to the chosen lag-1 autocorrelation coefficient.

Application to ideal and climatic time series
Figure 5b shows the application of the block bootstrap procedure to the diagonal slice s 1 = s 2 = s of the autobicoherence for the ideal case.The 95 % confidence intervals were also obtained using the ordinary bootstrap method.A pronounced peak at s = 30 was identified and represents the interaction between the primary frequency and its harmonic.By inspection of Fig. 5b, there is a clear difference between the widths of the confidence intervals obtained from the two bootstrapping procedures.For the ordinary bootstrap, the confidence intervals are narrow and the widths of the confidence intervals appear to be only weakly dependent on scale.On the other hand, the confidence intervals obtained using the block bootstrap procedure are wide, especially at large scales, and the width of the confidence intervals depends strongly on scale, increasing from small scales to large scales.It is also noted that, whereas the block bootstrap procedure has deemed no spurious peaks as significant, the ordinary bootstrap procedure deemed two the spurious peaks at s = 14 and s = 100 as significant.The implementation of the block bootstrap procedure can therefore enhance confidence in results, facilitating the investigation of a deeper physical understanding.
The application of the block bootstrap procedure to the diagonal slice s 1 = s 2 = s of the full autobicoherence spectrum of the QBO index is shown in Fig. 10.The 95% confidence intervals corresponding to the peaks (14, 14) and (28, 28) do not cross the 5 % significance bound and thus one has more confidence that those peaks are significant.All other peaks have been deemed insignificant.

Summary
Higher-order wavelet analysis together with significance testing procedures were used to detect nonlinearities embedded in an ideal time series and the QBO time series.The autobicoherence spectrum of the QBO index revealed phase coupling of the 28-month mode with a higher-frequency mode with period 14 months.A local autobicoherence spectrum of the QBO index showed that the strength of the nonlinearities varied temporally.Furthermore, the local biphase spectrum indicated that a statistical dependence among frequency components resulted in waveforms that were both skewed and asymmetric, indicating that the strength of negative QBO events were stronger than positive events, and that transitions between events occurred at different rates.The author has written a software package (Schulte, 2015) to implement all higher-order wavelet analysis methods presented in the paper.

Data availability
The data for the Quasi-biennial Oscillation used in this paper are freely available at the National Oceanic Atmospheric Administration's Earth System Research Laboratory Physical Science Division website available at: http://www.esrl.noaa.gov/psd/data/climateindices/list/.

Figure 1 .
Figure 1.(a) The QBO index and (b) the corresponding wavelet power spectrum.Contours enclose regions of 5 % statistical pointwise significance(Torrence and Compo, 1998).Light shading represents the cone of influence; i.e., the region in which edge effects cannot be ignored.

Figure 2 .
Figure 2. (a) a skewed time series and (b) its corresponding local biphase.The biphase close to zero indicates a nonlinear interaction resulting in a skewed oscillation.The biphase was calculated from the first three cosines in the summation described in the text.

Figure 3 .
Figure 3. (a) A sawtoothed time series and (b) its corresponding local biphase.The biphase close to 90 • indicates a nonlinear interaction resulting in an asymmetric waveform.The biphase was calculated from the first three cosines in the summation.

Figure 5 .
Figure 5. (a) Wavelet-based autobicoherence spectrum of the ideal time series.Thick contours enclose regions of 5 % pointwise significance after controlling the FDR.The diagonal line separates the spectrum into two symmetric regions.(b) The diagonal slice of the autobicoherence spectrum at s 1 = s 2 = s.The critical level for the test represented by the dotted line was calculated using Monte Carlo methods.

Figure 6 .
Figure 6.The wavelet-based autobicoherence spectrum of the QBO index for the period 1950-2013.Thick contours enclose regions of 5 % pointwise significance.

Figure 7 .
Figure 7. (a) The local autobicoherence and (b) local biphase corresponding to (30, 30) in the full autobicoherence spectrum shown in Fig. 5a.Biphases differing from 90 • indicate that the nonlinear interaction resulted in a waveform with skewness.

Figure 8 .
Figure 8. Same as Fig. 7 except at (28, 28) in the autobicoherence spectrum of the QBO index.Biphases differing from 90 • indicate that the nonlinear interaction resulted in a waveform with skewness.

Figure 10 .
Figure 10.Same as Fig. 5b except for the QBO index for the period 1950-2013.