2.1 The Algorithm

Suppose a function f( · ) is observed at a set of n irregularly spaced locations, $\underset{\_}{x}=(x_1,\cdots ,x_n)$. The proposed lifting scheme aims to decompose the data collected over the irregularly sampled grid, {(x_i, f_i = f(x_i))}ⁿ_{i = 1}, into a set of R smooth coefficients and (n − R) complex-valued detail coefficients, with R the desired resolution level. The quantity R is akin to the primary resolution level in classical wavelet transforms, see Hall and Patil (1996) for more details.

We propose to construct a new complex-valued transform that builds upon the LOCAAT paradigm of Jansen, Nason, and Silverman (2001, 2009), shown to efficiently represent local signal features in the fields of nonparametric regression (Nunes, Knight, and Nason 2006; Knight and Nason 2009) and spectral estimation (Knight, Nunes, and Nason 2012). We shall therefore refer to our proposed algorithm under the acronym $\mathbb{C}$-LOCAAT.

Similar to the real-valued LOCAAT algorithm, $\mathbb{C}$-LOCAAT can be described by recursively applying three steps: split, predict, and update, which we detail below. At the first stage (n) of the algorithm, the smooth coefficients are set as c_{n, k} = f_k, the set of indices of smooth coefficients is S_n = {1, …, n} and the set of indices of detail coefficients is D_n = ∅. The (irregular) sampling is described using the distance between neighboring observations, and at stage n we define the span of x_k as $s_{n,k}=\frac{x_{k+1}-x_{k-1}}{2}$. The sampling irregularity is intrinsically linked to the notion of wavelet scale, which in this context becomes continuous, as opposed to dyadic in the classical wavelet settings; this results in each coefficient having an associated scale across a continuum. This aspect will be discussed in detail following the introduction of the $\mathbb{C}$-LOCAAT algorithm.

In the split step, a point j_n to be lifted is chosen. Typically, points from the densest sampled regions are removed first, but other predefined removal choices are also possible (see Section 2.3). We shall often refer to the removal order as a trajectory.

In the predict step, the set of neighbors (J_n) of the point j_n is identified and used to estimate the value of the function at the selected point j_n. In contrast to real-valued LOCAAT algorithms, this is achieved using two prediction schemes, each defined by its respective filters, L and M. The filter L corresponds to estimation via regression over the neighborhood, as is usual in LOCAAT. To extract further information from the signal, our proposal is to construct the second filter (M) orthogonal on L, to ensure that it exploits further local signal information to the filter L. Section 2.2 discusses this in detail.

The prediction residuals from using the two filters are given by (1) $$\begin{array}{c}\hfill {\lambda }_{j_n}=l_{j_n}^n{c}_{n,j_n}-\sum _{i\in J_n}{l}_i^n{c}_{n,i},\end{array}$$(1) (2) $$\begin{array}{c}\hfill {\mu }_{j_n}=m_{j_n}^n{c}_{n,j_n}-\sum _{i\in J_n}{m}_i^n{c}_{n,i},\end{array}$$(2) where ${\left\{l_i^n\right\}}_{i\in J_n\cup \left\{j_n\right\}}$ and ${\left\{m_i^n\right\}}_{i\in J_n\cup \left\{j_n\right\}}$ are the prediction weights associated with L and M.

The complex-valued detail (wavelet) coefficient we propose is obtained by combining the two prediction residuals (3) $$\begin{array}{c}\hfill d_{j_n}={\lambda }_{j_n}+\mathrm{i}{\mu }_{j_n}.\end{array}$$(3) In the update step, the smooth coefficients ${\left\{c_{n,i}\right\}}_{i\in J_n}$ and spans of the neighbors ${\left\{s_{n,i}\right\}}_{i\in J_n}$ are updated according to filter L: (4) $$\begin{array}{ccc}\hfill c_{n-1,i}& =& c_{n,i}+b_i^n{\lambda }_{j_n},\hfill \\ \hfill s_{n-1,i}& =& s_{n,i}+l_i^n{s}_{n,j_n}\forall i\in J_n,\hfill \end{array}$$(4) where bⁿ_i are update weights. In practice, the update weights are chosen such that the mean of the series is preserved throughout the transform, thus preserving the characteristics of the original series (Jansen, Nason, and Silverman 2009). One such choice is to set $b_i^n=\left(s_{n,j_n}{s}_{n-1,i}\right)/\left({\sum }_{i\in J_n}{s}_{n-1,i}^2\right)$. The neighbors’ spans update accounts for the modification to the sampling grid induced by removing one of the observations. Updating according to the L filter only ensures that there is a unique coarsening of the signal for both the real and imaginary parts of the transform.

The observation j_n is then removed from the set of smooth coefficients, hence after the first algorithm iteration, the index set of smooth and detail coefficients are S_{n − 1} = S_n\{j_n} and D_{n − 1} = {j_n}, respectively. The algorithm is then iterated until the desired primary resolution level R has been achieved. In practice, the choice of the primary level R in LOCAAT lifting schemes is not crucial provided it is sufficiently low (Jansen, Nason, and Silverman 2009), with R = 2 recommended by Nunes, Knight, and Nason (2006).

After observations j_n, j_{n − 1}, …, j_{R + 1} have been removed, the function can be represented as a set of R smooth coefficients, ${\left\{c_{r-1,i}\right\}}_{i\in S_R}$, and (n − R) detail coefficients, ${\left\{d_k\right\}}_{k\in D_R}$ (D_R = {j_n, …, j_{R + 1}}). As in classical wavelet decompositions, the detail coefficients represent the high-frequency components of f( · ), while the smooth coefficients capture the low-frequency content in the data.

The lifting scheme can be easily inverted by recursively “undoing” the update, predict, and split steps described above for the first filter. Specifically, the update step is first inverted: $c_{n,i}=c_{n-1,i}-b_i^n{\lambda }_{j_n},\forall i\in J_n$, then the predict step is inverted by (5) $$\begin{array}{c}\hfill c_{n,j_n}=\frac{{\lambda }_{j_n}-{\sum }_{i\in J_n}{l}_i^n{c}_{n,i}}{l_{j_n}^n}\text{or}\end{array}$$(5) (6) $$\begin{array}{c}\hfill c_{n,j_n}=\frac{{\mu }_{j_n}-{\sum }_{i\in J_n}{m}_i^n{c}_{n,i}}{m_{j_n}^n}.\end{array}$$(6)

Undoing either predict (5) or (6) step is sufficient for inversion. As for real-valued lifting, inversion can also be performed via matrix calculations due to the transform linearity. However, using (5) for inversion is generally computationally faster, especially for large n.

Wavelet lifting scales. The notion of wavelet scale in this context becomes continuous and is intrinsically linked to the data sampling structure and trajectory (removal order) choice. Denote the lifting analog of the classical wavelet scale for a detail coefficient $d_{j_k}$ by ${\alpha }_{j_k}={log}_2\left(s_{k,j_k}\right)$, with low α-values corresponding to fine scales. To give lifting scales a similar interpretation to the classical notion of dyadic wavelet scale, we group wavelet functions of similar α-scales into discrete artificial levels {ℓ_i}^J*_{i = 1}, as proposed by Jansen, Nason, and Silverman (2009), for a chosen J*. The further use of artificial scales is discussed in Sections 2.3 and 3 (under the spectral estimation context) and in Appendix B (under the nonparametric regression context). Note that the usage of the same lifting trajectory for the two lifting branches (coupled with the one filter update) ensures that our proposed lifting transform generates a common scale for both real and imaginary parts. In other words, at each stage of the algorithm there is just one set of smooth coefficients associated with a unique set of scales.

2.2 Filter Construction

The proposed complex-valued lifting transform is illustrated schematically in Figure 1 in terms of two general prediction filters L and M. As already explained, we construct the second filter (M) orthogonal on L, thus ensure different signal content extraction.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 1. The complex-valued lifting scheme ($\mathbb{C}$-LOCAAT). Solid lines correspond to the steps of the standard LOCAAT lifting scheme whereas dotted lines indicate the extra prediction step required for the complex-valued scheme. After (n − R) applications, the function can be represented as a set of R smooth coefficients ${\left\{c_{r-1,i}\right\}}_{i\in S_{n-R}}$ and (n − R) detail coefficients ${\{{\lambda }_{j_k}+i{\mu }_{j_k}\}}_{k\in D_{n-R}}$, each associated with a particular scale ${\left\{{\alpha }_{j_k}\right\}}_{k\in D_{n-R}}$.

Display full size

For clarity of exposition, let us consider a LOCAAT scheme with a prediction step based upon two neighbors in a symmetrical configuration. The regression over the neighborhood generates prediction weights for the two neighbors, let us denote them by l₁ and l₃ (see Equation (1)); this prediction step can also be viewed as using a three-tap prediction filter (L) of the form (l₁, 1, l₃), which depends on the sampling of the observations $\underset{\_}{x}=(x_1,\cdots ,x_n)$ (Nunes, Knight, and Nason 2006). We determine the unique (up to proportionality) three-tap filter M that is orthogonal on L and ensures at least one vanishing moment. Hence, we can express the set of filter pairs as having the form $$\begin{array}{ccc}\hfill \mathbf{L}& =& (l_1,1,l_3),l_1,l_3>0\hfill \\ \hfill \mathbf{M}& =& (m_1,m_2,m_3),\hfill \end{array}$$and l₁m₁ + m₂ + l₃m₃ = 0 (i.e., L · M = 0) and l₁ + l₃ = 1, m₁ + m₃ = m₂ (i.e., ensure one vanishing moment). The solution to these constraints can be parameterized as $\mathbf{M}=(-\frac{1+l_3}{1+l_1}m,\frac{l_1-l_3}{1+l_1}m,m)$. The proportionality constant can be determined by bringing both filters L and M to the same scale through ‖L‖ = ‖M‖, which yields $m=\frac{l_1+1}{\sqrt{3}}$. Hence, the solution can be succinctly written as M = (Am, (1 + A)m, m) with $A=\frac{l_1-2}{l_1+1}$ and $m=\frac{l_1+1}{\sqrt{3}}$. This particular example of the lead filter L represents a prediction scheme using linear regression with two neighbors in a symmetrical configuration. This is a choice that has proved to be successful both for (real-valued) nonparametric regression (Nunes, Knight, and Nason 2006; Knight and Nason 2009) and for (real-valued) spectral estimation (Knight, Nunes, and Nason 2012).

Since L can be viewed as a prediction filter for a real-valued LOCAAT scheme, we can also employ the adaptive prediction filter choice of Nunes, Knight, and Nason (2006) in our proposed construction. The “best” local regression (order and neighborhood) is chosen at each predict step, subject to yielding minimizing the detail coefficients. Consequently, we obtain an adaptive complex-valued lifting transform, with the highly desirable flexibility of being able to adapt to the local characteristics of the data–see Appendix B in the supplementary material for an illustration of this adaptiveness in the nonparametric regression setting.

The orthogonality of the two filters M and L also mirrors the attractive properties of Fourier sinusoids, hence this choice results in an interpretable quantification of phase, which shall further be exploited according to the context—by phase alteration when denoising real-valued signals, or by ensuring phase preservation in the context of spectral estimation.

A further insight and justification of the proposed filter choice is provided in Appendix C in the context of coherence and phase estimation.

2.3 The Nondecimated Complex-Valued Lifting Transform

In the classical wavelet literature, the nondecimated wavelet transform (NDWT) (Nason and Silverman 1995) has properties that make it a better choice than the discrete wavelet transform (DWT) for certain classes of problems, see, for example, Percival and Walden (2006). The concept is akin to basis averaging, and has delivered successful results in both nonparametric regression and spectral estimation problems, not only in the classical wavelet setting (NDWT) but also for irregularly spaced data through the nondecimated lifting transform (NLT) (Knight and Nason 2009; Knight, Nunes, and Nason 2012).

In this section, we also exploit the benefits of this nondecimation paradigm for irregularly sampled data and to this end, we shall introduce the complex-valued nondecimated lifting transform ($\mathbb{C}$NLT). However, note that our use of the term “nondecimation” differs from the classical NDWT. Specifically, due to the irregular sampling structure, nondecimation cannot be performed via decomposing shifts of input data without data interpolation.

Although similar in spirit to the NLT, our transform hinges on the proposed complex-valued lifting scheme (Section 2.1) and therefore yields an overcomplete complex-valued data representation, extracting additional signal information. In particular, the $\mathbb{C}$NLT algorithm results in a wavelet transform that yields (complex-valued) wavelet coefficients at each grid point (x) and at multiple scales (α).

Next we shall describe our proposed univariate and bivariate $\mathbb{C}$NLT techniques. We shall show that in the nonparametric regression setting, our univariate proposal significantly outperforms current wavelet and nonwavelet denoising techniques (see Section 4 and Appendix B), while its bivariate extension allows for estimation of the dependence between pairs of series (see Section 3).

2.3.2 Bivariate $\mathbb{C}$NLT

We now consider the extension of $\mathbb{C}$NLT to the analysis of bivariate series.

Same irregular grid. Let us first assume we have observations {(x_i, f¹_i, f²_i)}ⁿ_{i = 1} on two functions f¹ and f², measured on the same x-grid. Apply the univariate $\mathbb{C}$NLT (Section 2.3.1) to each signal, using the same set of trajectories {T_p}^P_{p = 1} for both series.

The identical sampling grids result in an exact correspondence between the coefficients of each series, that is, for each coefficient of the first series there is a coefficient of the second series at exactly the same location and scale (see Figure 2^{Figure 3}(a)). In other words, after application of the $\mathbb{C}$NLT to both series, for each time point, x_k, we obtain two sets of complex-valued detail coefficients ${\left\{d_{x_k}^{1,p}\right\}}_{p=1}^P$ and ${\left\{d_{x_k}^{2,p}\right\}}_{p=1}^P$.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 2. Construction of bivariate $\mathbb{C}$NLT transform for time series observed on the same sampling grid (x refers to time here): (a) univariate $\mathbb{C}$NLT is applied using the same set of trajectories for both series and yields two sets of detail coefficients ${\left\{d_{x_k}^{1,p}\right\}}_{p,k}$ and ${\left\{d_{x_k}^{2,p}\right\}}_{p,k}$; (b) the $\mathbb{C}$NLT transform consists of combinations of coefficients from each series; (c) the detail coefficients are averaged within each scale.

Display full size

Figure 2. Construction of bivariate $\mathbb{C}$NLT transform for time series observed on the same sampling grid (x refers to time here): (a) univariate $\mathbb{C}$NLT is applied using the same set of trajectories for both series and yields two sets of detail coefficients ${\left\{d_{x_k}^{1,p}\right\}}_{p,k}$ and ${\left\{d_{x_k}^{2,p}\right\}}_{p,k}$; (b) the $\mathbb{C}$NLT transform consists of combinations of coefficients from each series; (c) the detail coefficients are averaged within each scale.

Different irregular grids. Let us now assume we have the data {(x¹_i, x²_i, f¹_i, f²_i)}ⁿ_{i = 1} on two functions f¹ and f², measured on the different x-grids.

As the scale associated with each detail coefficient is determined by the trajectory choice, we partition the x-grid into a set of artificial x-intervals {x^(j)}^T*_{j = 1}, where T* is chosen to provide the desired resolution level on the x-axis. As illustrated in Figure 3, the result can be visualized in terms of forming a grid over the area of the resulting detail coefficients.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 3. Construction of bivariate $\mathbb{C}$NLT transform for time series observed on different sampling grids (x refers to time here): (a) each series is lifted individually as described in Section 2.3; (b) the sets of coefficients in each grid square. The coefficients are sampled so that there is the same number in the grid square of each series; (c) the coefficients of each series are combined to form the appropriate bivariate quantities, producing one coefficient to represent each grid square.

Display full size

Figure 3. Construction of bivariate $\mathbb{C}$NLT transform for time series observed on different sampling grids (x refers to time here): (a) each series is lifted individually as described in Section 2.3; (b) the sets of coefficients in each grid square. The coefficients are sampled so that there is the same number in the grid square of each series; (c) the coefficients of each series are combined to form the appropriate bivariate quantities, producing one coefficient to represent each grid square.

Formally, for each artificial x-interval {x^(j)}^T*_{j = 1} and artificial scale {ℓⁱ}^J*_{i = 1}, the set of detail coefficients for each grid square (using trajectories {T_p}^P_{p = 1}) is given by (7) $$\begin{array}{c}\hfill D_{x^{\left(j\right)}}^1\left({\ell }^i\right)=G_m\left(d_{x_k^1}^{1,p}|{\alpha }_{x_k^1}^{1,p}\in {\ell }^i,x_k^1\in x^{\left(j\right)}\right)\end{array}$$(7) (8) $$\begin{array}{c}\hfill D_{x^{\left(j\right)}}^2\left({\ell }^i\right)=G_m\left(d_{x_k^2}^{2,p}|{\alpha }_{x_k^2}^{2,p}\in {\ell }^i,x_k^2\in x^{\left(j\right)}\right),\end{array}$$(8) where $d_{x_k^1}^{1,p}={\lambda }_{x_k^1}^{1,p}+\mathrm{i}{\mu }_{x_k^1}^{1,p}$ and $d_{x_k^2}^{2,p}={\lambda }_{x_k^2}^{2,p}+\mathrm{i}{\mu }_{x_k^2}^{2,p}$ are the complex-valued wavelet coefficients from f¹ and f², and G_m is a random sampling procedure selecting $m_{i,j}=\text{min}(\#\left(d^1\right),\#\left(d^2\right))$ coefficients. Recall that ${\alpha }_{x_k^1}^{1,p}$ and ${\alpha }_{x_k^2}^{2,p}$ represent the scales (log₂ of span) associated with the coefficients $d_{x_k^1}^{1,p}$ and $d_{x_k^2}^{2,p}$. Thus, for each artificial x-interval and scale, we obtain the same number of detail coefficients (although the exact coordinates of the coefficients may differ).

We term these constructions as the bivariate complex nondecimated lifting transform (bivariate $\mathbb{C}$NLT) on the same/different grid(s), as appropriate. Section 3 will discuss applications where the proposed bivariate $\mathbb{C}$NLT construction provides a framework for estimation of the dependence between pairs of series.

3.1 The Complex Lifting Periodogram

Recall that the $\mathbb{C}$NLT provides a set of detail coefficients and associated scales ${\{d_{x_k}^p,{\alpha }_{x_k}^p\}}_{p=1}^P$, where the scale associated with each detail coefficient ${\alpha }_{x_k}^p$ is a continuous quantity. In a spirit similar to that of Knight, Nunes, and Nason (2012), this information will allow a time-scale decomposition (typically termed the (wavelet) periodogram) of the variability in the data, with the crucial difference that the wavelets coefficients are now complex-valued and therefore contain more information. In constructing the periodogram, we use a set of discrete artificial scales, {ℓⁱ}^J*_{i = 1}, which partitions the range of the continuous lifting scales $\left\{{\alpha }_{x_k}^p\right\}$ for all p and k, with J* chosen to provide a desired periodogram “granularity.” Each scale ${\alpha }_{x_k}^p$ will fall into one unique level ℓⁱ for each p and observation x_k; let $P_{i,k}=\{p:{\alpha }_{x_k}^p\in {\ell }^i\}$ denote the set of trajectories such that x_k is associated with a scale in the set ℓⁱ, and n_{i, k} = |P_{i, k}| denote the size of the set. For each time point x_k, k = 1, …, n and artificial scale ℓⁱ, i = 1, …, J*, we introduce the complex lifting periodogram (also referred to in text as $\mathbb{C}$NLT periodogram) $$I_{x_k}\left({\ell }^i\right)=\frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\left|d_{x_k}^p\right|}^2=\frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\left({\lambda }_{x_k}^p\right)}^2+\frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\left({\mu }_{x_k}^p\right)}^2,$$where | · | denotes the complex modulus.

3.2 The Complex Lifting Cross-Periodogram

Similar to other complex wavelet transforms (Portilla and Simoncelli 2000; Selesnick, Baraniuk, and Kingsbury 2005), the complex-valued nature of the bivariate $\mathbb{C}\mathrm{NLT}$ coefficients (see Section 2.3.2) provides both local phase and spectral information. To estimate the dependence between pairs of time series, we first define the complex lifting cross-periodogram, the cross-spectral analog of the periodogram. As in Section 2.3.2, our discussion will be split based on whether the data have been sampled over the same or different grids.

Bivariate time series observed on the same grid. For each time point x_k, k = 1, …, n and artificial scale ℓⁱ, i = 1, …, J*, define the complex lifting cross-periodogram (also referred to as $\mathbb{C}$NLT cross-periodogram) for series observed on the same grid as (9) $$\begin{array}{ccc}\hfill I_{x_k}^{(1,2)}\left({\ell }^i\right)& =& \frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{d}_{x_k}^{1,p}\overline{d_{x_k}^{2,p}},\hfill \end{array}$$(9) where $d_{x_k}^{1,p}={\lambda }_{x_k}^{1,p}+\mathrm{i}{\mu }_{x_k}^{1,p}$ and $d_{x_k}^{2,p}={\lambda }_{x_k}^{2,p}+\mathrm{i}{\mu }_{x_k}^{2,p}$ are the detail coefficients from f¹ and f². The $\mathbb{C}$NLT cross-periodogram consists of combinations of coefficients from each series and provides information about the relationship between the signals. Note that unlike the $\mathbb{C}$NLT periodogram, the cross-periodogram is complex-valued.

Similar to classical Fourier cross-spectrum methodology (see, e.g., Priestley 1983), the $\mathbb{C}$NLT cross-periodogram can be separated into its real and imaginary parts to define the $\mathbb{C}$NLT co-periodogram and the $\mathbb{C}$NLT quadrature periodogram, respectively, resulting in $$\begin{array}{ccc}\hfill c_{x_k}\left({\ell }^i\right)& =& \frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\lambda }_{x_k}^{1,p}{\lambda }_{x_k}^{2,p}+\frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\mu }_{x_k}^{1,p}{\mu }_{x_k}^{2,p},\hfill \\ \hfill q_{x_k}\left({\ell }^i\right)& =& \frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\mu }_{x_k}^{1,p}{\lambda }_{x_k}^{2,p}-\frac{1}{n_{i,k}}\sum _{p\in P_{i,k}}{\lambda }_{x_k}^{1,p}{\mu }_{x_k}^{2,p}.\hfill \end{array}$$These quantities, together with the individual lifting spectra of each process, can be used to calculate the measures of phase and coherence between the two series f¹ and f²: (10) $$\begin{array}{c}\hfill {\varphi }_{x_k}\left({\ell }^i\right)={tan}^{-1}\left(\frac{-q_{x_k}\left({\ell }^i\right)}{c_{x_k}\left({\ell }^i\right)}\right),\end{array}$$(10) (11) $$\begin{array}{c}\hfill {\rho }_{x_k}\left({\ell }^i\right)=\frac{\sqrt{c_{x_k}{\left({\ell }^i\right)}^2+q_{x_k}{\left({\ell }^i\right)}^2}}{\sqrt{I_{x_k}^{\left(1\right)}\left({\ell }^i\right)I_{x_k}^{\left(2\right)}\left({\ell }^i\right)}}.\end{array}$$(11) The $\mathbb{C}$NLT cross-periodogram provides a measure of the dependence between series, but its magnitude is affected by the individual $\mathbb{C}$NLT periodograms of the signals. Hence as in the regularly sampled setting, it is preferable to normalize this quantity, providing a coherence measure that satisfies $0\le {\rho }_{x_k}\left({\ell }^i\right)\le 1$ (as in (11)). This is similar to the coherence measure for regularly sampled signals introduced by Sanderson, Fryzlewicz, and Jones (2010). The $\mathbb{C}$NLT phase as defined in (10) provides an indication of any time lag between the signals. Several examples examining the coherence and phase between signal pairs are given in Section 3.3.

Bivariate time series observed on different grids. Closer to real data scenarios, we now consider time series that were sampled over different irregular grids, with one such real data example being discussed in Section 3.3.3. To obtain the cross-spectral quantities, we combine the appropriate sets of detail coefficients for each grid, corresponding to f¹ and f², that is, $D_{x^{\left(j\right)}}^1\left({\ell }^i\right)$ and $D_{x^{\left(j\right)}}^2\left({\ell }^i\right)$ introduced in Equations (7) and (8). For each artificial time period, x^(j), j = 1, …, T* and artificial scale ℓⁱ, i = 1, …, J*, we define the complex lifting cross-periodogram for series observed on different irregular grids as (12) $$\begin{array}{c}\hfill I_{x^{\left(j\right)}}^{(1,2)}\left({\ell }^i\right)=\frac{1}{n_{i,j}}\sum _{s=1}^{n_{i,j}}\text{order}{\left\{D_{x^{\left(j\right)}}^1\left({\ell }^i\right)\right\}}_s\text{order}{\left\{\overline{D_{x^{\left(j\right)}}^2\left({\ell }^i\right)}\right\}}_s,\end{array}$$(12) where n_{i, j} is the number of pairs in the grid square defined at time x^(j) and scale ℓⁱ, and order{D}_s indicates the sth time-ordered detail.

If the sampling schemes coincide for the two series ({x¹_k}_k ≡ {x²_k}_k) and the same trajectories are used to generate the details ${\left\{d_{x_k}^{1,p}\right\}}_{p,k}$, respectively, ${\left\{d_{x_k}^{2,p}\right\}}_{p,k}$, then Equations (9) and (12) coincide, except for the quantities being also averaged over the defined artificial time period. The co- and quadrature periodograms may be obtained in the same fashion as above, and subsequently used to yield the lifting phase and coherence in this setting.

Figures 2 and 3 provide a visual representation for the complex lifting cross-periodogram construction under the assumption of the same, respectively, different sampling grids.

We now make some remarks about the proposed periodogram constructions.

Scale interpretation. The relationship between artificial scale (ℓⁱ) and classical Fourier frequency can be described in terms of the scale, which maximizes the coherence for a Fourier wave of period T. Defining $\rho \left({\ell }^i\right)=\frac{1}{n}{\sum }_{k=1}^n{\rho }_{x_k}\left({\ell }^i\right)$, the design of the filters outlined in Section 2.2 is such that ${\ell }^i={argmax}_{j\in \{1,\cdots ,J^{*}\}}\rho \left({\ell }^j\right)=T/3$.

We emphasize that this relationship is dictated by the choice of filter pairs: the $\mathbb{C}$NLT periodogram and co-periodogram (as defined above) are composed of the sum of the wavelet coefficients from the two schemes, while the quadrature periodogram contains products of the coefficients. Hence, to ensure that the resulting estimates are interpretable, the two filters are specified so that combinations of coefficients (either through multiplication or summation) provide the same scale-frequency relationship (see Sanderson 2010, sec. 5.3 and 6.2.1). The provided mapping between wavelet lifting scale and Fourier frequency can be used to compare our results to those of classical Fourier-based methods (see Section 3.3 next).

Periodogram smoothing over time. As is customary, the $\mathbb{C}$NLT periodogram will be smoothed over time using simple moving average smoothing, that is, we compute ${\tilde{I}}_{x_k}\left({\ell }^i\right)=\frac{1}{\#\left(M_k^i\right)}{\sum }_{j\in M_k^i}{I}_{x_j}\left({\ell }^i\right)$, where Mⁱ_k = {j: x_k − Mⁱ < x_j ⩽ x_k + Mⁱ} and Mⁱ denotes the width of the averaging window, permitted to take different values for each scale, lⁱ.

3.3 Examples

We now illustrate the proposed methodology by application to both simulated and real irregular time series. The results were produced in the R statistical computing environment (R Core Team 2013), using modifications to the code from the adlift package (Nunes and Knight 2012) and the nlt package (Knight and Nunes 2012).

3.3.1 Simulated Data

Signals sampled on the same irregular sampling grid. In this example, the methods of Section 3.2 are applied to bivariate series observed on the same sampling grid: {(x_k, f¹_k, f²_k)}²⁰⁰_{k = 1}, where $$\begin{array}{ccc}\hfill f_k^1& =& sin\left(\frac{2\pi x_k}{10}\right)+sin\left(\frac{2\pi x_k}{30}\right)+sin\left(\frac{2\pi x_k}{70}\right)+{\zeta }_k^1,\hfill \\ \hfill f_k^2& =& sin\left(\frac{2\pi (x_k-\tau )}{30}\right)+{\zeta }_k^2,\hfill \end{array}$$where τ = 0 for x_k < 200 and τ = 6 for x_k ⩾ 200, and the quantities ζ¹_k and ζ²_k are independent, identically distributed Gaussian variables with mean zero and variance 0.2². The observations are irregularly sampled such that (x_{k + 1} − x_k) ∈ {n/10: n = 10, 11, …, 30} and $\frac{1}{(n-1)}{\sum }_{k=1}^{n-1}(x_{k+1}-x_k)=2$.

Estimates for coherence and phase are computed using the complex-valued lifting scheme using a sample of P = 1500 randomly sampled trajectories, discretizing using J* = 20 artificial scales and smoothing over time using a window of width Mⁱ = 60, ∀ i. The coherence estimate (Figure 4, right) provides a clear visualization of the dependence between the two series, with a peak occurring at scale log ₂(30/3) = 3.3 (equivalent to a Fourier period of 30). The time lag that is introduced halfway through the second signal is also clearly captured by the phase estimate (Figure 5, left), which is approximately zero for the first half of the series, then shows a marked increase for the second half.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 4. Coherence estimation for data observed on the same irregular sampling grid: using the real-valued bivariate lifting scheme (left); using the complex-valued lifting scheme (right). Scale gets coarser from bottom upward.

Display full size

Figure 4. Coherence estimation for data observed on the same irregular sampling grid: using the real-valued bivariate lifting scheme (left); using the complex-valued lifting scheme (right). Scale gets coarser from bottom upward.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 5. Phase estimation using the complex-valued lifting scheme: data observed on the same irregular sampling grid (left); data observed on different irregular sampling grids (right); data observed on the same regular grid (bottom). Scale gets coarser from bottom upward.

Display full size

Figure 5. Phase estimation using the complex-valued lifting scheme: data observed on the same irregular sampling grid (left); data observed on different irregular sampling grids (right); data observed on the same regular grid (bottom). Scale gets coarser from bottom upward.

For comparison, the estimated coherence using a real-valued bivariate scheme (Sanderson 2010) is also reported (Figure 4, left). It is interesting to note that although this method also clearly estimates a dependence for the first half of the series, it does not continue to detect it following the time delay. This again emphasizes the advantage of using a second filter, present in the complex-valued lifting transform.

Signals sampled on different irregular sampling grids. The methods described in Section 3.2 are now demonstrated by revisiting the same simulated data example, but with the two series observed on different irregularly spaced sampling grids: {(x¹_k, x²_k, f¹_k, f²_k)}²⁰⁰_{k = 1}. Aside from the sampling, the series satisfy the same properties as previously described.

The estimates were obtained using a discretization of J* = 15 artificial scales and T* = 60 artificial time points, while a smoothing window of width Mⁱ = 60 was applied at all scales as in the previous example. The resulting estimated coherence and phase are shown in Figure 6, respectively, Figure 5, right. It is interesting to note that while estimates broadly agree with those corresponding to sampling using the same (irregular) grid (Figure 4 and Figure 5, left), the price to pay for the different sampling schemes is the reduced clarity of the estimator. This is point is further reinforced by the phase estimate corresponding to a regular sampling situation, see Figure 5, bottom.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 6. Coherence estimation for data observed on different irregular sampling grids. Scale gets coarser from bottom upward.

Display full size

Figure 6. Coherence estimation for data observed on different irregular sampling grids. Scale gets coarser from bottom upward.

Coherence and phase analysis comparison with Fourier-based methods. For comparison with established Fourier-based techniques, we also performed coherence analysis of stationary, regularly sampled vector autoregressive (VAR) processes, as well as phase analysis of the signals described above. For regularly sampled stationary processes, we compared our estimates to the well-behaved Fourier estimates, while in the presence of sampling irregularity/nonstationarity, we compared our method to the short-time Fourier transform (STFT) and the Lomb–Scargle method. For brevity, we do not include the coherence and phase comparison plots here, but they can be found in Appendix A of the supplementary material.

Specifically, in the supplementary material we illustrate the coherence estimates obtained through both a classical Fourier-based approach and our lifting-based method on two bivariate VAR processes. The resulting estimates agree very well, with the lifting-based estimate displaying a slight depreciation when compared to the well-behaved Fourier estimates, suited for regular sampling and stationary process behavior. However, in general if the data are believed to be amenable to be analyzed with standard methodology, Fourier-based estimation should be preferred to the proposed method that was specifically designed to offer a solution for the challenging situations that include irregular sampling.

As already highlighted, traditional methods do not readily handle data that feature both potential nonstationarities and irregular sampling, thus STFT required further intervention while the Lomb–Scargle method failed to account for nonstationarity. Thus to obtain the desired phase analysis, we mapped the irregular data to a regular grid (by, e.g., interpolation) and then used STFT to capture the nonstationary time-frequency content of the data. The Lomb–Scargle analysis naturally dealt with the sampling irregularity, but assumed stationarity and therefore it did not provide time-localization information. The phase estimation plots of the STFT method exhibit little resolution in time or frequency, possibly due to the spectral blurring induced by the overlapping windows in the STFT as noted by Shumway and Stoffer (2013). Furthermore, for signals sampled over different irregular grids, the method creates additional blurring in the phase plot. By contrast, the Lomb–Scargle method is able to deal naturally with the irregular sampling structure of the signals, but it does not contain any time-phase information. In addition, there is no marked distinction in frequency where the phase is large, unlike for that of our complex lifting method (see Figure 5). These features yet again highlight the appeal of our technique.

3.3.2 Simulated Data With Varying Time Delay

The next example explores the effect of increasing the time delay between two signals. For each value of τ = 1, …, 15, the series {(x_k, f¹_k, f²_k)}²⁰⁰_{k = 1} are simulated following $$\begin{array}{ccc}\hfill f_k^1& =& sin\left(\frac{2\pi x_k}{30}\right)+{\zeta }_k^1,\hfill \\ \hfill f_k^2& =& sin\left(\frac{2\pi (x_k-\tau )}{30}\right)+{\zeta }_k^2,\hfill \end{array}$$where (x_{k + 1} − x_k) ∈ {n/10: n = 10, 11, …, 30} and $\frac{1}{(n-1)}{\sum }_{k=1}^{n-1}(x_{k+1}-x_k)=2$, ζ¹_k and ζ²_k are independent, identically distributed Gaussian variables with mean zero, variance 0.2².

Just as in the classical (Fourier) analysis, it is interesting to inspect the coherence and phase across frequencies (here, scales) to relate the common behavior of the two series and possible time delays, respectively. The estimated coherence and phase corresponding to the increasing τ = 1, …, 15 are shown in Figure 7. To give an overall sense of the coherence and phase magnitude over time, the estimates are averaged over the full time range to give $\rho \left({\ell }^i\right)=\frac{1}{200}{\sum }_{k=1}^{200}{\rho }_{x_k}\left({\ell }^i\right)$. We used P = 750 randomly sampled trajectories and discretized using J* = 20 artificial scales.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 7. (a) Coherence and (b) phase between f¹ and f² (Section 3.3.2) as a function of scale and $\tau \in \overline{0,15}$. For τ ≠ 0 the coherence and (absolute) phase are greatest at scale log ₂(30/3).

Display full size

Figure 7. (a) Coherence and (b) phase between f¹ and f² (Section 3.3.2) as a function of scale and $\tau \in \overline{0,15}$. For τ ≠ 0 the coherence and (absolute) phase are greatest at scale log ₂(30/3).

When τ = 0, the coherence is 1 and the phase is 0. For τ ≠ 0 the coherence is greatest at a scale of log ₂(30/3), corresponding to the period of variation (T = 30) in the data. The coherence intensity and response over scale are affected by the magnitude of the time delay. The coherence is lowest at time delays around 7.5 (T/4), and at these shifts the peak at scale log ₂(30/3) is also more pronounced. At τ = 15 (T/2), the signals are sign reversed versions of each other and, again, the observed coherence is 1 at all scales. The phase is also greatest at scale log ₂(30/3). The phase response varies as a function of time delay and alternates between positive and negative values, with |φ(ℓⁱ)| maximized at ℓⁱ = T/4. This is displayed in Figure 8, which shows the estimated phase at scale log ₂(30/3) as a function of time delay.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 8. Estimated phase between f¹ and f² (Section 3.3.2) at scale log ₂(30/3) (equivalent to a Fourier period of 30), as a function of τ.

Display full size

Figure 8. Estimated phase between f¹ and f² (Section 3.3.2) at scale log ₂(30/3) (equivalent to a Fourier period of 30), as a function of τ.

For completeness, we also provide a direct comparison with classical Fourier coherence and phase estimation when the signals are regularly sampled (see Figure 9). While the overall behavior is similar for both the classical and $\mathbb{C}\mathrm{NLT}$ methods, the Fourier method displays less variability across coherence estimates with the changing time delay (Figure 9, left), as well as more localized phase-frequency information (Figure 9, right). However, in general, note that if the data are believed to be amenable to be analyzed with standard methodology, this should be preferred to the proposed $\mathbb{C}\mathrm{NLT}$ method that was specifically designed to offer a solution for the challenging situations that include irregular sampling.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 9. (a) Coherence and (b) phase between f¹ and f² (Section 3.3.2) as a function of frequency and $\tau \in \overline{0,15}$ using classical Fourier methods.

Display full size

Figure 9. (a) Coherence and (b) phase between f¹ and f² (Section 3.3.2) as a function of frequency and $\tau \in \overline{0,15}$ using classical Fourier methods.

3.3.3 Financial Time Series

In this section, we demonstrate the use of the proposed complex-valued lifting transform through an application to financial data consisting of prices of all trades on 1 March 2011 (in normal trading hours) for two IT companies, Baidu and Google, both traded on the NASDAQ stock exchange. Comparison of the two companies is of interest as the main product of both is a search engine, but they are based in different geographical regions.

Often several trades per second occur and in this case the last quoted value for each second is selected. Thus, the finest sampling interval is one second, but as there are seconds with no trades, the time series are not equally spaced. For the analysis, we consider the returns of each series—for Google, the series contains 7984 observations with an average sampling distance of 2.93 sec and range 1 to 48; for Baidu, the series contains 6535 observations with an average sampling gap of 3.58 sec and range 1 to 52.

The data were analyzed using the methodology described in Section 3.2 using J* = 15 artificial scales and T* = 390 artificial time intervals (each time interval has a width of 60 sec). The estimates were smoothed over time using a window width of M¹ = 60 min at the finest scale and increasing by a factor of 1.05, to provide a larger smoothing window for each subsequent scale. The coherence estimate is shown in Figure 10(a) for scales up to 10. The main feature of the resulting coherence estimate is an increased coherence around scale 6, corresponding to a Fourier frequency of T ≈ 3 min. The magnitude of the coherence at this scale is seen to be more pronounced toward the end of the day. There is also a period of higher coherence observed in the middle of the day, at low wavelet scales (corresponding to high-frequency information).

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 10. Coherence between Google and Baidu using methods from Section 3.2: (a) computed on different irregular sampling grids; (b) computed using 1 min averages. Scale gets coarser from bottom upward.

Display full size

Figure 10. Coherence between Google and Baidu using methods from Section 3.2: (a) computed on different irregular sampling grids; (b) computed using 1 min averages. Scale gets coarser from bottom upward.

Complex-Valued Wavelet Lifting and Applications

All authors

Jean Hamilton http://orcid.org/0000-0003-3326-9842, Matthew A. Nunes http://orcid.org/0000-0002-4719-2690, Marina I. Knight http://orcid.org/0000-0001-9926-6092 & Piotr Fryzlewicz

https://doi.org/10.1080/00401706.2017.1281846

Published online:

26 May 2017

Figure 11. Ethanol data and estimates. Small circles = data; solid line = estimate. Top left: smoothing spline with cross-validated smoothing parameter; top right: multiple observation adaptive lifting using $\mathbb{R}$-lift with heteroscedastic variance computation and EbayesThresh posterior median thresholding; bottom: $\mathbb{C}$-AP1S with heteroscedastic variance computation and level-dependent soft thresholding.

Display full size

Figure 11. Ethanol data and estimates. Small circles = data; solid line = estimate. Top left: smoothing spline with cross-validated smoothing parameter; top right: multiple observation adaptive lifting using $\mathbb{R}$-lift with heteroscedastic variance computation and EbayesThresh posterior median thresholding; bottom: $\mathbb{C}$-AP1S with heteroscedastic variance computation and level-dependent soft thresholding.

One usual treatment of such irregular data would be to consider it in terms of the 1 min average returns. The estimated coherence using the aggregated data is shown in Figure 10(b), where J* = 10 artificial scales and T* = 78 artificial time intervals (each representing a range of 5 min) were used. Notice that finer behavior details are erased, reflecting the coarser sampling rate of the averaged data, and that spurious coherence is unsurprisingly induced by aggregation.