Some Aspects of the Discrete Wavelet Analysis of Bivariate Spectra for Business Cycle Synchronisation

The paper considers some of the issues emerging from the discrete wavelet analysis of popular bivariate spectral quantities such as the coherence and phase spectra and the frequency-dependent time delay. The approach utilised here is based on the maximal overlap discrete Hilbert wavelet transform (MODHWT). Firstly, via a broad set of simulation experiments, we examine the small and large sample properties of two wavelet estimators of the scale-dependent time delay. The estimators are the wavelet cross-correlator and the wavelet phase angle-based estimator. Our results provide some practical guidelines for the empirical examination of short- and medium-term lead-lag relations for octave frequency bands. Further, we point out a deficiency in the implementation of the MODHWT and suggest using a modified implementation scheme, which was proposed earlier in the context of the dual-tree complex wavelet transform. In addition, we show how MODHWT-based wavelet quantities can serve to approximate the Fourier bivariate spectra and discuss issues connected with building confidence intervals for them. The discrete wavelet analysis of coherence and phase angle is illustrated with a scale-dependent examination of business cycle synchronisation between 11 euro zone countries. The study is supplemented by a wavelet analysis of the variance and covariance of the euro zone business cycles. The empirical examination underlines the good localisation properties and high computational efficiency of the wavelet transformations applied and provides new arguments in favour of the endogeneity hypothesis of the optimum currency area criteria as well as the wavelet evidence on dating the Great Moderation in the euro zone.


Introduction
Wavelet analysis is a kind of frequency study that enables us to examine local signal properties efficiently. It is a relatively new mathematical construct with a broad range of applications in statistics, data compression and image processing. But this approach has also found its place in modern time series analysis as it makes it possible to analyse time series that are subject to structural breaks, outliers, local trends and changing cyclical patterns or those that show other transient characteristics. The distinguishing feature of this technique among other time-frequency methods is an endogenously varying time window, i.e. the ability to analyse short oscillations with narrow time windows and longer cycles with wider windows. Thus, wavelet methodology is thought to constitute the next logical step in spectral analysis, one that elaborates on the time localisation properties of frequency methods. The methodology is known to have a significant impact on, e.g., geophysics, oceanography and medicine. However, it is much less popular among economic sciences, with business cycle studies becoming one of for all (discrete) time units and octave frequency bands (or dyadic scales). The MODHWT-based methodology encompasses the analysis performed with the usual maximal overlap discrete wavelet transform (MODWT) and goes beyond it producing quantities which directly correspond to the bivariate Fourier spectra. In our theoretical considerations we start with pointing out the need to modify the implementation scheme of the MODHWT in a way suggested in the context of the dual-tree complex wavelet transform of Kingsbury and Selesnick (see, e.g., Selesnick et al. 2005). Further, we show how the wavelet quantities can serve to approximate the Fourier cross-spectra and discuss issues connected with constructing confidence intervals (CIs) for them. Finally, via a broad set of simulation experiments we examine the small and large sample properties of two wavelet estimators of the scale-dependent time delay -a quantity measuring the causal distance between time series on a scale-by-scale basis. The estimators are: the wavelet cross-correlator (WCC) and the wavelet phase angle (WPA)-based estimator. Our results provide some practical guidelines for the empirical examination of short-and medium-term lead-lag relations for octave frequency bands, pointing out the better small sample performance of the WPA-based estimator, especially in the case of low signal-to-noise ratios (SNRs). In the empirical analysis of business cycle variability and synchronisation in the euro zone we utilize the two 'continuous discrete' wavelet transforms: the maximal overlap DWT and the MODHWT. From the point of view of an economist willing to study business cycles, the MODWT and MODHWT offer the following: − A model-free (non-parametric) approach to examining the frequency characteristics of economic processes, i.e. short-, medium-and long-run features in the series. In particular, because of their non-parametric nature, wavelets enable us to examine non-linear processes without the loss of information. − Good time-frequency resolution, and because of this, efficiency in terms of the computations needed to extract the features. This enables the precise examination of a time-varying frequency content of time series in an efficient way.
− The decomposition of the variance and covariance of stationary processes according to octave frequency bands. 2 In particular, the wavelet variance gives a simplified alternative to the spectral density function, which uses just one value per octave frequency band. The same is true for the wavelet co-and quadrature spectra, which give piecewise constant approximations to the appropriate Fourier cross-spectra on a scale-by-scale basis (see section 2.4). − Precise timing of the shocks causing and influencing business cycles. − Low computational complexity. 3 − The examination of trended, seasonal and integrated time series without prior transformations . In particular, we do not need to deseasonalise the data, as seasonal components are left automatically in lower decomposition levels, unless one is interested in examining very short cycles less than two years in length. In addition, there is no need for the prior elimination of deterministic and stochastic trends because wavelet filtering usually embeds enough differencing operations. − The efficient estimation of short-term lead-lag relations for different frequency bands. − Global and local (short-term) measures of association for business cycle components such as wavelet correlations and cross-correlations, wavelet coherence and wavelet phase angle.
Recent studies on business cycle synchronisation within the euro zone (see, e.g., de Haan et al. 2008;Gonçalves et al. 2009; and references therein) usually provide evidence in favour of the endogeneity hypothesis of the optimum currency area criteria as stated in Frankel and Rose (1998), according to which (intraindustry) trade intensification and monetary integration lead to more correlated business cycles. Our empirical examination covering 11 euro zone member countries tries to contribute to the debate by looking at synchronisation patterns _________________________ 2 See Percival (1995) for the variance. The covariance case is examined by Whitcher (1998) (see also Whitcher et al. 2000).
3 The conventional discrete wavelet transform (DWT) can be computed using an algorithm that is faster than is the well-known fast Fourier transform (FFT) -Mallat's pyramid algorithm based on a mirror filters cascade and downsampling by two, which requires only O(N) multiplications. By contrast, the computational complexity of the MODWT is O(Nlog 2 N), which is exactly the same as the FFT (Percival and Walden 2000: 159), while the MODHWT consumes twice as many operations.
alone, decomposed on a scale-by-scale basis. The study documents a rise in synchronisation between business cycles after the first steps towards European integration were taken in the second half of the 1980s. In addition, changes in business cycle variability are examined, thus providing new evidence on dating the Great Moderation and staying in agreement with the hypothesis of an early start of the process (Blanchard and Simon 2001).
The structure of the paper is as follows. In the next section, we briefly introduce the wavelet transform in its conventional and non-decimated (maximal overlap) versions as well as the wavelet analysis of variance and covariance. Next, we present the MODHWT and discuss more deeply the wavelet bivariate spectral analysis, its connections with the Fourier analysis as well as certain implementation problems. Section 3 presents the results of the simulation analysis comparing two wavelet methods of examining lead-lag relations for octave frequency bands, while Section 4 summarises our empirical findings. Finally, the last section offers brief conclusions.

Wavelet Analysis
Wavelet transformation consists in decomposing a signal into the shifted and scaled versions of a basic function ) (x ψ , called the mother wavelet. There are different kinds of this decomposition depending on the wavelet transform applied. The continuous wavelet transform (CWT) enables us to recognise local features in the data, especially in the case of signals defined over the entire real axis, although it results in an excessive redundancy of information. The DWT provides a parsimonious representation of the data and is particularly useful in discrete time series processing, especially in noise reduction and information compression. The MODWT removes certain deficiencies of the conventional discrete transformation by considering all time units, while -similarly to the DWT -octave frequency bands are analysed.

Conventional and Maximal Overlap DWTs
The conventional DWT of a real-valued function is defined as follows: where , and the wavelet daughters are dyadically shifted and scaled versions of the mother wavelet, i.e.: is usually defined via another function (the scaling function or father wavelet), ) (x φ that, applied to the signal after shifting and scaling analogous to (2), produces another set of coefficients in the form: These are known as scaling coefficients. For a given j, the wavelet coefficients are computed as differences in moving averages of the previous scaling coefficients and are associated with scale , while their squares contribute to the decomposition of the energy of the signal on the time-frequency plane. By contrast, the level j scaling coefficients are moving averages of scale .
These two types of coefficients give the multiresolution decomposition of the original function in the form: The functions and are known as approximations (smooths) and details. The highest level approximation represents the smooth, lowfrequency component of the signal, while the details , , …, are associated with oscillations of length In filtering notation, the DWT is defined via quadrature mirror filters: the lowpass (scaling) filter and the high-pass (wavelet) filter .
When processing discrete signals we consider a vector of length in the form . Then the highest possible decomposition level is J and the number of wavelet and scaling coefficients of the conventional DWT for level is , respectively. By contrast, the MODWT produces the same number of wavelet and scaling coefficients at each decomposition level ( t j W , and t j V , , accordingly) as it does not use downsampling by two. The coefficients are appropriately scaled in order to retain variance preservation. They are given as follows: where and are the j-th level wavelet and scaling filters of length obtained by convolving together the following j filters (Percival and Walden 2000: 95-96):  g  j   g  ,  ,  ,  ,  g  ,  ,  ,  ,  ,  g  ,  ,  ,  ,  g   g  ,  ,  g  ,  ,  ,  g  ,  , is a low-pass filter with the cut-off frequency . In the notation above we assume that . We further use also: For further considerations, we provide relationships between the transfer functions of the above filters with the following correspondence: Then, we have (Percival and Walden 2000: 154): The same relationships also hold for the transfer functions of , which we further denote as Among the most popular real wavelet and scaling filters are the compactly supported orthonormal Daubechies filters: the extremal phase (dL) and least asymmetric (laL) filters. The two families are characterised by the smallest filter length L for a given number of vanishing moments (VMs). 5 In addition, the extremal phase scaling filters have the fastest build-up of the energy sequence, while the least asymmetric filters are approximately linear phase.

Wavelet Analysis of Variance and Covariance
For the stochastic process , the time-dependent wavelet variance is defined as: Assuming that (12) does not depend on time, 6 we arrive at the variance decomposition across scales in the form: _________________________ 5 Roughly speaking, VMs are responsible for eliminating artefacts because of the wavelet function as well as for the degree of approximation to an ideal bandpass filter, and they make it possible to interpret the filters as generalised differences of adjacent observations with the number of embedded difference operations equal to the number of VMs (see Daubechies 1992: 153;Mallat 1998: 166;Percival and Walden 2000: 483). The number of VMs for the Daubechies filters equals half the filter length. 6 The assumption is also fulfilled for non-stationary processes if they are integrated of order d and the width of the wavelet filter, L, is sufficient to eliminate non-stationarity. In the case of Daubechies wavelet filters, the condition is: L ≥ 2d (see, e.g., Percival and Walden 2000: 304). In order to have , we further assume that L > 2d.

www.economics-ejournal.org
The wavelet variance at level j, corresponding to scale , , explains the variation of oscillations of a length approximately in the interval . Similarly, the wavelet covariance and wavelet correlation are introduced. For the stochastic processes and , the time-varying wavelet covariance is defined as: As in the case of the variance decomposition (13), if the wavelet covariance does not depend on time, it produces a decomposition of the covariance between and across scales t X t Y j λ : Next, let us define the (time invariant) wavelet correlation coefficient for scale j λ via: The quantity above measures the strength and direction of linear dependence between two processes for a given decomposition level j (scale j λ ). Finally, the wavelet cross-covariance and cross-correlation are given as: As mentioned in section 2.1, MODWT-based estimators have generally better statistical properties compared with their DWT-based counterparts. Firstly, MODWT coefficients produce better estimates of the wavelet variance in terms of efficiency. Secondly, they give an estimator of the wavelet covariance whose variance does not depend on the true time lag between time series. Furthermore, decimation by two affects the lag-resolution of the DWT-based estimators of the wavelet cross-covariance and cross-correlation, so they should not be used in practice (see Percival and Walden 2000: 308-310;Gençay et al. 2002: 252-253). For these reasons, we concentrate on estimation using the MODWT coefficients.
An unbiased estimator of the wavelet variance is defined as: where t j W , are the MODWT wavelet coefficients, is the length of the wavelet filter for scale is the number of wavelet coefficients unaffected by the extrapolation method at the end of the sample.
Estimates of wavelet covariance and wavelet correlation are computed via the following formulas: while an unbiased estimate of the wavelet cross-covariance is obtained via: Methods of constructing CIs for the quantities described in this section are discussed by Percival (1995) and Whitcher (1998) (see also Percival and Walden 2000;Serroukh et al. 2000;Whitcher et al. 2000;Gençay et al. 2002). All the above estimators may be based on only a portion of wavelet coefficients, which provides estimates of the local versions of the wavelet quantities. A good time resolution is exactly what the non-decimated DWT offers and -together with certain simplifications in obtaining global spectral estimates -is the most important characteristic of the approach presented here. The same holds for the MODHWT-based quantities described in §2.4.

MODHWT
The MODHWT makes use of a recently introduced class of filters based on Hilbert wavelet pairs (HWPs) and utilises a non-decimated (maximal overlap) version of the dual-tree complex wavelet transform of Kingsbury (1998Kingsbury ( , 2001 and Selesnick (2001). 7 The approach was advocated by Whitcher and Craigmile (2004) (see also Whitcher et al. 2005). The filters in a HWP are approximate Hilbert transforms of each other and, as in the case of the usual DWT, form a basis for the collection of orthogonal bandpass filters. This time, however, the approximate analyticity of the filters enables us to compute the quantities that directly correspond to the appropriate bivariate Fourier spectra. 8 Let and be conjugate quadrature mirror filters. The father and mother wavelets are obtained via: Now consider another pair of such filters: and that define another couple of father and mother wavelets: and . We say that is the Hilbert transform of if: For an introduction, see Selesnick et al. (2005). where and are the Fourier transforms of and , respectively. This means that the wavelets are π out of phase with each other. The following theorem was proved by Selesnick (2001): 9 If transfer functions of two scaling filters fulfil the condition: , then the corresponding wavelets are a Hilbert transform pair.
The condition (25) says that the digital filter should be a half-sample delayed version of , i.e.
. As a half-sample delay cannot be implemented with finite impulse response filters, only approximate solutions are available.
In following sections, we use mainly the HWP filters introduced by Selesnick (2002). All the filters we apply below have the following property: the filters in the Hilbert pair are of the same length and have the same squared gain functions. In Selesnick's so-called 'common factor approach', firstly an all-pass filter with an approximately constant fractional group delay is constructed and then orthonormal filters are built via a linear system of equations and spectral factorisation. Under a specified degree of approximation (L) to the half-sample delay, the design procedure produces short filters with a given number of VMs (K). The length of each HWP(K, L) filter equals ) ( 2 L K + . In our study, we apply mid-phase solutions for HWP(3, 3), HWP(4, 2), HWP(3, 5) and HWP(4, 4) and denote them 'kKlL'. An alternative approach introduced by Kingsbury (2001) produces the socalled Q-shift (quarter-shift) filters, which are approximately linear phase and whose wavelets in the Hilbert pair are mirror images of each other. For comparison purposes, in our simulation analysis we also use the 6-tap Q-shift filter of Kingsbury (2001) with one VM and the 12-tap Q-shift filter of Tay et al. (2006) with five VMs and denote them 'kin' and 'tkp12', respectively.
and are not much affected either. As such, the HWPs can be seen as short-term versions of the cosine and sine waves forming the classic Fourier transform.
The MODHWT consists in a simultaneous application of a pair of wavelet (and scaling) filters in their non-decimated (maximal overlap) forms. As a result, two sequences of coefficients are obtained, which are the real and imaginary parts of the final wavelet coefficients. In other words, the following filters are used: . These filters produce the complex wavelet and scaling coefficients in the following form: www.economics-ejournal.org As discussed in Selesnick et al. (2005), the simplest approach towards invertible analytic wavelet transform is via DWT post-processing. However, this results in a slightly higher computational complexity as in such a case we operate on two parallel complex wavelet transforms. Furthermore, performing the Hilbert transform as the first one is not recommended since then we lose the possibility of optimising it for all scales simultaneously. The dual-tree complex wavelet transform of Kingsbury and Selesnick is based on two real orthogonal wavelet filters with the Hilbert transform built into them. Thanks to this the Hilbert transform automatically adapts to the wavelet scales. This feature makes the approach particularly attractive compared with other time-frequency methods producing instantaneous amplitudes, phases and frequencies such as the classic demodulation method (see, e.g., Granger and Hatanaka 1964;Priestley 1981) or the modern Hilbert-Huang transform (see Huang and Shen 2005).

Wavelet Analysis of Coherence and Phase Angle
This section includes our theoretical considerations and starts with a detailed description of the relationship between the MODHWT-based wavelet spectra and the appropriate Fourier quantities. The wavelet analogues of Fourier spectral characteristics of bivariate time series were introduced by Whitcher and Craigmile (2004). Let X t j W , and Y t j W , be the complex-valued wavelet coefficients obtained via filtering and . Assuming that the wavelet filters applied have enough VMs to eliminate any deterministic trend components in the series, the time-varying wavelet spectrum of for scale denote the time-varying wavelet cospectrum and quadrature spectrum (quad-spectrum), respectively. If the wavelet cospectrum and quadrature spectrum do not depend on time, it is possible to relate them to the appropriate Fourier quantities. Let the process be covariance stationary with an absolute summable cross-covariance sequence. We will denote its cross- . 10 The wavelet cospectrum for scale j λ is then: where ,l j h . As in our case the two squared gain functions in (30) are identical, we obtain: . Further, assuming that the wavelet filter is long enough to be considered a good approximation to an ideal bandpass filter, we have: As the quantity: _________________________ 10 For wavelet filters with enough VMs, the discussion concerning wavelet co-and quadrature spectra can be directly generalised to the case of non-stationary processes with stationary backward differences. To this end, we consider two integrated processes: whose differences of order and , respectively, are jointly stationary. Then, following Whitcher and where . Note, however, that the generalisation does not apply to the part of our analysis that utilises the complex scaling coefficients www.economics-ejournal.org . If it is possible to assume that the Fourier cross-spectrum is piecewise constant over the octave frequency bands, estimators of may serve to consistently estimate the Fourier cospectrum. j C 11 In any case however, the wavelet quantities discussed here provide piecewise constant approximations to their Fourier counterparts and summarise the information included in the cross-spectrum in a way similar to the wavelet variance in the univariate spectral analysis (see Percival 1995).
In order to obtain similar results for the wavelet quadrature spectrum, we recall the analyticity property of the HWP. The condition (24) is equivalent to: Making use of this and utilising the approximations: (see Percival and Walden 2000: 476), we obtain: _________________________ 11 Note that such an assumption is valid only for rather special kinds of relationships such as the 'fixed angle lag' relationship (Granger and Hatanaka 1964: 98) with a constant amplitude spectrum over the frequency bands of interest or a linear regression without delay.
Unfortunately, it turns out that the first approximation assumed in (34) is of little use for the first several decomposition levels, as the approximately analytic mother wavelets and are most helpful in describing the behaviour of the associated wavelet filters at level j as . In order to obtain more practical results, one is advised to apply different orthogonal quadrature mirror filters in the first stage of Mallat's pyramid algorithm, namely filters approximately satisfying the condition: Selesnick et al. 2005). This condition is different (and easier to implement) than is the half-sample delay requirement: To see that it solves the approximation problem let us consider in more detail the implementation of the MODHWT. We will now discriminate between the level 1 and level j (j = 2, 3, …) filters. We maintain the previous notation for the transfer functions of the remaining filters, while the transfer functions of the first level scaling and wavelet filters will be denoted as: , respectively. For the first level scaling filters, we have: www.economics-ejournal.org Using this and (11a), for the first level wavelet filters, we obtain: 12 The imaginary part of the second stage wavelet filter is given as: If the level 2 scaling filters approximately satisfy the half-sample delay condition, we further obtain the following approximation for : while the appropriate relation for 0 < f is obtained via complex conjugation. Therefore, at the second decomposition level we do not need to substitute for (34) and the approximation is valid. A similar relationship holds for all the subsequent stages. To see this let us consider the j-th decomposition level and the frequencies satisfying . We start by writing: For the first factor in (39) is given as: and the whole expression is then: _________________________ 12 Here we assume that the level one real and imaginary filters have the same (even) length. We change this assumption further in our computations by considering the imaginary filters to be of length L + 1, where L is even and equals the length of the real parts of the filters. This, however, does not change the result that follows.
The first stage complex scaling filter is easily obtained via the translation of any real scaling filter by one sample and using it as the imaginary part of the resulting filter. Then, the complex wavelet filter is computed via the quadrature mirror relationship applied separately to these two parts. Because the transfer functions at the first stage of the wavelet decomposition obviously do not satisfy the analyticity property, all the quantities computed with the help of the wavelet quad-spectrum should basically be interpreted starting from the second level. However, this does not cause a problem for business cycle studies, which are typically based on monthly or quarterly data.
As in the case of the cospectrum, the wavelet quad-spectrum enables us to compute the average value of its Fourier analogue in the interval ] , [ Finally, we arrive at the following approximation to the Fourier cross-spectrum: www.economics-ejournal.org where is the average value of the Fourier cospectrum in the interval . Then, we obtain: where in the last equality it is assumed that the value of the Fourier quad-spectrum in the interval ] , 0 [ 1 2 1 + J is constant and equal to .
1 + J Q Next, as in Whitcher and Craigmile (2004), we consider the time-varying wavelet cross-amplitude spectrum: the time-varying wavelet phase spectrum (wavelet phase angle -WPA): and the time-varying wavelet magnitude squared coherence: 13 denote the time-varying wavelet spectra equal to twice the wavelet variance, . The Schwartz inequality for complex random variables guaranties that

_________________________
13 We refer to this as the wavelet coherence or wavelet coherence spectrum. Its square root is called the wavelet coherency.
www.economics-ejournal.org Let us consider a simple example of a stationary bivariate process -a linear regression with delay in the form: where and t X t η are stationary processes, uncorrelated with each other at all leads and lags. Then, the Fourier quantities take the following values: Furthermore, we can use also the time delay: as well as the . Assuming that the wavelet transform produces a bandpass white noise, the wavelet co-, quad-and amplitude spectra are the following: where the last approximation takes place for high enough decomposition levels j. 14 The wavelet coherence is given as: and for higher decomposition levels this can be approximated via: Assuming that the individual Fourier spectra are approximately piecewise constant over octave frequency bands, for high enough scales the wavelet coherence should provide an acceptable approximation to the appropriate Fourier quantity. The wavelet phase spectrum is: In order to obtain a wavelet estimator of the parameter τ, we also introduce a quantity which we call the wavelet time delay. This is defined as: where is the centre frequency of the octave band, computed as the arithmetic mean of the upper and lower cut-off frequencies, i.e. π ≤ x , the last approximation will work for (j = 2, 3, …) -comp. Percival and Walden (2000: 344).
and in the example above for we get: A slightly more general situation arises in the case of the so-called time delay estimation problem described as follows. We consider two spatially separated sensor measurements, and , of an unobserved signal . Let us assume that the processes satisfy: where , t S Xt η and Yt η are stationary and mutually uncorrelated at all leads and lags. Then: From this we obtain that the phase spectrum is as in the previous example, i.e. The estimators of (46)-(48) and (56) are obtained by replacing the wavelet cospectrum and quad-spectrum as well as the wavelet individual spectra with their estimates computed via smoothing in time. This smoothing is particularly necessary for the estimation of the wavelet coherence. In Whitcher and Craigmile (2004), a simple two-sided moving average is suggested, and this is the approach taken here as well. Figures 2 and 3 present the mean estimates of the wavelet coherence, phase spectrum and time delay for samples ranging from 10 to 500 wavelet coefficients not affected by the boundary, obtained with the k4l2 HWP filter for the linear regression with delay (49) with 1 = τ . Figure 2 illustrates the case, when the first stage filters are different and fulfil the one-sample delay condition. The la12 www.economics-ejournal.org Daubechies filter was applied in the real part of the first stage complex filter. For comparison purposes, the appropriate results obtained without this modification are also presented (see Figure 3). As we can see, the modified procedure gives acceptable results starting from scale two, while the simplified method introduces more bias in both the coherence and time delay estimations, especially at the second decomposition level. Another observation concerns the small sample bias of the wavelet coherence estimator, which clearly increases with the scale.  Figure 2. Estimates of the wavelet coherence, wavelet phase spectrum and wavelet time delay with the modified method; the first stage filters are la12 and its one-sample shifted version and the complex filter for the higher levels is k4l2; figures present the theoretical Fourier quantities (thin dotted lines) together with the mean estimates of the corresponding wavelet quantities obtained with 500 replications for samples consisting of 10, 20, …, 500 wavelet coefficients unaffected by the boundary for the linear regression model with delay τ = 1, α = 1 and t and t X η being two independent AR(1) processes with autoregressive parameters 0.8 and unit error variances (thick solid lines). To construct CIs for the wavelet coherence the multivariate process in the form: is considered together with the function: where , is the (1-α/2)-quantile of the standard normal distribution. where 2 α ς To construct confidence intervals for the WPA we assume that 0 ) ( ≠ j XY C λ and take: Then, applying the delta method, we arrive at: where the large sample variance, , is computed as previously by utilising an appropriate vector of partial derivatives equal to: Then, the confidence intervals are given in the form similar to (62). Multiplying them by a constant will produce approximate CIs for the wavelet time delay. Finally, assuming that the wavelet gain (57) is positive, similar reasoning to the presented above will provide approximate CIs for this quantity.
Estimates of the large sample variance of the wavelet spectra estimators can be obtained via nonparametric kernel methods. We examined the properties of two kernel estimators: one based on the popular Bartlett kernel and the other based on the truncated (rectangular) kernel, which however can provide negative values of the estimates. The speed of convergence of the kernel estimators to the experimentally driven true values depends on the kind of the wavelet spectrum, the type of the bivariate data-generating process and the scale of the analysis. We concentrated on a linear regression with delay, for which different stationary AR(1) models for and t X t η were considered. For lower scales -λ = 2 (4) -good results in terms of unbiasedness were obtained, even in samples of length as small as N = 50 (75) wavelet coefficients, while for λ = 8 acceptable approximations started with N = 200, and even more data were necessary for λ = 16. The experimentally determined best values of the truncation parameters M are both sample size-and scale-dependent. 15

Wavelet Estimation of the Time Delay -a Simulation Study
In this section, we summarise the results of our simulation experiments examining the statistical properties of two wavelet estimators of the time delay parameter. We concentrate on model (49) and analyse mainly the small sample performance of the estimators in order to recommend a method for examining short-and mediumterm lead-lag relations for octave frequency bands. In particular, such a method might be of interest in business cycle studies, as it should be useful when analysing changing patterns of business and growth cycle synchronisation. The estimators compared are: the wavelet cross-correlator (WCC) that is based on maximising the values of the cross-covariance estimates, i.e.: and the estimator of the wavelet time delay (56), i.e.:

_________________________
15 For exactly the same DGP and data lengths as in the case of the estimates in Figures 2 and 3, acceptable results in larger samples were obtained with for the rectangular kernel and for the Bartlett kernel, where j = 2, 3, … denotes the decomposition level. More details of this examination are available upon request. Owing to a programming error, the conclusions from this part of our simulations differ slightly compared with those reported in a previous version of the paper.
which we call further the WPA-based delay estimator. It is worth stressing that we do not maximise the absolute value of the cross-covariance in (65) since a large negative covariance would be interpreted as an anti-phase relationship. The Cramér-Rao lower bound on the variance of any unbiased time delay estimator was derived as (see Carter 1987): where is the Fourier coherence of the processes under study (in our case, the wavelet coefficients (67) predicts that the variance of an optimal estimator decreases with the value of coherence, the signal bandwidth and the centre frequency. For the WCC, it is known that for jointly stationary processes and large enough data samples the Cramér-Rao lower bound is automatically achieved in the case of the signal and noise processes with spectra that are flat over the same range of frequencies and zero outside this range (see Scarbrough et al. 1981, Carter 1987. If the spectra of and t X t η are relatively featureless within octave frequency bands, the MODWT wavelet coefficients of and will be approximately bandpass white noises and the WCC becomes efficient asymptotically. However, the actual performance of the asymptotically efficient estimator can be much worse, especially for low SNRs (see the simulation results in Scarbrough et al. 1981, Carter 1987. Regarding the estimators based on the Fourier phase angle, it has been proven that they are fully consistent with other asymptotically optimal methods after regression analysis is applied to the phase data (see Piersol 1981).
When discussing the properties of the estimators (65) and (66), it is worth underlining that the WPA-based method enables us to estimate delays that are not integer multiples of the sampling period In addition, it guarantees a maximal time localisation limited only by the length of the applied filter. By contrast, the WCC makes it possible to estimate delays that are longer than half of the centre period, 0 , it is asymptotically unbiased, even for the first stage of analysis, and it can be based on shorter filters, as the approximately analytic complex wavelet filters / 1 j f are usually longer than are the real filters with similar squared gain functions. The length of the wavelet filter plays a crucial role in empirical examination, since it directly influences the number of wavelet coefficients that are unaffected by the extrapolation method at the ends of the sample and, therefore, determines the maximal number of decomposition levels as well as the precision of estimation. In our simulations, the following data generating process was used: , and the following wavelet filters: la12 for the WCC and la12 (first stage) + k4l2 (higher stages) for the WPA-based estimator. 16 The search range for the WCC was τ + ±10 . In each case, 1000 replications were run. 17 In the presentation, we also include the outcomes for the first decomposition level, largely because to some extent they are comparable to the other stages. The findings resulting from the experiments are summarised below: 18 • For wide ranges of SNRs and scales, the WPA-based estimator is better than the WCC in small samples (see Figure 5). For the majority of outcomes, the relative efficiency of the two methods defined as ) RMSE( ) RMSE( WCC WPA _________________________ 16 These wavelet filters were chosen because of their popularity and also to guarantee maximal similarity in the implementation of the WCC and WPA methods: la12 and k4l2 are of the same length (L = 12) and have similar squared gain functions. However, basically the same results were obtained for, e.g., la8 + k3l3, la12 + k4l4 and real part of k4l2 + k4l2. Selesnick's HWP(K, L) filters outperformed the two Q-shift filters that were also considered, i.e. kin and tkp12. 17 All computations, including the empirical part, were executed in Matlab. Numerical codes are available by e-mailing the author. 18 We also comment shortly on other experiments we performed. More detailed results are available upon request.
www.economics-ejournal.org increases with the scale, the sample size and partially also with the SNR. In large samples (see Figure 6), the WCC dominates the WPA-based estimator or their performance is similar. In small samples, the relative efficiency of these methods depends largely on the search range for the WCC, although for similar ranges of delays for both methods the WPAbased estimator gives better results.
• In larger samples, the root mean square errors for both estimators increase with the scale. It is generally advisable to assume smoothing windows with a length proportional to the scale (comp. Cohen and Walden 2010). In addition, we can see that even in large samples in the case of high decomposition levels and low SNRs both estimates can be very imprecise.   • For low SNR, both estimators show small sample bias, although with opposite signs: the WPA-based estimator -towards 0, while the cross-correlator -in the opposite direction (see Figure 4). This suggests that when an estimate of the time delay parameter is needed, the usual biased estimator of the crosscovariance might be preferred in small samples. Our other experiments also demonstrate that the small sample bias of the WCC largely depends on whether the search range for the WCC is symmetric around the true value of the delay. • In a different set of simulation experiments, in which the WPA-based estimator (66) is replaced with its integer-valued counterpart: the outcomes in small and moderate samples are even more in favour of the WPA-based estimator as (66*) removes the bias caused by the centre frequency imprecision.
• Other experiments not presented here indicate that if the bivariate data generating process is jointly stationary, including observations affected by the periodic extension at the end of the sample can give better small sample results with both methods at higher scales, where the number of affected coefficients is large. This, however, does not take place for non-stationary processes with a stationary delay, which exhibit a slowly varying local mean. In such a case, a periodic extension with level adjustment should provide better results. • All the above observations are unchanged across different wavelets, although the outcomes obtained with the WPA-based method depend on the analytic properties of the HWP filters. Good analytic wavelets, however, produce almost identical results (for example, k4l2, k3l3). Different values of β and γ (including the non-stationary case) do not change the conclusions either.
The main finding of this section is that in business cycle studies, which are typically based on relatively short time series, the WPA-based methodology seems to be particularly attractive and can be used at least as a supplementary method. It is worth stressing that, in addition to its good localisation properties, the WPAbased estimator is also simple and efficient computationally. For these reasons, we believe it can be recommended for empirical analysis on business cycle synchronisation.

Empirical Examination
In the empirical study we use quarterly GDP volume estimates from the OECD Quarterly National Accounts (measure: VOBARSA) covering the period from the first quarter of 1960 till the second quarter of 2010 (202 observations) for the following 11 countries: Austria, Belgium, Finland, France, Germany, Greece, Ireland, Italy, the Netherlands, Portugal and Spain. In addition, the OECD GDP volume for the euro area (16 countries) is used, which covers the shorter period from 1995 till the end of the sample (62 observations). The examination is divided into two parts. In the first part, the local wavelet variance analysis is performed, while the second deals with the local and global wavelet analysis of synchronisation.

Business Cycle Variability
Our examination of business cycle variability was performed with the help of the d4 Daubechies wavelet filter of length 4, which guarantees very good localisation properties. The results are presented in Figures 7 and 8. 1960 1972 1985 1997  Firstly, we notice very similar patters of volatility changes across countries in our sample, except for Finland (scales 4 and 8) and Ireland (all scales). Furthermore, the contributions of different scales to the total variance as well as the estimates of the wavelet variance alone vary little across the economies. For some countries, we observe a systematic decline in the variance at all decomposition levels, which started at the beginning of our sample, as is the case for Germany (except for the more volatile period around the reunification as well as for the highest scale) and Spain. It is seen that the oil price shocks of 1973 and 1979 were captured almost entirely by the shortest components of business cycle fluctuations. Thanks to this, the scale 8 wavelet variance provides a clearer view of the Great Moderation, revealing that the process might have started well before the mid-1980s (comp. Aguiar-Conraria and Soares 2010 for similar evidence for the United States obtained with the continuous wavelet methodology). Moreover, the most recent perturbations (the financial crisis of [2007][2008][2009] are becoming apparent at the lowest decomposition level as can be seen in the case of the euro area GDP.

Business Cycle Synchronisation
Business cycle synchronisation in the euro area was examined with the help of local wavelet correlations, global and local wavelet coherences and global and local wavelet time delays. Figure 9 presents the running wavelet correlations for scales 4, 8 and 16 computed after the MODWT based on d4 wavelet had been applied to the observations. We decided to treat Germany as the reference country because of the high correlation between its scale 4 wavelet coefficients and those for the euro zone compared with the appropriate correlation for France and the euro zone (see Figure 9). The most interesting finding is that for the majority of countries in the sample we observe a systematic increase in the strength of the instantaneous relationships between the business cycles of the examined countries starting from the second half of the 1980s, especially for scales 4 and 8 corresponding to oscillations with period lengths below 4 and 8 years, respectively. The change in the patterns of the dynamic correlations agrees with the introduction of the Single European Act, which was signed in 1986 and came into effect in 1987. It may also be observed that for the highest scale considered (cycles of 8 years and above) there often is the opposite tendency, i.e. a cyclical divergence. For Greece, this seems to take place also at scale 8.  Figure 11. Running wavelet coherence and wavelet time delay for scales 4 and 8 corresponding to oscillations with period lengths 8-16 (2-4 years) and 16-32 (4-8 years); the first stage filters are la12 and its one-sample shifted version and the complex filter for the higher levels is k4l2; data windows consist of 30 non-boundary wavelet coefficients for scale 4, and 40 -for scale 8, circularly shifted to align them to the real data; step = 1. the majority of countries in the sample at the third decomposition level (for shorter cycles). For Austria and the Netherlands, high correlations are also present at the fourth level. These two countries' cycles show the strongest relationship with the German business cycle. Positive delays mean that the German cycle is behind the other countries' cycles, as is the case for the long French cycle and the short Irish cycle. Instantaneous dependencies with the German business cycle take place for countries such as Belgium, the Netherlands, Greece and Italy as well as for the lower decomposition levels for almost all countries in the sample. The local analysis in Figure 11, performed in windows of a scale-dependent size, reveals that the shorter cycles are becoming more synchronised from the middle of our sample, while in the case of the longer cycles the lead-lag relations are quite stable over time (partially because only a couple of coefficient windows were available), although the wavelet coherences are rising.
The overall conclusion from the real and complex wavelet analysis is that the synchronisation between euro zone business cycles started to rise after the first important steps toward European integration were taken. This is in line with the endogeneity hypothesis of the optimum currency area criteria as stated by Frankel and Rose (1998). Finally, Figure 12 presents a comparison of local wavelet time delay estimates obtained using the WCC and WPA methods. It turns out that these methods produce largely similar results.

Conclusions
The non-decimated real and complex discrete wavelet transforms provide a summary of the evolutionary spectral and cross-spectral properties of processes under scrutiny with high computational efficiency, good localisation properties and without an excessive redundancy of information that occurs when using the continuous wavelet methodology. These features together with a fresh look at an old problem seem to be the main reasons why the approach might be worth considering in business cycle examination. The paper discussed some of the questions arising in the discrete wavelet analysis of popular bivariate spectral quantities such as the amplitude, coherence and phase spectra and the frequency-dependent time delay. In particular, we showed how the wavelet bivariate spectra can serve to approximate the corresponding Fourier quantities and discussed certain implementation issues. Our simulation study of the properties of two wavelet estimators of the time delay parameter pointed out the practical relevance of the wavelet phase angle-based estimator suggested here, which can be used at least as a supplementary method of examining short-and medium-term lead-lag relations for octave frequency bands.  Figure 12. Running wavelet estimates of time delay for scales 4 and 8 corresponding to oscillations with period lengths 8-16 (2-4 years) and 16-32 (4-8 years); the solid blue line is the result obtained with the WPA-based estimator and the dashed red line -with the WCC; in the WPA-based method the first stage filters are la12 and its one-sample shifted version and the higher level filter is k4l2; the WCC is based on la12; the search ranges for the WCC are [-6, 6] and [-12, 12] for scales 4 and 8, respectively; data windows consist of 30 wavelet coefficients unaffected by circularity for scale 4, and 40 -for scale 8, the coefficients were circularly shifted to align them to the real data; step = 1; the numbers on the horizontal axis are the mid-points of the subsamples.
The complex discrete wavelet methodology was illustrated with an examination of business cycle synchronisation in the euro zone. The study was also supplemented with a wavelet analysis of the variance and covariance of European business cycles. The empirical examination gives some new arguments in favour of the endogeneity hypothesis of the optimum currency area criteria as well as the early start of the Great Moderation in Europe.