arXiv:1607.05167v2 [math.ST] 10 Aug 2017

Two-step wavelet-based estimation for mixed Gaussian fractional
processes ∗†‡
Patrice Abry
Laboratoire de Physique,
Université de Lyon,
ENS de Lyon,
Université Claude Bernard,
CNRS, F-69342 Lyon

Gustavo Didier
Mathematics Department
Tulane University

Hui Li
Mathematics Department
Tulane University
August 14, 2017

Abstract
A mixed Gaussian fractional process {Y (t)}t∈R = {P X(t)}t∈R is a multivariate stochastic process obtained by pre-multiplying a vector of independent, Gaussian fractional process
entries X by a nonsingular matrix P . It is interpreted that Y is observable, while X is a hidden process occurring in an (unknown) system of coordinates P . Mixed processes naturally
arise as approximations to solutions of physically relevant classes of multivariate fractional
SDEs under aggregation. We propose a semiparametric two-step wavelet-based method for
estimating both the demixing matrix P −1 and the memory parameters of X. The asymptotic normality of the estimators is established both in continuous and discrete time. Monte
Carlo experiments show that the finite sample estimation performance is comparable to that
of parametric methods, while being very computationally efficient. As applications, we model
a bivariate time series of annual tree ring width measurements, and establish the asymptotic
normality of the eigenstructure of sample wavelet matrices.

1

Introduction

Numerous data sets from a wide range of applications in science, technology and engineering have
been analyzed by means of fractional processes or models. Examples include natural systems
(hydrodynamic turbulence, Mandelbrot (1974); geophysics, Foufoula-Georgiou and Kumar (2014);
heart rate variability, Ivanov et al. (1999); infraslow – i.e., below 1Hz – brain activity, Ciuciu et al.
(2014)) and artificial systems (e.g., Internet traffic, Taqqu et al. (1997), Fontugne et al. (2017)).
Self-similar processes form a subclass of fractional processes that has been widely studied and
∗
The first author was partially supported by grant ANR-16-CE33-0020 MultiFracs. The second author was
partially supported by the prime award no. W911NF-14-1-0475 from the Biomathematics subdivision of the Army
Research Office, USA. The second author’s long term visits to ENS de Lyon were supported by the school.
†
AMS Subject classification. Primary: 62M10, 60G18, 42C40.
‡
Keywords and phrases: fractional stochastic process, multivariate, operator self-similarity, demixing, wavelets.

used in applications. A univariate stochastic processes Z = {Z(t)}t∈R is called self-similar when
it satisfies the scaling relation
L

{Z(ct)}t∈R = {cH Z(t)}t∈R ,

c > 0,

(1.1)

L

for some Hurst exponent H > 0, where = denotes the equality of finite dimensional distributions.
In particular, fractional Brownian motion (fBm) is the only Gaussian, self-similar, stationary
increment stochastic process (e.g., Embrechts and Maejima (2002), Taqqu (2003)). The probability theory and statistical methodology for univariate self-similar and related processes is now
voluminous (e.g., Mandelbrot and Van Ness (1968), Taqqu (1975, 1979), Dobrushin and Major
(1979), Granger and Joyeux (1980), Hosking (1981), Fox and Taqqu (1986), Dahlhaus (1989),
Beran (1994), Robinson (1995a, 1995b), Beran et al. (2013), Bardet and Tudor (2014), Clausel et
al. (2014b), Pipiras and Taqqu (2017), to name a few).
In modern applications, however, data sets are often multivariate, since several natural and
artificial systems are monitored by a large number of sensors. Accordingly, the literature on multivariate fractional processes has been expanding at a fast pace. The contributions include Hosoya
(1996, 1997), Lobato (1997), Marinucci and Robinson (2000), Shimotsu (2007), Becker-Kern and
Pap (2008), Robinson (2008), Hualde and Robinson (2010), Nielsen (2011), Sela and Hurvich
(2012) and Kechagias and Pipiras (2015a, 2015b), in the time and Fourier domains, and Wendt
et al. (2009), Amblard and Coeurjolly (2011), Amblard et al. (2012), Coeurjolly et al. (2013),
Achard and Gannaz (2016), Frecon et al. (2016), Abry and Didier (2017), in the wavelet domain
(see also Marinucci and Robinson (2001), Robinson and Yajima (2002), Nielsen and Frederiksen
(2011), Shimotsu (2012) on the related fractional cointegration literature in econometrics).
In this paper, we propose a new semiparametric statistical method for a subclass of multivariate fractional processes, i.e., those of the form
{Y (t)}t∈R = {P X(t)}t∈R ,

(1.2)

where P is a nonsingular matrix and
{X(t)}t∈R = {(X1 (t), . . . , Xn (t))T }t∈R

(1.3)

is a vector of independent Gaussian fractional processes. The process Y = {Y (t)}t∈R is assumed
observable. On the other hand, X = {X(t)}t∈R can be interpreted either as a hidden process
whose components get scrambled by a mixing matrix parameter P , or as one occuring in a different
system of coordinates (see Remark 2.4 on nonsquare matrices P ). One key statistical challenge is
to retrieve the fractional information (e.g., on Hurst exponents or memory parameters) contained
in X. If, for example, X is a vector of (independent) fBm entries
{X(t)}t∈R = {(Bh1 (t), . . . , Bhn (t))T }t∈R ,

0 < h1 ≤ . . . ≤ hn < 1,

(1.4)

where hi , i = 1, . . . , n, denote the individual Hurst exponents, then the univariate-like statistical
analysis of each entry of Y will often generate estimates that are undetermined convex combinations of Hurst exponents or, at large scales, estimates of the largest Hurst exponent (c.f. Abry
and Didier (2017), Introduction).
It has been shown (Tsai et al. (2017)) that processes of the form (1.2) naturally arise as
approximations to solutions of physically relevant classes of multivariate fractional SDEs under
aggregation (this is recapped in Section 2.1). In addition, it is well known that many real data sets
– e.g., tree ring widths, economic output, river flows, or rainfall – are obtained through aggregation
2

over a certain time interval, which points to the usefulness of the model (1.2). Multivariate
fractional processes of the form (1.2) are also closely related to the so-named operator self-similar
(o.s.s.) random processes and fields (Laha and Rohatgi (1981), Hudson and Mason (1982)), a
topic that has attracted much attention recently (e.g., Maejima and Mason (1994), Mason and
Xiao (2002), Biermé et al. (2007), Xiao (2009), Guo et al. (2009), Didier and Pipiras (2011, 2012),
Clausel and Vedel (2011, 2013), Li and Xiao (2011), Dogan et al. (2014), Puplinskaitė and Surgailis
(2015), Didier et al. (2017b, 2017a)). In the context of o.s.s. and related processes, the estimation
of the matrix P is itself of great interest, since it makes up the system of coordinates of the Hurst
matrix (see Example 2.1).
The class (1.2) further provides an extension to the framework of fractional processes of the sonamed mixed processes from the blind source separation literature in signal processing, the latter
being well-established in traditional settings such as that of ARMA-like signals (e.g., Belouchrani
et al. (1997), Cardoso (1998), Pham and Cardoso (2001), Moreau (2001), Yeredor (2002), Parra
and Sajda (2003), Stone (2004), Ziehe et al. (2004), Choi et al. (2005), O’Grady et al. (2005),
Fevotte and Godsill (2006), Li et al. (2009), Common and Jutten (2010)).
In the preliminary study Didier et al. (2015), presented without proofs, the hidden process X
is given by (1.4) and a demixing estimator is proposed for P that draws upon the diagonalization
of sample covariance matrices. In this paper, we consider the broad framework where each (independent) entry of X in (1.3) is a continuous time fractional process with stationary increments of
some order, possibly zero (i.e., X is stationary). In addition, it is not assumed that, entrywise, X
is exactly self-similar as in (1.1) (see (2.9), (2.10) and (2.13) and the discussion in Example 2.1).
We construct a semiparametric two-step wavelet-based method for the estimation of the demixing
matrix P −1 and the individual memory parameters d1 , . . . , dn that can be summed up as follows.
−1 by jointly diago(S1) demixing step (change of coordinates): generate an estimator Pd
j
nalizing two wavelet variance matrices (i.e., W (2 ) at two different octaves j; see (3.3)) of
the mixed process Y ;

(S2) memory parameter estimation step: estimate d1 , . . . , dn by applying univariate wavelet
−1 Y (Veitch and Abry (1999), Bardet
b = Pd
regression to each entry of the demixed process X
(2002), Moulines et al. (2007a, 2007b, 2008)).
The use of a wavelet framework has the benefit of computational efficiency (Daubechies (1992),
Mallat (1999)), while being a natural choice for stochastic systems with stationary increments
of arbitrary order. In fact, for a large enough number of vanishing moments Nψ (see (2.17)),
wavelet coefficients {D(2j , k)}k∈Z ∈ Rn are stationary in the shift parameter k at every octave
j (see (2.17), (3.1) and Remark 2.6). In addition, basing step (S1) on wavelet variance matrices
−1 is consistent and asymptotically normal (Theorem
of Y ensures that the demixing estimator Pd
3.2). The latter property does not generally hold for estimators based on sample covariance
matrices; indeed, it is in part a consequence of the quasi-decorrelation property of the wavelet
transform (Flandrin (1992), Wornell and Oppenheim (1992), Masry (1993), Bardet and Tudor
(2010), Clausel et al. (2014a)). The estimator of the vector of Hurst parameters generated at step
(S2) is also consistent and jointly asymptotically normal (Theorem 3.3). With a view toward
hypothesis testing, the consistency and asymptotic normality of the estimators generated at both
steps (S1) and (S2) are shown to hold under mild assumptions even in the presence of equal
Hurst parameters (Corollary 3.1). Moreover, under the more realistic assumption that Y in
(1.2) is observed in discrete time, the asymptotic properties of the proposed estimators do not
qualitatively change (Theorems 4.2, 4.3 and Corollary 4.1).

3

We conducted broad Monte Carlo experiments for instances where X is made up of independent fractional Brownian motion components. In dimension 4, the results show that the
performance of the proposed two-step estimation method is similar to that for univariate estimators of Hurst parameters over finite samples. Moreover, notwithstanding its semiparametric and
hence more general nature, the method’s performance is comparable to that of fully parametric
Whittle-type maximum likelihood estimation in terms of mean squared error, while bearing the
advantage of being computationally very fast. In addition, an application of the two-step method
to a bivariate data set from bristlecone pine tree rings from California shows that the latter can
be reasonably modeled by means of the mixed form (1.2).
It should be noted that the two-step nature of the estimation method makes it rather flexible.
Although step (S2), as proposed, involves applying entrywise a univariate wavelet estimator, in
principle the wavelet-based demixing technique in step (S1) can be combined with any other
univariate method such as Whittle, local Whittle or spectral log-regression estimation (see, for
instance, Bardet et al. (2003)).
This paper is organized as follows. In Section 2, we lay out the notation, assumptions and
theoretical background of the paper. Section 3 contains the main mathematical results of the
paper, including the properties of wavelet analysis, assuming measurements in continuous time.
In particular, in Sections 3.2 and 3.3, we construct steps (S1) and (S2) of the two-step estimation
method, respectively. In Section 4, we extend the two-step estimation method to the context of
discrete time measurements. Section 5 contains all Monte Carlo studies. In Section 6, we provide
two applications. We analyze and model the aforementioned tree ring data set, and establish the
asymptotic normality of the eigenstructure of the sample wavelet variance matrix at fixed scales,
which is of independent interest. All proofs can be found in the Appendix, together with auxiliary
results.

2

Preliminaries

The dimension of the mixed process Y is denoted by n ≥ 2 throughout the paper.
We shall use the following matrix notation. M (m, n, R) is the vector space of all m × n
real-valued matrices, whereas M (n, R) is a shorthand for M (n, n, R). GL(n, R) is the general
linear group (invertible matrices), O(n) is the orthogonal group of matrices O such that OO∗ =
I = O∗ O, where ∗ represents the matrix adjoint and T is reserved for vector transpose. S(n, R),
S≥0 (n, R) and S>0 (n, R) are, respectively, the space of symmetric, the cone of symmetric positive
semidefinite and the cone of symmetric positive definite matrices. The symbol 0 represents a
vector or matrix of zeroes. A block-diagonal matrix with main diagonal blocks P1 , . . . , Pn or m
times repeated diagonal block P is represented by
diag(P1 , . . . , Pn ),

diagm (P),

(2.1)

respectively. The symbol k · k represents a generic matrix or vector norm. The lp entrywise norm
of the matrix A is denoted by
kAklp = k(ai1 ,i2 ) i1 =1,...,m klp =
i2 =1,...,n

m X
n
X

|ai1 ,i2 |p

i1 =1 i2 =1

The Fourier transform of any function f ∈ L2 (R) is defined by
Z
fb(x) =
f (t)e−ixt dt.
R

4

1/p

.

(2.2)

For S = (si1 ,i2 )i1 ,i2 =1,...,n ∈ M (n, R), let
vecS (S) = (s11 , s21 , . . . , sn1 , s22 , s32 , . . . , sn2 , . . . , snn ),
vecD (S) = (s11 , s22 , . . . , snn ),

vec(S) = (s11 , . . . , sn1 , s12 , . . . , sn2 , . . . , snn ).

(2.3)

In other words, the operator vecS (·) vectorizes the lower triangular entries of S, vecD (·) vectorizes
the diagonal entries of S, and vec(·) vectorizes all the entries of S. Note that the expressions in
(2.3) are defined as row vectors; this will make the notation simpler in several statements. When
establishing bounds, C denotes a positive constant whose value can change from one inequality
to the next.

2.1

Aggregation and mixed processes

Recent work (Chan and Tsai (2010), Tsai et al. (2017)) has established the connection between
aggregation and the emergence of mixed processes. We sketch the basic idea for the reader’s
convenience. A natural multivariate extension of Langevin-type dynamics is given by the SDE
t ≥ 0,

dY (t) = ΦY (t)dt + ΣdBh (t),

−Φ, Σ ∈ S>0 (n, R),

(2.4)

where Bh (t) = (Bh1 (t), . . . , Bhn (t))T is a vector of independent fBm entries {Bhi (t)}t≥0 with
Hurst parameters
0 < hi < 1, i = 1, . . . , n.
(2.5)
The solution of (2.4) can be written a.s. as
Z

Φt

t

Y (t) = e Y (0) +

eΦ(t−u) ΣdBh (u),

t ≥ 0,

(2.6)

0

which generalizes the univariate fractional Ornstein-Uhlenbeck process (Cheridito et al. (2003),
Prakasa Rao (2010)). Consider the case where the continuous time process {Y (t)}t≥0 defined by
(2.6) is digitalized by aggregation over interval 4, i.e.,
Yz4 =

Z

z4

Y (u)du,

z ∈ N ∪ {0}.

(z−1)4

Then, as 4 → ∞,
L

−diag(4−h1 , . . . , 4−hn )Σ−1 ΦYz4 → (Bh1 (z) − Bh1 (z − 1), . . . , Bhn (z) − Bhn (z − 1))T ,

(2.7)

L

where → denotes convergence of the finite dimensional distributions. Therefore, for large 4, the
aggregate process Yz4 can be approximated by the mixed process
z ∈ N ∪ {0}.

Yez := P Xz ,

(2.8)

Recall that fractional Gaussian noise (fGn) is the increment process of fBm. In (2.8), Xz is a vector
of n independent fGn entries with Hurst parameters (2.5) and P = −Φ−1 Σdiag(4h1 , . . . , 4hn ).
Note that the process (2.8) is a particular case of (1.2), with the latter restricted to discrete time.

5

2.2

Assumptions

Unless otherwise stated, we will make the following assumptions on Y throughout the paper.
Assumptions (A1), (A2) and (A3) describe, respectively, the covariance structure of the hidden
process X, the conditions on the mixing matrix P and the regularity properties of high frequency
components.
Assumption (A1): the observed process has the mixed form (1.2), where Xi , i = 1, . . . , n, in (1.3)
is either a Ni -th (Ni ≥ 1) order (covariance) stationary process with harmonizable representation
n Z eitx − PNi −1 1 (itx)l
o
l=0 l!
−(di −Ni )
e
, Ni − 1/2 ≤ di < Ni + 1/2,
{Xi (t)}t∈R =
|x|
g
(x)
B(dx)
i
(ix)Ni
t∈R
R
(2.9)
or a (covariance) stationary process (i.e., Ni = 0) with harmonizable representation
nZ
o
eix − 1 −di
e
eitx
{Xi (t)}t∈R =
, −1/2 ≤ di < 1/2.
(2.10)
|x| gi (x)B(dx)
ix
t∈R
R
By convention, the so-named memory parameters are ordered as
−1/2 ≤ d1 < d2 < . . . < dn .

(2.11)

e
e
e
In (2.9) and (2.10), B(dx)
is a Gaussian random measure satisfying B(−dx)
= B(dx)
and
2
e
E|B(dx)| = dx.
Assumption (A2):
P ∈ GL(n, R),

kp·l k = 1,

pll ≥ 0,

l = 1, . . . , n.

(2.12)

Assumption (A3): the C-valued functions gi (x) in (2.9) and (2.10) are bounded and satisfy
||gi (x)|2 − |gi (0)|2 | < L|x|β ,

L > 0,

i = 1, . . . , n,

(2.13)

for any x ∈ (−δ, δ) for some small δ > 0. In (2.13), β ∈ (1, 2] and satisfies
β + 1 < 2d1 + 2α

(2.14)

α > 1.

(2.15)

for some

Example 2.1 If the high frequency functions gi (x) are constant and Ni − 1/2 < di < Ni + 1/2,
i = 1 . . . , n, then the observed process Y satisfies the so-named operator self-similarity property.
L
In other words, {Y (ct)}t∈R = {cH Y (t)}t∈R , c > 0, where H = P diag(h1 , . . . , hn )P −1 is the Hurst
matrix with Hurst eigenvalues
1
hi = di − , i = 1 . . . , n,
(2.16)
2
and cH is defined by the matrix exponential
exp{log c H} =

∞
X
(log c H)k
k=0

k!

.

If, in addition, Ni = 1, i = 1, . . . , n, then Y is an operator fractional Brownian motion, namely,
a Gaussian, operator self-similar, stationary increment process (Mason and Xiao (2002), Didier
and Pipiras (2011, 2012)).
6

Example 2.2 The framework provided by assumptions (A1–3) is quite general. For example,
one arbitrary entry Xi , i = 1, . . . , n, of the hidden process X can be a fBm, a fGn, or a fractional
Ornstein-Uhlenbeck process. These processes are associated, respectively, with the high frequency
C
function instances gi (x) ≡ C (Ni = 1), gi (x) ≡ C (Ni = 0), and gi (x) = eixix−1 λ+ix
(Ni = 0) for
Cix
d
−ix
−d
some λ > 0. The instance Ni = 0 and gi (x) = eix −1 |x| (1 − e ) 1[−π,π) corresponds, in discrete
time, to FARIMA(0, d, 0) (e.g., Taqqu (2003)).
Remark 2.1 For a fixed i, in the boundary cases di = Ni − 1/2 the finiteness of second moments
in (2.9) and (2.10) implies that the high frequency function gi (x) must decay fast enough as
x → ∞ so as to make up for the lack of integrability of the power law. On the range di < −1/2,
see Remark 4.3 in Section 4.2.
Remark 2.2 Note that the assumption (2.11) that memory parameters are pairwise distinct is
lifted in Section 3.4.
Remark 2.3 In (2.11), one incurs no loss of generality by assuming that the memory parameters
are disposed in ascending order. This fact can be easily illustrated in dimension n = 2. Suppose
that the mixed process has the form Y (t) = P (Xd2 (t), Xd1 (t))T , where Xdi (t), i = 1, 2, are
independent fractional processes defined in (2.9) or (2.10) with parameters d1 < d2 . Let


0 1
.
R=
1 0
Then, Y (t) = P R (Xd1 (t), Xd2 (t))T , whence P R can be treated as the mixing matrix with unit
vector columns.
Remark 2.4 Mathematically speaking, it is natural to ask how useful it is to consider the model
(1.2) with a full rank matrix P ∈ M (m, n, R), where m 6= n. However, both cases m > n and
m < n fall outside the scope of this paper. When m > n, the observed process Y is improper,
namely, its finite dimensional distributions are contained in a proper subspace of Rn for some
t 6= 0 (even if, in addition, the high frequency functions gi , i = 1, . . . , n, are constant, Y cannot be
operator self-similar: see Example 2.1 or Hudson and Mason (1982)). When m < n, the spectral
densities involved are potentially much more complicated, with added power laws. Either situation
calls for the construction of particular methods.
Remark 2.5 Assumption (A3) is typical in a semiparametric estimation setting (e.g., Robinson
(1995a) and Moulines et al. (2007b, 2008)). Note that larger values of β correspond to greater
smoothness of the functions gi , i = 1, . . . , n, around the origin.
In Section 3, we will implicitly make the following assumptions on the underlying wavelet
basis, hence they will be omitted from statements.
Assumption (W 1): ψ ∈ L1 (R) is a wavelet function, namely,
Z
Z
2
ψ (t)dt = 1,
tq ψ(t)dt = 0, q = 0, 1, . . . , Nψ − 1,
R

Nψ ≥ Nn + 1,

(2.17)

R

for some number Nψ of vanishing moments, where Nn is as in (2.9) or (2.10).
Assumption (W 2):
supp(ψ) is a compact interval.
7

(2.18)

Assumption (W 3): for α as in (2.15),
b
sup |ψ(x)|(1
+ |x|)α < ∞.

(2.19)

x∈R

b
Under (2.17), (2.18) and (2.19), ψ is continuous, ψ(x)
is everywhere differentiable and its first
Nψ − 1 derivatives are zero at x = 0 (see Mallat (1999), Theorem 6.1 and the proof of Theorem
7.4). The condition (W 1) is equivalent to asserting that the first Nψ − 1 derivatives of ψb vanish
at the origin. This implies, using a Taylor expansion, that
|ψb(l) (x)| = O(|x|Nψ −l ),

l = 0, 1 . . . , Nψ ,

x → 0.

(2.20)

Example 2.3 If ψ is a Daubechies wavelet with Nψ vanishing moments, supp(ψ) = [0, 2Nψ − 1]
(see Mallat (1999), Proposition 7.4).
Remark 2.6 Assumption (W 1) requires using a number of vanishing moments Nψ larger than
the unknown integration order Nn . In practice, though, the latter parameter is rarely greater
than 2, so the requirement is easily met even for low values of Nψ .
Remark 2.7 Section 4, on measurements in discrete time, requires a slightly different set of
assumptions on the wavelet basis (see Section 4.1).

3

Wavelet-based estimation: continuous time

In Section 3.1, we establish basic as well as the asymptotic properties of the wavelet transform
of the process Y at fixed scales. Sections 3.2 and 3.3 contain the main mathematical results of
the paper. In the former and in the latter, respectively, the demixing step (S1) and the postdemixing Hurst parameter estimation step (S2) are laid out in full detail, and their asymptotic
properties are shown. Note that (S1) only involves wavelet analysis at fixed scales, while (S2)
generally requires taking a coarse scale limit a(ν)2j → ∞, due to the lack of exact self-similarity in
(2.9) and (2.10). Recall that, throughout this section, we are implicitly assuming that conditions
(W 1–3) hold.

3.1

Wavelet analysis at fixed scales: properties and asymptotic theory

For a wavelet function ψ ∈ L2 (R) with a number Nψ of vanishing moments, the vector wavelet
transform of Y is naturally defined as
Z
n
j
R 3 D(2 , k) =
2−j/2 ψ(2−j t − k)Y (t)dt, j ∈ N ∪ {0}, k ∈ Z,
(3.1)
R

provided the integral in (3.1) exists in an appropriate sense. It will be convenient to make the
change of variable z = 2−j t − k, and reexpress
Z
j
j/2
D(2 , k) = 2
ψ(z)Y (2j z + 2j k)dz.
R

The wavelet domain process {D(2j , k)}k∈Z is stationary in k (Proposition 3.1). The wavelet
spectrum (variance) at scale j is the positive definite matrix
ED(2j , k)D(2j , k)∗ = ED(2j , 0)D(2j , 0)∗ =: EW (2j ),
8

(3.2)

and its natural estimator, the sample wavelet variance, is the random matrix
Kj
1 X
W (2 ) =
D(2j , k)D(2j , k)∗ ,
Kj
j

Kj =

k=1

ν
,
2j

j = j1 , . . . , j m ,

(3.3)

for a total of
ν available (wavelet) data points.

(3.4)

The next proposition describes some properties of the wavelet coefficients (3.1) as well as the
general form of the wavelet spectrum (3.2).
Proposition 3.1 Under the assumptions (A1 − 2), let D(2j , k) and EW (2j , k) be as in (3.1) and
(3.2), respectively. Then,
(P 1) the wavelet transform (3.1) is well-defined in the mean square sense, and ED(2j , k) = 0;
d

(P 2) (stationarity for a fixed scale) {D(2j , k + h)}k∈Z = {D(2j , k)}k∈Z , h ∈ Z;
(P 3) the wavelet spectrum (3.2) can be expressed as

Z
 
x
∗
∗
2
b
EW (2j ) = 2jD
|ψ(x)|
|x|−D G j |x|−D dx 2jD .
2
R

(3.5)

In (3.5),
G(x) = P diag(|g1∗ (x)|2 , . . . , gn∗ (x)|2 )P ∗ ,

(3.6)

D = P diag(d1 , . . . , dn )P −1 ,

(3.7)

where, in (3.6),
(
gi∗ (x) =

gi (x) sin(x/2)
x/2 , di < 1/2;
gi (x),
di ≥ 1/2,

i = 1 . . . , n;

(P 4) the wavelet spectrum has full rank, namely, det EW (2j ) 6= 0, j ∈ N;
By a standard calculation, the wavelet variance (3.5) can be recast as
EW (2j ) = P E(2j )1/2 diag(22jd1 , . . . , 22jdn )E(2j )1/2 P ∗ ,

(3.8)

where
j

Z

E(2 ) = diag
R

2
b
|ψ(y)|
|y|−2d1 g1∗



y
2j



2

Z
dy, . . . ,

2
b
|ψ(y)|
|y|−2dn gn∗

R



y
2j



2


dy .

(3.9)

The following theorem establishes the asymptotic distribution of the vectorized sample wavelet
spectrum at a fixed set of octaves.
Theorem 3.1 Suppose Y = {Y (t)}t∈R satisfies the assumptions (A1 – 3). Let j1 < . . . < jm be
a fixed set of octaves. Then,
p
T
d
Kj (vecS (W (2j ) − EW (2j ))
→ N n(n+1) ×m (0, F ), ν → ∞
(3.10)
j=j1 ,...,jm

2

(see (2.3) on the notation vecS ). In (3.10), the matrix F ∈ S( n(n+1)
m, R) has the form F =
2
0
0
0
(Gjj )j,j =1,...,m , where each block Gjj ∈ M (n(n + 1)/2, R) is described in Proposition B.1.
9

3.2

Wavelet-based demixing (step (S1))

The joint diagonalization of two matrices is a well-known problem. For the case of symmetric
matrices, its description and full characterization can be stated as follows (see Theorem 4.5.17,
(b), in Horn and Johnson (1985)). Suppose C0 and C1 are symmetric and C0 is nonsingular.
Then, there are a nonsingular S ∈ M (n, R) and complex diagonal matrices Λ0 and Λ1 such that
C0 = SΛ0 S ∗ ,

C1 = SΛ1 S ∗ ,

(3.11)

if and only if the matrix C0−1 C1 is diagonalizable (in its Jordan form). In light of this, we can
cast a joint diagonalization algorithm in the form of pseudocode.
Pseudocode for exact joint diagonalization (EJD)
Input: C0 , C1 are symmetric matrices and the former is positive definite;
−1/2

Step 1: set W = C0

so that C0−1 = W ∗ W ;

Step 2: compute Q ∈ O(n) in the spectral decomposition W C1 W ∗ = Q∗ D1 Q;
Step 3: compute the demixing matrix B := QW ;
Step 4: stop and exit.

Example 3.1 In view of (3.8), it is clear that C0 = EW (2J1 ), C1 = EW (2J2 ), J1 < J2 , can
be jointly diagonalized, where the underlying process is defined in (1.2) under the assumptions
(A1-2). In addition,


C0−1 C1 = (P ∗ )−1 diag(22(J2 −J1 ) d1 , . . . , 22(J2 −J1 ) dn )E(2J1 )−1 E(2J2 ) P ∗ .
This expression constitutes a diagonal Jordan decomposition, whence (3.11) holds.
Remark 3.1 Steps 1–4 of the EJD algorithm should not be confused with steps (S1) and (S2)
of the proposed wavelet-based estimation method).
bν of a demixing matrix is defined next.
The proposed wavelet-based estimator B
Definition 3.1 ((S1) demixing step, continuous time) Consider two octaves 0 ≤ J1 < J2
for which
diag(22(J2 −J1 ) d1 , . . . , 22(J2 −J1 ) dn )E(2J1 )−1 E(2J2 )

has pairwise distinct diagonal entries. (3.12)

bν is the output of the EJD algorithm when
For ν ∈ N, the wavelet-based demixing estimator B
setting
C0 = W (2J1 ) and C1 = W (2J2 ).
(3.13)
In Theorem 3.2, stated next, we establish the consistency and asymptotic normality of the
estimator put forward in Definition 3.1. The result involves characterizing the set of solutions
provided by the EJD algorithm. In view of (3.8), this relies on reexpressing
1/2

(C0 =) EW (2J1 ) = P E(2J1 )

1/2

diag(22J1 d1 , 22J1 d2 , . . . , 22J1 dn )E(2J1 )
10

P ∗ =: RR∗ ,

1/2

(C1 =) EW (2J2 ) = P E(2J2 )

1/2

diag(22J2 d1 , 22J2 d2 , . . . , 22J2 dn )E(2J2 )

P ∗ =: RΛR∗ ,

(3.14)

where
1/2

Λ := diag(22(J2 −J1 ) d1 , . . . , 22(J2 −J1 ) dn )E(2J1 )−1 E(2J2 ),
(3.15)
and then making use of the matrix polar decomposition of R. Then, consistency and asymptotic
normality stem from obtaining the behavior of the sample counterparts W (2J1 ) and W (2J2 ) vis-àvis (3.14) by means of Proposition B.1 and Theorem E.1, plus the Delta method when developing
limits in distribution.
R := P E(2J1 )

diag(2J1 d1 , . . . , 2J1 dn ),

Theorem 3.2 For j ∈ N, let E(2j ) be as in (3.9). Also let
I = {Π ∈ M (n, R) : Π has the form diag(±1, . . . , ±1)}.

(3.16)

(i) Then,
MEJD = {Π diag(2−J1 d1 , . . . , 2−J1 dn )E(2J1 )−1/2 P −1 , Π ∈ I}

(3.17)

is the set of matrix solutions produced by the EJD algorithm when setting
C0 = EW (2J1 ) and C1 = EW (2J2 );

(3.18)

bν }ν∈N and some
(ii) in addition, assume condition (3.12) holds. For some estimator sequence {B
matrix Π ∈ I,
P
bν →
B
Π diag(2−J1 d1 , . . . , 2−J1 dn )E(2J1 )−1/2 P −1 ,

ν → ∞;

bν }ν∈N as described in (ii) satisfies
(iii) an estimator sequence {B
√
d
bν − Π diag(2−J1 d1 , . . . , 2−J1 dn )E(2J1 )−1/2 P −1 ))T →
ν(vec(B
N (0, ΣF (J1 , J2 ))

(3.19)

(3.20)

for some matrix Π ∈ I, where the covariance matrix ΣF (J1 , J2 ) is a function of F , and F
is defined in Theorem 3.1, with m = 2.
Remark 3.2 Note that, for J1 < J2 , since E(2J1 )−1 E(2J2 ) → I as J1 , J2 → ∞, then, under
(2.11), condition (3.12) always holds for large enough J1 , J2 .
Remark 3.3 As shown in the proof of Theorem 3.2, in (3.19) and (3.20), the factor Π ∈ I
stems from the spectral decomposition W C1 W ∗ = Q∗ D1 Q in Step 2 of the EJD algorithm. A
convenient choice is a matrix Π such that the main diagonal entries of Q ∈ O(n) are all positive.
bν−1 has a limit in probability of the form B −1 :=
Remark 3.4 By (3.19), any sequence B
P diag(β1 , . . . , βn ), |βi | =
6 0, i = 1, . . . , n, i.e., involving a non-identifiability factor post-multiplying
the mixing matrix P . However, note that D = P diag(d1 , . . . , dn )P −1 = B −1 diag(d1 , . . . , dn )B,
i.e., the columns of B −1 consist of (non-unit) eigenvectors of the memory matrix D. Consequently,
b := B
bν−1 diag(db1 , . . . , dbn )B
bν is a natural estimator of the latter, where db1 , . . . , dbn are univariate
D
(e.g., wavelet-based) estimators of the individual Hurst exponents obtained from the demixed
process.
Nevertheless, producing a direct estimator of P is straightforward. Just normalize each column
b −1 and multiply it by −1 if necessary as to arrive at a matrix Pb with
of the matrix estimator B
ν
positive diagonal entries (cf. (2.12)). This procedure is used in Section 5.
Remark 3.5 More precisely, the covariance matrix in the limit (3.20) can be written as
ΣF (J1 , J2 ) = A3 Σ2 A∗3 , where Σ2 and A3 are given by expressions (B.25) and (B.26), respectively. It is clear that the expression for ΣF (J1 , J2 ) is quite intricate, and the construction of
theoretical confidence intervals is a matter for future investigation (cf. Wendt et al. (2017)).
11

3.3

Wavelet-based estimation of memory parameters after demixing/changing
the coordinates (step (S2))

Throughout this section, a scaling factor a(ν) is assumed to be a dyadic sequence such that
a(ν)
ν
+
→ 0,
ν
a(ν)1+2β

ν→∞

(3.21)

where β satisfies (2.14) (see Remark 3.7 below on the choice of a(ν) in practice).
We start off with the output of step (S1) of the proposed two-step method (Section 3.2). Let
bν be the demixing matrix described in (3.19). Then, the demixed process is defined by
B
t ∈ R,

b
bν Y (t),
X(t)
:= B

(3.22)

of which only ν (wavelet) data points are available (c.f. (3.4)). For j ∈ N, let
WXb (a(ν)2j ),

EWX (a(ν)2j ),

(3.23)

b and the wavelet variance of the hidden process X, respecbe the sample wavelet variance of X
tively. Proposition B.2 in the Appendix establishes the asymptotic normality of WXb (a(ν)2j ) when
centered at EWX (a(ν)2j ). So, we are now in a position to define an estimator for the vector of
memory parameters dT = (d1 , . . . , dn ) of the hidden process X.
Definition 3.2 ((S2) Memory parameter estimation step, continuous time) Let
WXb (·)ii0 ,

EWX (·)ii0 ,

i, i0 = 1, . . . , n,

(3.24)

be the (i, i0 )-th entries of the matrices WXb (·) and EWX (·), respectively. Consider the regression
weight vectors
i T
wi = (w1i , . . . , wm
) ,
(3.25)
where

m
X
l=1

wli = 0,

2

m
X

jl wli = 1,

i = 1, . . . , n.

(3.26)

l=1

The wavelet-based estimator of the memory parameters d1 , . . . , dn in (2.11) is obtained by regressing the main diagonal terms WX (a(ν)2j )ii on the scale indices a(ν)2j , j = j1 , . . . , jm , i.e.,


 Pm

jl
1
db1
b (a(ν)2 )11 )
l=1 wl log2 (WX


. 
..
b=
(3.27)
d
.
 ..  := 
.
Pm
j
n
l
dbn
b (a(ν)2 )nn )
l=1 wl log2 (WX
b is provided in the following theorem.
The asymptotic distribution of the estimator d
b T = (db1 , . . . , dbn ) be the estimator defined by (3.27). Then,
Theorem 3.3 Let d


 
db1
d1


r
ν
 ..   ..  d
 .  −  .  → N (0, W), ν → ∞.
a(ν)
dn
dbn
In (3.28),
W = diag((w1 )T V (h1 )w1 , . . . , (wn )T V (hn )wn ),
12

(3.28)

the weight vectors wi , i = 1, . . . , n satisfy (3.26), and the matrix V (h) = {Vk1 ,k2 (h)}k1 ,k2 =1,...,m is
defined entrywise by
Vk1 ,k2 (d) =
where K(d) =

R

R |ψ(x)|

b

4πb4d−1
jk ,jk

 jk1 x 
−4d b 2
x
ψ
bjk1 ,jk2
22(jk1 +jk2 )d K 2 (d) R

2 |x|−2d dx

1

Z

2

2

 2jk2 x  2
dx,
ψb
bjk1 ,jk2

(3.29)

and bjk1 ,jk2 = gcd(2jk1 , 2jk2 ).

Remark 3.6 Theorem 3.3 shows that the individual memory estimators db1 , . . . , dbn are asympb estimated from the demixed
totically independent. In fact, the joint asymptotic distribution of d,
b
process X, is equal to that of the joint entrywise wavelet-based estimators of d1 , . . . , dn obtained
from the hidden process X (see Remark B.1). In other words, asymptotically, the demixing step
(S1) washes out the effect of the mixing matrix P on the estimation procedure.
Remark 3.7 In practice, the choice of a(ν) involves a statistical compromise. A large value
of a(ν) with respect to ν implies a relatively small bias, but also a relatively large variance.
Simulation results suggest the ratio ν/a(ν)2j should be no less than 23 .
Remark 3.8 Removing the condition (2.11) can alter the limits (3.28). For example, suppose
there are two blocks of equal memory parameters
d1 = . . . = dn1 < dn1 +1 = . . . dn ,

n1 , n − n1 ≥ 2,

and the high frequency functions gi (x) are identically constant for i = 1, . . . , n. Then, in Step 2
P
cC
b1 W
c∗ →
of the EJD algorithm, W
diag(22d1 (J2 −J1 ) In1 , 22dn (J2 −J1 ) In−n1 ). Thus, the eigenvectors
cC
b1 W
c ∗ do not have a limit in probability. In this case, the demixed process takes the form
of W
b
b
X(t)
= ADX(t)
(see expression
(B.29) for the definition of the matrix D), where the random
!
c
c
A1 A2
b=
matrix A
satisfies
c3 A
c4
A
c1 ∈ O(n1 ),
A

c4 ∈ O(n − n1 ),
A

√
c2 = OP (1/ ν),
M (n1 , n − n1 , R) 3 A

√
c3 = OP (1/ ν),
M (n − n1 , n1 , R) 3 A

c1 and A
c4 do not have a limit in probability. Therefore, we can write
and A
!
!
c1 X1 (t) + A
c2 X2 (t)
c1 X1 (t)
A
A
b
X(t)
=
=
c3 X1 (t) + A
c4 X2 (t)
c4 X2 (t) + oP (1),
A
A
where



X1 (t) = 

D(1, 1)X1 (t)
..
.
D(n1 , n1 )Xn1 (t)



,




D(n1 + 1, n1 + 1)Xn1 +1 (t)


..
X2 (t) = 
,
.
D(n, n)Xn (t)

c1 X1 (t)
and D(i, i) is the i-th diagonal entry of D, i = 1, . . . , n. Thus, each entry of the processes A
c4 X2 (t) has memory parameter d1 and dn , respectively. Even though we cannot retrieve the
and A
mixing matrix, we can still estimate the memory parameters and obtain an asymptotically normal
distribution. However, corresponding to each block of parameters, the estimators among each set
db1 , . . . , dbn1 and dbn1 +1 , . . . , dbn are asymptotically dependent (though independent across sets).
13

3.4

On the case of blocks of equal memory parameters

With a view toward hypothesis testing, we also consider the case where some, or all, memory
parameters d1 , . . . , dn are equal. In light of Remark 3.8, we will need make some change to our
assumptions. However, to attain consistency and asymptotic normality in steps (S1) and (S2),
it suffices to add minor constraints on the high frequency functions gi (x), i = 1, . . . , n, and hence
replace (A1) and (A3) with the following assumptions.
Assumption (A10 ): the observed process Y has the mixed form (1.2), where each component
Xi , i = 1, . . . , n, of the hidden process in (1.3) has the form (2.9) or (2.10), and the memory
parameters can be ordered as
−1/2 < d1 = . . . = dn1 < dn1 +1 = . . . = dn2 < . . . < dnp +1 = . . . = dn .
Assumption (A30 ): In addition to satisfying (A3), the high frequency functions gi (x), i =
1, . . . , n, are such that the matrix diag(22(J2 −J1 ) d1 , . . . , 22(J2 −J1 ) dn )E(2J1 )−1 E(2J2 ) has pairwise
distinct diagonal entries.
Corollary 3.1 Suppose the mixed process Y satisfies assumptions (A10 ), (A2) and (A30 ). Then,
the conclusions of Theorem 3.1, Theorem 3.2 and Theorem 3.3 hold.

4

Wavelet-based estimation: discrete time

In practice, only observations in discrete time are available, which renders the computation of
the theoretical wavelet coefficients D(2j , k) impossible. In this section, we study the asymptotic
performance of the two-step wavelet-based methodology under the assumption that only ν wavelet
data points from a discrete time sample
{Y (k)}k∈Z

(4.1)

of (1.2) are available (c.f. (3.4)). In Section 4.1, we lay out the notation and assumptions.
In Section 4.2, we develop the asymptotic distribution of the two-step wavelet-based method
estimators.

4.1

Notation and assumptions

Throughout this section, we suppose the wavelet approximation coefficients stem from Mallat’s
pyramidal algorithm, under a multiresolution analysis of L2 (R) (MRA; see Mallat (1999), chapter
7). Accordingly, we need to replace (W 2) with the following more restrictive condition.
Assumption (W 20 ): the scaling and wavelet functions ϕ ∈ L2 (R) and ψ ∈ L2 (R), respectively,
are compactly supported, integrable and ϕ(0)
b
= 1.
We also add the following condition.
Assumption (W 4): the function
X

k m ϕ(· − k)

k∈Z

is a polynomial of degree m for all m = 0, . . . , Nψ − 1.

14

Remark 4.1 The Daubechies scaling and wavelet functions satisfy (W 1), (W 20 ) and (W 3-4)
(Moulines et al. (2008), page 1927).
Throughout this section, we assume that the conditions (W 1), (W 20 ) and (W 3 − 4) hold. In
particular, conditions (W 1) and (W 4) imply that
Z
X
ψ(2−j t)
ϕ(t + l)lm dt = 0, j ≥ 0, m = 0, . . . , Nψ − 1.
(4.2)
R

4.2

l∈Z

Asymptotic theory for the two-step wavelet-based method (steps (S1) and
(S2))

Given (4.1), we initialize the algorithm with the vector-valued sequence
Rn 3 e
a0,k := Y (k),

k ∈ Z,

also called the approximation coefficients at scale 20 = 1. At coarser scales 2j , Mallat’s algorithm
is characterized by the iterative procedure
X
X
e
aj+1,k =
hk0 −2k e
aj,k0 , dej+1,k =
gk0 −2k e
aj,k0 , j ∈ N, k ∈ Z,
k0 ∈Z

k0 ∈Z

where the filter sequences {hk }k∈Z , {gk }k∈Z are called low- and high-pass MRA filters, respectively. Due to (W 20 ), only a finite number of filter terms is non-zero, which is convenient for
computational purposes (Daubechies (1992)). The normalized wavelet coefficients are defined by
e j , k) := 2−j/2 dej,k .
Rn 3 D(2
Let
f (2j ) = ED(2
e j , 0)D(2
e j , 0)∗ ,
EW

(4.3)

Kj
1 Xe j
j
f
e j , k)∗
W (2 ) =
D(2 , k)D(2
Kj

(4.4)

k=1

be the wavelet variance matrix and its sample counterpart, respectively, where Kj is as in (3.3).
The following theorem is the discrete time analogue of Theorem 3.1 and establishes the asymptotic
distribution of the wavelet variance matrices at fixed octaves.
Theorem 4.1 Let {Y (k)}k∈Z be the sequence (4.1). Let j1 < . . . < jm be a fixed set of octaves.
Then,
p
T
d
f (2j ) − EW
f (2j )))
Kj (vecS (W
→ N n(n+1) ×m (0, Fe),
(4.5)
j=j1 ,...,jm

2

as ν → ∞ (see (2.3) on the notation vecS ). In (4.5), the matrix Fe ∈ S( n(n+1)
m, R) has the form
2
e jj 0 )j,j 0 =1,...,m , where each block G
e jj 0 ∈ M (n(n + 1)/2, R) is described in Proposition B.1.
Fe = (G
f (2j ) can be recast as
Note that EW

where

f (2j ) = P Λ
e j P ∗,
EW

(4.6)

Z

Z
2
−2d1 ∗
2
2
−2dn ∗
2
e
Λj = diag
|Hj (x)| |x|
|g1 (x)| dx, . . . , |Hj (x)| |x|
|gn (x)| dx

(4.7)

R

R

(see Proposition C.1 in the Appendix). As in continuous time, expression (4.6) indicates that an
eν of P −1 can be generated by jointly diagonalizing W
f (2J1 ) and W
f (2J2 ), for J1 6= J2 .
estimator B
15

Definition 4.1 ((S1) demixing step, discrete time) Consider two octaves 0 ≤ J1 < J2 for
which
e −1
eJ Λ
(4.8)
Λ
2 J1 has pairwise distinct diagonal entries.
eν is the output of the EJD algorithm when
For ν ∈ N, the wavelet-based demixing estimator B
setting
f (2J1 ) and C1 = W
f (2J2 ).
C0 = W
(4.9)
As a consequence of Theorem 4.1 and by following the same argument as in the proof of
eν .
Theorem 3.2, we obtain the limiting distribution of B
Theorem 4.2 Assume condition (4.8) holds. Then,
√
d
eν − Π Λ
e −1/2 P −1 ))T →
N (0, ΣFe (J1 , J2 )),
ν(vec(B
J1

(4.10)

e J is defined by (4.7), for some matrix
where Λ
1
Π ∈ {Π ∈ M (n, R) : Π has the form diag(±1, . . . , ±1)}.
In (4.10), the covariance matrix ΣFe (J1 , J2 ) is a function of Fe, and Fe is defined in Theorem 4.1
with m = 2.
Remark 4.2 As in continuous time (see Remark 3.2), the condition (4.8) is not restrictive (see
Proposition C.4).
Let
k ∈ Z,

e
eν Y (k),
X(k)
:= B

(4.11)

be the demixed process, of which ν (wavelet) data points are available (see (4.1)). As with its
f e (a(ν)2j ) is
continuous time counterpart WXb (a(ν)2j ) (see (3.23)), the sample wavelet variance W
X
e W
fX (a(ν)2j ) (see Proposition C.3 in the
asymptotically normal when centered at the matrix DE
e We are now in a position to define
Appendix, and also expression (C.9) for the definition of D).
T
the estimators of the memory parameters d = (d1 , . . . , dn ).
Definition 4.2 ((S2) Memory parameter estimation step, discrete time) For i, i0 =
f e (a(ν)2j ). The
1, . . . , n, let WXb (a(ν)2j )ii0 be the (i, i0 )-th entry of the sample wavelet variance W
X
wavelet-based estimator of the memory parameters d1 , . . . , dn in (2.11) is obtained by regressing
the terms WXb (a(ν)2j )ii on the scale indices a(ν)2j , j = j1 , . . . , jm , i.e.,


 Pm
1
f e (a(ν)2jl )11 )
de1
l=1 wl log2 (WX


..
. 
e=
d
,
 ..  := 
.
Pm
n
j
e
f e (a(ν)2 l )nn )
dn
w log (W


l=1

l

2

X

i )T , i = 1, . . . , n, satisfy (3.26).
where the weight vectors wi = (w1i , . . . , wm

e is established.
In the following theorem, the asymptotic normality of the estimator d

16

(4.12)

e T = (de1 , . . . , den ) be the estimator defined by (4.12). Suppose the scaling
Theorem 4.3 Let d
factor a(ν) satisfies
a(ν)
ν
→ 0, ν → ∞,
+
ν
a(ν)1+2β∗
where


β∗ =

min{β, 2d1 + 2}, −1/2 < d1 < 1/2;
min{β, 2d1 },
d1 ≥ 1/2.

(4.13)

Then,


r

 
de1
ν
 ..  
 . −
a(ν)
den



d1

d
..  →
f
N (0, W),
. 
dn

ν → ∞,

(4.14)

where
f = diag((w1 )T Ve (d1 )w1 , . . . , (wn )T Ve (dn )wn ),
W
the weight vectors wi , i = 1, . . . , n satisfy (3.26), Ve (h) is a m × m matrix whose (l, l0 )-th entry is
min(jl ,jl0 )

4π22d|jl −jl0 |2
Vel,l0 (d) =
K(d)

Z
|x|<π

D|jl −jl0 | (x, d) is defined in (C.12), and K(d) =

R

|D|jl −jl0 | (x; d)|2 dx,

l, l0 = 1, . . . , m,

2
−2d dx.
b
R |ψ(x)| |x|

The next result is the discrete time analogue of Corollary 3.1, i.e., for the case where some,
or all, memory parameters are equal. Note that the assumptions on the process Y do not change
from continuous to discrete time.
Corollary 4.1 Suppose the underlying mixed process Y (k) satisfies (A10 ), (A2) and (A30 ). Then,
the conclusions of Theorem 4.1, Theorem 4.2 and Theorem 4.3 hold.
Remark 4.3 All results in continuous time hold if we allow the components Xi (t) to be stationary
with memory parameter di < −1/2. In discrete time, all results hold if we assume that supp gi =
[−π, π), since in this case the second and the third terms, respectively, in expressions (C.13) and
(C.15) are identically zero.

5
5.1

Monte Carlo studies
Performance over finite samples

We studied the performance of the two-step wavelet-based method over finite samples assuming
the hidden process X is made up of 4 independent fractional Brownian motion components observed in discrete time. For notational simplicity, denote X := Bh , Y := BH (see Example 2.1).
Recall that, in this case, the relation (2.16) holds between the memory parameters and the individual Hurst exponents. We simulated R = 500 sample paths of with sizes ranging from n = 210
to 220 (results are reported for the smallest and largest sample size only) with individual Hurst
parameters h = diag(0.2, 0.4, 0.6, 0.8) and mixing matrix


0.6834 −0.7142 0.6960 −0.1165
 −0.0096 0.4539 −0.0908 0.7740 

(5.1)
P =
 0.4771 −0.2345 0.3359 −0.4243 
0.5525 −0.4784 −0.6281 0.4553
17

(see also Remark 5.2 on the choice of P ). The entrywise Hurst exponents are denoted by hX,i ,
hY,i , i = 1, . . . , n, whereas hX,i
e , i = 1, . . . , n, denotes the Hurst exponents of the demixed sequence
d
−1
−1 .
e = P Y for normalized demixing matrix estimates Pd
X
The results consist of comparisons of the Monte Carlo log-averages of the sample wavelet
fX (2j )ii i, hlog2 W
fY (2j )ii i and hlog2 W
f e (2j )ii i (h·i denotes for Monte Carlo average)
variance hlog2 W
X
for each of the n = 4 components for the sample sizes 220 and 210 (Figures 1 and 4); boxplots for
b
hX,i − hi , b
hY,i − hi and b
hX,i
e − hi , i = 1, 2, 3, 4 (Figures 2 and 5); and boxplots for the 16 entries
d
of P −1 P − I (Figures 3 and 6). Following the procedure described in Remark 3.4, the columns
of Pb were adjusted as to eliminate the non-identifiability factor. In all cases, the sample wavelet
variance matrices were computed based on Daubechies wavelet filters with Nψ = 2 vanishing
moments. Using a different wavelet with Nψ ≥ 2 yields similar conclusions.
fY (2j )ii i
In Figures 1 and 4, as expected for the mixed data Y all components of hlog2 W
fX (2j )ii i. After demixing, all components
display patent departures from the original data hlog2 W
f e (2j )ii i remarkably superimpose those of hlog2 W
fX (2j )ii i, with the possible exception
of hlog2 W
X
of a few coarse scales for h = 0.2 and 0.4. In addition, the boxplots in Figures 2 and 5 show
b
that the Monte Carlo distributions for b
hX,i
e − hi resemble those of hX,i − hi , which illustrates
−1 is very well estimated
the successful demixing of Y . Figures 3 and 6 further indicate that Pd
with negligible biases. In all comparisons, as expected the observed estimator properties improve
significantly when passing from the relatively small sample size 210 to the large sample size 220 ,
hence reflecting the asymptotic statement of Theorem 3.2, (iii). In addition, simulation results
not displayed also show that the standard deviation of the estimates decreases with the sample
√
size according to the scaling ratio C/ ν for some C > 0, as anticipated.
Remark 5.1 Theorem 4.2 leaves open the question of how to optimally choose the octaves J1 <
J2 . For multiple choices of wavelet octaves, namely, J1 = 1 (which involves the largest number
of sum terms in (3.3)) and J2 = 2, . . . , 6, Table 1 shows the performance of the individual Hurst
exponents’ estimators in terms of Monte Carlo bias, standard deviation and (square root) mean
squared error. For sample sizes 220 and 210 , the results indicate that for low values of the Hurst
exponents, the use of two widely separated wavelet octaves produces better results in terms of
mean squared error, whereas for large values of the Hurst exponents the choice of octaves has
little impact on the estimation.

Remark 5.2 Simulation studies not included show that the choice of the mixing matrix (5.1)
does not substantially affect the finite sample results. Moreover, the demixing estimator is very
robust with respect to the condition number of the mixing matrix P . The distributions of the
estimated scalar Hurst eigenvalues after demixing are barely affected for condition numbers of the
order of at least 105 .

5.2

Two-step wavelet-based and maximum likelihood estimation: a comparative study

Due to its wide applicability and well-known asymptotic properties, maximum likelihood estimation is a natural choice and the associated methodology in a multivariate framework has been
constructed by several authors (see references in the Introduction). In this section, we conduct
Monte Carlo experiments to compare the statistical and computational finite sample performances
18

h

J1 , J2

0.20

1,2
1,3
1,4
1,5
1,6
1,2
1,3
1,4
1,5
1,6
1,2
1,3
1,4
1,5
1,6
1,2
1,3
1,4
1,5
1,6

0.40

0.60

0.80

b
h
(220 )
0.25
0.22
0.22
0.21
0.21
0.40
0.40
0.39
0.40
0.39
0.59
0.59
0.59
0.59
0.59
0.79
0.79
0.79
0.79
0.79

bias

sd

0.05
0.02
0.02
0.01
0.01
-0.00
-0.00
-0.01
-0.00
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01
-0.01

0.04
0.03
0.03
0.02
0.02
0.02
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01
0.01

√

MSE

0.06
0.04
0.03
0.03
0.03
0.02
0.02
0.02
0.01
0.01
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02

b
h
(210 )
0.31
0.25
0.24
0.23
0.22
0.45
0.41
0.40
0.40
0.40
0.60
0.58
0.58
0.58
0.58
0.76
0.77
0.77
0.77
0.77

bias

sd

0.11
0.05
0.04
0.03
0.02
0.05
0.01
0.00
0.00
-0.00
-0.00
-0.02
-0.02
-0.02
-0.02
-0.04
-0.03
-0.03
-0.03
-0.03

0.10
0.08
0.08
0.08
0.08
0.08
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07

√

MSE

0.14
0.10
0.09
0.09
0.08
0.10
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.07
0.08
0.07
0.07
0.07
0.07

Table 1:
Choice of scales 1,000 Monte Carlo runs, sample sizes 220 and 210 , h =
(0.2, 0.4, 0.6, 0.8).

H = 0.4
log2 wavelet spectrum

original
mixed
demixed

30
20
10
0

0

5

10
15
log2 scale
H = 0.6

20

40

log2 wavelet spectrum

log2 wavelet spectrum

log2 wavelet spectrum

H = 0.2
40

original
mixed
demixed

30
20
10
0

0

5

10
15
log2 scale

20

60
original
mixed
demixed

40

20

0

0

5

10
15
log2 scale
H = 0.8

20

60
original
mixed
demixed

40

20

0

0

5

10
15
log2 scale

20

Figure 1: Scaling log W·,· (2j ) vs. j for each of the n = 4 components based on the wavelet
variance scales 21 and 22 . The plots were produced by means of 500 Monte Carlo runs of sample
size 220 , with parameter values h = (0.2, 0.4, 0.6, 0.8) and Nψ = 2.

19

H=0.2

0.5
0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

-0.1

-0.1
1

2

H=0.4

0.5

3

1

original mixed demixed

H=0.6

0.5

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

-0.1

-0.1
2

3

H=0.8

0.5

0.4

1

2

original mixed demixed

3

1

original mixed demixed

2

3

original mixed demixed

Figure 2: Boxplots based on the wavelet variance scales 21 and 22 for i = 1, 2, 3, 4, b
hX,i − hi
(hidden, left), b
hY,i − hi (mixed, middle) and b
hX,i
−
h
(demixed,
right),
for
each
of
the
n = 4
i
e
components, sorted by ascending order in terms of h. The plots were produced by means of 500
Monte Carlo runs of sample size 220 , with parameter values h = (0.2, 0.4, 0.6, 0.8) and Nψ = 2.

0.02
0.01
0
−0.01
1

2

3

4

1

2

3

4

1

2

3

4

1

2

3

4

0.01
0
−0.01
−3

x 10
5
0
−5

0.01
0
−0.01

−1 P − I.
Figure 3: Boxplots based on the wavelet variance scales 21 and 22 for the 16 entries of Pd
d
−1
The (i1 , i2 )-th boxplot denotes the (i1 , i2 )-th entry of P P − I. The plots were produced by
means of 500 Monte Carlo runs of sample size 220 , with parameter values h = (0.2, 0.4, 0.6, 0.8)
and Nψ = 2.

20

H = 0.4
log2 wavelet spectrum

original
mixed
demixed

10

5

0

0

2

4
6
log2 scale
H = 0.6

8

15

log2 wavelet spectrum

log2 wavelet spectrum

log2 wavelet spectrum

H = 0.2
15

original
mixed
demixed

10

5

0

0

2

4
6
log2 scale

8

15
original
mixed
demixed

10

5

0

0

2

4
6
log2 scale
H = 0.8

8

20
original
mixed
demixed

15
10
5
0

0

2

4
6
log2 scale

8

Figure 4: log W·,· (2j ) vs. j for each of the n = 4 components based on the wavelet variance scales
21 and 22 . The plots were produced by means of 500 Monte Carlo runs of sample size 210 , with
parameter values h = (0.2, 0.4, 0.6, 0.8) and Nψ = 2.

H=0.2

0.5

0

-0.5

0

1

2

-0.5

3

original mixed demixed

1

2

3

original mixed demixed

H=0.6

0.5

H=0.8

0.5

0

-0.5

H=0.4

0.5

0

1

2

-0.5

3

original mixed demixed

1

2

3

original mixed demixed

Figure 5: Boxplots based on the wavelet variance scales 21 and 22 for i = 1, 2, 3, 4, b
hX,i − hi
b
b
(hidden, left), hY,i − hi (mixed, middle) and hX,i
e − hi (demixed, right), for each of the n = 4
components, sorted by ascending order in terms of h. The plots were produced by means of 500
Monte Carlo runs of sample size 210 , with parameter values h = (0.2, 0.4, 0.6, 0.8) and Nψ = 2.

21

of two-step wavelet-based and maximum likelihood (ML) estimation. For the sake of illustration,
we opt for Whittle-type estimation for fitting a mixed bivariate operator fractional Gaussian
noise. This involves reexpressing the likelihood function in the Fourier domain and using some
approximations. For the reader’s convenience, we provide a brief description of the method; for
more details see, for instance, Hosoya (1996, 1997), Robinson (2008) and Tsai et al. (2017).
In (1.2), suppose X is a vector of two independent fractional Gaussian noise entries with
Hurst parameters hi , i = 1, 2. Then, the (negative) Whittle log-likelihood function of Y can be
approximated by
l(h1 , h2 , A) = 2T log | det A| +

T 
X



e
log |2(1 − cos xi ) det G(xi ; h1 , h2 ) |

i=1



T
X
∗ −1
−1 −1
e
+
tr (A ) {2(1 − cos xi )G(xi ; h1 , h2 )} A IY (xi ) ,

(5.2)

i=1

e h1 , h2 ) :=
where A := P diag(e(h1 ), e(h2 )), e(hi ) := {Γ(2hi + 1) sin(πhi )/2π}1/2 , G(x;
P
1
−2h
−2h
e
e
e
i
i
diag(R(x, h1 ), R(x, h2 )), R(x, hi ) := 4πhi {(2πM − x)
+ (2πM + x)
}+ M
k=−M |x +
−2hi −1 for some large integer M , T := [(ν − 1)/2], I (x) := J (x)J (x)∗ /(2πν), J (x) :=
2kπ|
Y
Y
Y
Y
Pν
Y
exp(itx),
and
x
=
2πi/ν
are
the
Fourier
frequencies.
The
(Whittle)
ML
estimator
is
t
i
t=1
defined by
θb := argminθ l(θ).
(5.3)
In (5.3), l(·) is given by (5.2), and we write θb = (b
h1 , b
h2 , Pb). The estimator (5.3) was implemented
in Matlab using the function fminsearch.m to minimize l(h1 , h2 , A) with respect to the unknown
parameters h1 , h2 and A.
For the simulation study, we picked the parameter values


0.78 0.62
.
(5.4)
(h1 , h2 ) = (0.3, 0.9), P =
0.62 0.78
Monte Carlo averages for the two-step wavelet-based and ML estimators for the parameters h1 ,
h2 and P are reported in Table 2.
The simulation study shows that the semiparametric two-step wavelet-based and the parametric Whittle-type ML methods display comparable
finite sample performances as measured
√
by Monte Carlo bias, standard deviation and MSE. In fact, the former method estimates h1
and P slightly more accurately, whereas the latter does better with h2 . However, the two-step
wavelet-based method is far more computationally efficient. In fact, the ML estimator requires
minimizing (5.2) with respect to n + n2 unknown parameters, which can be numerically very
difficult in higher dimension n. As shown in Table 3, the computational time per realization of
ML grows rapidly as a function of the path size ν, and the ratio between computational times for
the two methods grows exponentially fast. Furthermore, our computational studies indicate that
the minimization procedure required by ML is somewhat sensitive to the initial guess.
In all fairness, the computational performance of ML can be surely improved by replacing
the all-purpose fminsearch.m with a special optimization algorithm. Nevertheless, this computational study illustrates the fact that the potential numerical hurdles in the construction of
viable maximum likelihood estimation for mixed fractional processes are significantly more stringent than those for the proposed two-step wavelet-based method. In addition, the computational
robustness of the latter with respect to the sample path size is striking.

22

1
0
−1
1

2

3

4

1

2

3

4

1

2

3

4

1

2

3

4

1
0
−1

0.5
0
−0.5

1
0
−1

−1 P − I.
Figure 6: Boxplots based on the wavelet variance scales 21 and 22 for the 16 entries of Pd
−1 P − I. The plots were produced by
The (i1 , i2 )-th boxplot denotes the (i1 , i2 )-th entry of Pd
10
means of 500 Monte Carlo runs of sample size 2 , with parameter values h = (0.2, 0.4, 0.6, 0.8)
and Nψ = 2.

method
ML

two-step wavelet method

parameter
h1
h2
p1,1
p1,2
p2,1
p2,2
h1
h2
p1,1
p1,2
p2,1
p2,2

bias
0.1479
-0.0358
-1.1076
4.5516
4.6128
-1.1042
0.0522
-0.1125
-0.0207
0.0182
0.0196
-0.0170

sd
0.1705
0.0761
1.2362
1.2604
1.2244
1.1772
0.0954
0.0919
0.2592
0.3841
0.2462
0.3686

√

MSE
0.2257
0.0841
1.6598
4.7229
4.7725
1.6140
0.1088
0.1452
0.2600
0.3845
0.2469
0.3690

Table 2: Biases, standard deviations and (square root) mean squared errors over 100 replications
with sample size ν = 210 from the two-step wavelet-based and ML methods for the parameters
h1 , h2 and P = (pij )i,j=1,2 as in (5.4).

23

6

Applications

We now provide two applications of the method constructed above.
In Section 6.1, we illustrate the two-step wavelet-based method by fitting a bivariate series of annual tree ring measurements from bristlecone pine trees in California. The data
can be found in the Time Series Data Library, which is available on the website DataMarket (https://datamarket.com/data/list/?q=provider:tsdl). The so-named White Mountain
and Methuselah pine tree data sets are provided by C. W. Ferguson, E. Schulman and H. C. Fritts,
and by D. A. Graybill, respectively. In Section 6.2, we draw upon the results in Section 3.1 to
establish the asymptotic normality of the eigenstructure of the sample wavelet variance matrix
at fixed scales. This is of independent interest because sample wavelet variance matrices do not
generally follow a Wishart distribution. This results from the presence of residual correlation
after the application of the wavelet transform.

6.1

Modeling tree ring data

Many tree ring data sets exhibit long range dependence properties (Tsai and Chan (2005)). Annual
tree ring width measurements can be modeled as aggregates of the underlying continuous time
growth rate process over time intervals between two consecutive sampling time points. Assuming
reasonable physical models, the latter, in turn, can be approximated by a mixed fractional process,
as explained in Section 2.1. Although the full data set covers the period 5142 BC – 1962 AD,
we focus instead on the subperiod 4141 BC – 1962 AD, since preliminary wavelet-based analysis
revealed stationarity in the latter. The time series are displayed in Figure 7, top plots.
Data analysis is conducted both in the time and wavelet domains. We examine the data by
means of sample autocorrelation and cross-correlation functions (ACFs and CCFs, respectively),
f (2j )11 and log2 W
f (2j )22 (see (4.4)) as functions of
main diagonal wavelet scaling plots log2 W
log2 2j = j, as well as the so-named sample wavelet coherence function w
b12 (2j ), j = j1 , . . . , jm .
The latter is a wavelet version of the CCF and can also be used to check the cross-correlation in
bivariate data. For each j, the associated term is defined by
q
j
j
f
f (2j )11 W
f (2j )22
w
b12 (2 ) = W (2 )12
W
(see Whitcher et al. (2000)).
Because it is well known that spurious cross-correlation may occur as a result of the presence
of fractional memory in each time series, it is pivotal to pre-whiten the data (e.g., Cryer and Chan
(2008), Section 11.3). The corresponding sample ACFs, shown on the lower panel in Figure 7,
suggest that the time series have long memory. This is confirmed by wavelet analysis, as displayed

sample path size
28
210
212
214

time in seconds (per realization)
ML
two-step wavelet
2.5
0.0035
22.0
0.0050
216.0
0.0100
2495.0
0.0120

time ratio
(ML/two-step wavelet)
720
4400
21600
213870

Table 3: Computational performance: Whittle-type ML and two-step wavelet based methods,
dimension n = 2.

24

2.0
0.0
1000

3000

5000

0

1000

3000
Time(years)

White Mountains

Methuselah

5000

0.4
0.0

0.4

ACF

0.8

0.8

Time(years)

0.0

ACF

1.0

Methuselah

2.0
1.0
0.0

White Mountains

0

0

50

100

150

200

0

Lag

50

100

150

200

Lag

Figure 7: Upper: Time series plots of tree ring measurements; Lower: Sample autocorrelations of
tree ring measurements.

25

CCF

0.0

0.2

0.4

0.6

0.8

x&y

−200

−100

0

100

200

100

200

Lag

CCF

0.0

0.2

0.4

0.6

0.8

x&y

−200

−100

0
Lag

Figure 8: Upper: Sample cross-correlation between tree ring measurements, after pre-whitening.
Lower: Sample cross-correlation of the demixed data, after pre-whitening. The dashed lines
√
correspond to the threshold ±1.96/ ν at 5% significant level.

26

f (2j )11 and log2 W
f (2j )22 suggest scaling behavior
in Figure 9 (left plot). Indeed, both log2 W
with Hurst parameters that clearly depart from 1/2, i.e., long memory. Moreover, the fact that
both curves resemble each other (namely, close Hurst parameter values) can be explained as the
preponderance of one of the two underlying scaling laws (see the discussion in the Introduction).
The upper panel in Figure 8 displays the sample cross-correlation (for pre-whitened data). It
reveals that the sequences are contemporaneously strongly correlated but not cross-correlated at
any nonzero lag values. This is confirmed by the wavelet coherence function (Figure 9, right plot),
which shows significant and nearly constant correlation across all scales.
The demixing step (S1) of the proposed wavelet-based method yields the following estimated
demixing matrix


0.9112 −0.7827
d
−1
P =
.
0.1467 1.1922
−1 to the original data. Inspection of
Demixed ring tree time series are computed by applying Pd
the sample cross-correlation function for the demixed tree ring data (after pre-whitening) reveals
that the proposed wavelet-based method successfully decorrelated the data (lower panel in Figure
8). This is further confirmed by the wavelet coherence function (Figure 9, right plot), which
evidences near zero correlations at all scales but a few of the coarsest. In addition, both functions
f (2j )11 and log2 W
f (2j )22 (for demixed data) still display scaling behavior. However, the
log2 W
Hurst exponents seem quite distinct and bounded away from 1/2. This is confirmed by the
proposed estimation method. After demixing, the memory parameter estimation step (S2) yields
the parameter estimates b
h1 = 0.65, b
h2 = 0.93 (using scales (j1 , j2 )=(3,7)), and b
h1 = 0.65, b
h2 =
0.96 (using scales (j1 , j2 )=(3,9)) (recall that, in this case, the relation between the Hurst and
memory parameters h and d, respectively, is given by (2.16)). In other words, there is little
sensitivity of the parameter estimates to the choice of octave range. Table 4 further reports a
Monte Carlo study of the sample mean and sample standard deviation of h\
2 − h1 for the case
h1 = h2 = h. The difference between the estimated Hurst parameters for the demixed tree ring
data is b
h2 − b
h1 = 0.96 − 0.65 = 0.31 > 1.645 × sd(h\
2 − h1 ), which lies far outside the confidence
interval. In other words, there is evidence for the hypothesis h1 < h2 in the demixed ring tree
data. Note that this could not have been detected had we skipped step (S1), i.e., if Hurst exponent
estimation had been conducted directly on the original data.

true h
h=0.7

h=0.8

parameter
b
h1
b
h2
\
h2 − h1
b
h1
b
h2
\
h2 − h1

mean
0.6985
0.7229
0.0244
0.7957
0.8229
0.0272

sd
0.0183
0.0176
0.0185
0.0191
0.0195
0.0202

Table 4: wavelet estimation: (j1 , j2 )=(3,9), sample size=6000, number of Monte Carlo
runs=1000.

27

2
1

p
log2 W1,2 (2j )/ (W1,1 (2j )W2,2 (2j ))

1

log2 W (2j ), log2 WXb (2j )

0

-1

-2

-3

-4

W11
W22
WXb,11
WXb,22

-5

-6

-7

-8
0

1

2

3

4

5

6

7

8

9

10

11

j = log2 (scale)

0.8

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

Mixed
Demixed

-0.8

-1
0

1

2

3

4

5

6

7

8

9

10

11

j = log2 (scale)

f (2j )ii (wavelet variances) versus j for bivariate tree ring data.
Figure 9: Left: log2 W
f (2j )22 show scaling
f (2j )11 and log2 W
Before the demixing step (S1) (black), both functions log2 W
behavior with similar Hurst parameter values clearly departing from 1/2. This confirms the presf (2j )22
f (2j )11 and log2 W
ence of long memory. After the demixing step (S1), the functions log2 W
still display scaling behavior, yet with quite distinct Hurst exponents, and clearly departing from
1/2. Right: wavelet coherence function. Before the demixing step (S1) (black), the wavelet
coherence function shows significant (and nearly equivalent) correlations across all scales. After
the demixing step (S1) (red), it shows nearly zero correlation at all scales, which is evidence of
successful demixing.

28

6.2

Asymptotic theory for the eigenstructure of sample wavelet variance matrices

In order to state Theorem 6.1 below, consider the matrix spectral decompositions
b j Lj O
b∗ ,
W (2j ) = O
j

EW (2j ) = Oj Λj Oj∗ ,

bj , Oj ∈ O(n),
O

(6.1)

bj , Oj have columns o
bj,·i , oj,·i , respecwhere Lj := diag(lj,1 , . . . , lj,n ), Λj := diag(λj,1 , . . . , λj,n ), O
tively, for i = 1, . . . , n, and
lj,1 ≤ . . . ≤ lj,n ,

λj,1 ≤ . . . ≤ λj,n ,

bj,1i ≥ 0,
o

oj,1i ≥ 0,

j = j1 , . . . , j m .
(6.2)
In other words, the eigenvalues appearing on the main diagonal entries of Lj and Λj are ordered
bj are all nonnegative, which
from smallest to largest, and the entries on the first row of Oj and O
makes these orthogonal matrices identifiable. Following Magnus and Neudecker (1980), p. 427,
we recall the definition of the so-named duplication matrix D ∈ M (n2 , 21 n(n + 1), R). It consists
of the (unique) operator D that performs the transformation
D(vecS (A))T = (vec(A + A∗ − dg(A)))T ,

i = 1, . . . , n,

A = (ai1 i2 )i1 ,i2 =1,...,n ∈ M (n, R),

(6.3)

where dg(A) := diag(a11 , . . . , ann ). Moreover, for S ∈ S(n, R) with ordered eigenvalues λ1 <
. . . < λn and their respective normalized eigenvectors o·1 , . . . , o·n , we further define the operator


(oT·1 ⊗ oT·1 )D


..


.


T
T


(o·n ⊗ o·n )D


J (S) = 
,
(6.4)
T
+

 (o·1 ⊗ (λ1 In − S) )D 


..


.
+
T
(o·n ⊗ (λn In − S) )D (n+n2 )×n(n+1)/2
where we can apply the relation
vec(A)D = vecS (A + A∗ − dg(A))

(6.5)

(see Lemma 3.7, (i), in Magnus and Neudecker (1980)). The proof of Theorem 6.1 relies on
Proposition B.1, Theorem E.1 (on the weak convergence of eigenvalues and eigenvectors) and the
Delta method.
Theorem 6.1 Let {W (2j )}j=j1 ,...,jm be a set of sample wavelet variance matrices (see (3.3)).
Suppose
EW (2j ) has pairwise distinct eigenvalues, j = j1 , . . . , jm ,
(6.6)
bj , Oj be as in (6.1). Then,
and let F be as in (3.10). Let the matrices Lj , Λj , O
p
T
p
bj − Oj )
Kj vecD (Lj − Λj ), Kj vec(O

d

j=j1 ,...,jm

→ Nn(n+1)m (0, JF J ∗ ),

ν → ∞,

(6.7)

where J = diag(J1 , . . . , Jm ) and Ji , i = 1, . . . , m, is given by J (S) in (6.4) with S := EW (2ji ).
f (2j ).
Remark 6.1 Note that the conclusion in Theorem 6.1 also holds when replacing W (2j ) by W
Remark 6.2 The conclusions of Theorem 6.1 may not hold when the condition (6.6) is not in
place. Proposition F.1 and Example F.1 in Appendix F illustrate this fact in a particular case.
29

A

Asymptotic theory for the wavelet variance of univariate
Gaussian fractional processes

In this section, we establish the asymptotic normality of the wavelet variance of univariate Gaussian fractional processes (n.b.: the framework of Moulines et al. (2007a, 2007b, 2008) is for discrete
time processes). Throughout the section, we assume the underlying wavelet function ψ ∈ L2 (R)
satisfies the conditions (W 1–3), the underlying process {X(t)}t∈R has the form (2.9) or (2.10),
and satisfies assumption (A3). The main result, Theorem A.1, is used in the proof of Proposition
B.2.
The wavelet transform of the univariate process X is defined by
Z
2−j/2 ψ(2−j t − k)X(t)dt, j ∈ N ∪ {0}, k ∈ Z.
d(2j , k) =
R

The wavelet variance at octave j and its natural estimator, the sample wavelet variance, are
denoted by, respectively,
Ew(2j ) := Ed2 (2j , 0),
(A.1)
and

Kj
1 X 2 j
w(2 ) :=
d (2 , k),
Kj
j

Kj =

k=0

ν
,
2j

j = j1 , . . . , jm .

(A.2)

Let ν be the total number of available (wavelet) data points. Throughout this section, we take a
sequence of scaling factor {a(ν)}ν∈N satisfying (3.21).
The following lemma will be used in the subsequent proposition.
Lemma A.1 For any two fixed octaves j, j 0 ∈ N,
Z
j0
j
3−4d
b
b
|x|−4d |ψ(a(ν)2
x)|2 |ψ(a(ν)2
x)|2 ||g ∗ (x)|2 − |g(0)|2 |2 dx = 0,
lim a(ν)
ν→∞

(A.3)

R

where

(
∗

g (x) =

g(x) sin(x/2)
x/2 , d < 1/2;
g(x),
d ≥ 1/2.

Proof: By assumption (A3),
||g ∗ (x)|2 − |g(0)|2 | < C|x|β ,

|x| < δ.

(A.4)

We can break up the integral on the left-hand side of (A.3) into
Z
−4d+3
j
j0
b
b
a
|x|−4d |ψ(a(ν)2
x)|2 |ψ(a(ν)2
x)|2 ||g ∗ (x)|2 − |g(0)|2 |2 dx
|x|<δ

+a−4d+3

Z

0

j
j
b
b
|x|−4d |ψ(a(ν)2
x)|2 |ψ(a(ν)2
x)|2 ||g ∗ (x)|2 − |g(0)|2 |2 dx.

(A.5)

|x|≥δ

We first consider the integration domain |x| < δ. By (A.4) and a change of variable, the first term
in the sum (A.5) is bounded by
Z
−4d+3
j
j0
b
b
Ca(ν)
|x|−4d |ψ(a(ν)2
x)|2 |ψ(a(ν)2
x)|2 |x|2β dx
|x|<δ

30

= Ca(ν)2−2β

Z

0

b j x)|2 |ψ(2
b j x)|2 dx
|x|−4d+2β |ψ(2

|x|<a(ν)δ
2−2β

Z

b j x)|2 |ψ(2
b j 0 x)|2 dx.
|x|−4d+2β |ψ(2

≤ Ca(ν)

R

R
b j x)|2 |ψ(2
b j 0 x)|2 dx < ∞, and
However, (2.14), (2.19) and (2.20) imply that R |x|−4d+2β |ψ(2
a(ν)2−2β → 0 as ν → ∞. So,
Z
j
j0
3−4d
b
b
|x|−4d |ψ(a(ν)2
x)|2 |ψ(a(ν)2
x)|2 ||g ∗ (x)|2 − |g(0)|2 |2 dx → 0,
a(ν)
|x|<δ

as ν → ∞. On the other hand, turning to the integration domain |x| ≥ δ, (2.19) implies that the
second term in the sum (A.5) is bounded by
Z
|x|−4d−4α dx → 0,
Ca(ν)3−4d−4α
|x|≥δ

as ν → ∞. This shows (A.3). 
Proposition A.1 For j, j 0 ∈ N, let w(a(ν)2j ) be as in (A.2). Then,
a(ν)−4d
4d−1 j+j 0

→ 4πb

2

ν
0
Cov(w(a(ν)2j ), w(a(ν)2j ))
a(ν)

4

Z

−2d

|x|

|g(0)|

R

 2j x 
ψb
b

2

 2j 0 x  2
b
dx,
ψ
b

ν → ∞,

(A.6)

0

where b = gcd(2j , 2j ).
Proof: The main argument is similar to the proof of Proposition 3.1 in Wendt et al. (2017), so
we just outline the main steps for the reader’s convenience.
0
It suffices to consider the subsequence ν = a(ν)2j+j ν∗ . By (3.5), the left-hand side of (A.6)
can be reexpressed as
j0

j

2X
ν∗ 2X
ν∗
0
−4d 1
a(ν)
Cov(d2 (a(ν)2j , k), d2 (a(ν)2j , k 0 ))
ν∗
0
k=1 k =1
j0

j

2
2X
ν∗ 2X
ν∗ 
−4d 1
j
j0 0
= 2a(ν)
Ed(a(ν)2 , k)d(a(ν)2 , k )
ν∗
0
k=1 k =1

0
2j ν

= 2a(ν)

−4d+2 j+j 0

2

2j ν

Z
2
0
1 X∗ X∗
ia(ν)(2j k−2j k)x
−2d ∗
2b
j0
j
b
e
|x| |g (x)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
ν∗
R
0
k=1 k =1

j+j 0

= 2a(ν)−4d+2 2

0
2
2j ν 2j ν  Z
0
1 X∗ X∗
ia(ν)(2j k−2j k)x
−2d ∗
2b
j0
j
b
e
|x| |g (x)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
ν∗
R
0

k=1 k =1

Z
2 
0
ia(ν)(2j k−2j k)x
−2d
2b
j0
j
b
−
e
|x| |g(0)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
R

31

j0

j

2
2 ν 2 ν Z
1 X∗ X∗
j
j0
j0
j x)ψ(a(ν)2
b
b
eia(ν)(2 k−2 k)x |x|−2d |g(0)|2 ψ(a(ν)2
+2a(ν)
2
x)dx ,
ν∗
R
k=1 k0 =1
(A.7)
where the first equality is a consequence of the Isserlis theorem. We now show that
−4d+2 j+j 0

j0

j

2
2X
ν∗ 2X
ν∗  Z
0
ia(ν)(2j k−2j k)x
−2d ∗
2b
−4d+2 1
j0
j
b
e
|x| |g (x)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
a(ν)
ν∗
R
0
k=1 k =1

Z
2 
0
ia(ν)(2j k−2j k)x
−2d
2b
j0
j
b
−
e
|x| |g(0)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
→ 0,

ν → ∞.

(A.8)

R

The summation in (A.8) can be reexpressed as (for the details, see the proof of Proposition 3.1,
(iv) in Wendt et al. (2017))
2
X ξr (ν∗ )  Z
j0
j x)ψ(a(ν)2
b
b
a(ν)−4d+2
x)dx
eia(ν)rx |x|−2d |g ∗ (x)|2 ψ(a(ν)2
ν∗
R
r∈Π(ν∗ )

2 
Z
j0
ia(ν)rx
−2d
2b
j
b
=: Θ1 .
−
e
|x| |g(0)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx

(A.9)

R
0

In (A.9), Π(ν∗ ) = gcd(a(ν)2j , a(ν)2j )Z ∩ Bjj 0 (ν∗ ), Bjj 0 (ν∗ ) is the range for r such that the pairs
0
0
(k, k 0 ) satisfying 2j k − 2j k 0 = gcd(2j , 2j )w for some w ∈ Z in the region
0

1 ≤ k 0 ≤ 2j ν∗ ,

1 ≤ k ≤ 2j ν∗ ,
and

ξr (ν∗ )
0
→ gcd(2j , 2j ), ν → ∞.
ν∗
By Parseval’s theorem, the sequences
2 
 Z
j0
ia(ν)rx
−2d ∗
2b
j
b
e
|x| |g (x)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx

(A.10)

r∈Z

R

and
 Z
e

ia(ν)rx

−2d

|x|

j0
j x)ψ(a(ν)2
b
b
x)dx
|g(0)| ψ(a(ν)2
2

2 
r∈Z

R

are summable. Moreover, by (A.10), for large enough ν,
2
X  Z
j j0
ia(ν)rx
−2d ∗
2b
−4d+2
j0
j
b
Θ1 < (gcd(2 , 2 ) + 1)a(ν)
e
|x| |g (x)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
R

r∈Π(ν∗ )

Z
2 
ia(ν)rx
−2d
2b
j0
j
b
−
e
|x| |g(0)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
.
R

X  Z

0

= (gcd(2j , 2j )+1)a(ν)−4d+2

R

r∈Π(ν∗ )

Z
·

e

ia(ν)rx

−2d

|x|

0

j
j x)ψ(a(ν)2
b
b
eia(ν)rx |x|−2d (|g ∗ (x)|2 +|g(0)|2 )ψ(a(ν)2
x)dx

∗

2

2

j0

j x)ψ(a(ν)2
b
b
(|g (x)| − |g(0)| )ψ(a(ν)2
x)dx

R

32





j

j0

−4d+2

≤ (gcd(2 , 2 )+1)a(ν)

 X Z

e

ia(ν)rx

−2d

|x|

j0
j x)ψ(a(ν)2
b
b
(|g (x)| +|g(0)| )ψ(a(ν)2
x)dx
∗

2

2

R

r∈Π(ν∗ )

 X Z
2 1/2
ia(ν)rx
−2d
∗
2
2 b
j0
j
b
·
e
|x| (|g (x)| − |g(0)| )ψ(a(ν)2 x)ψ(a(ν)2 x)dx
,

(A.11)

R

r∈Π(ν∗ )

where the last inequality is a consequence of Cauchy-Schwarz inequality. By Parseval’s theorem,
the first summation term on the right-hand side of (A.11) is bounded by
Z

1/2
j
2 b
j0
2
b
||g (x)| + |g(0)| | |ψ(a(ν)2 x)| |ψ(a(ν)2 x)| dx

−4d

∗

|x|

2

2 2

R
2d−1/2

Z

−4d

|x|

≤ Ca(ν)

b j x)|2 |ψ(2
b j 0 x)|2 dx
|ψ(2

1/2

R

≤ Ca(ν)2d−1/2 .
Turning back to (A.11), this implies that
−2d+3/2

Θ1 ≤ Ca(ν)

 X Z

ia(ν)rx

e

−2d

|x|

2 1/2
j x)ψ(a(ν)2
b
b
x)dx
(|g (x)| − |g(0)| )ψ(a(ν)2
∗

2

j0

2

R

r∈Π(ν∗ )


1/2
Z
−4d+3
−4d
∗
2
2 2 b
j
2 b
j0
2
≤ C a(ν)
|x| ||g (x)| − |g(0)| | |ψ(a(ν)2 x)| |ψ(a(ν)2 x)| dx
→0
R

as ν → ∞. The last inequality is a consequence of Parseval’s theorem, and the limit follows from
Lemma A.1. This proves (A.8), as desired. Consider the last term in the sum (A.7). By an
analogous procedure, we obtain, as ν → ∞,
j0

2a(ν)

−4d+2 j+j 0

2

j

2
2 ν 2 ν Z
0
1 X∗ X∗
ia(ν)(2j k−2j k)x
−2d
2b
j0
j
b
e
|x| |g(0)| ψ(a(ν)2 x)ψ(a(ν)2 x)dx
ν∗
R
0
k=1 k =1
j0

j+j 0 +1

=2

j

2
2 ν 2 ν Z
0
1 X∗ X∗
i(2j k−2j k)x
j0
−2d
2b j
b
e
|x| |g(0)| ψ(2 x)ψ(2 x)dx
ν∗
R
0
k=1 k =1

j

j0

j+j 0

→ 2 gcd(2 , 2 )2

∞ Z
X

= 2b

2

∞ Z
X
z=−∞
0

e

|x|

−2d

b j x)ψ(2
b j 0 x)dx
|g(0)| ψ(2
2

2

R

z=−∞
4d−1 j+j 0

0

i gcd(2j ,2j )zx

e

izx

−2d

|x|

b j x/b)ψ(2
b j 0 x/b)dx
|g(0)| ψ(2
2

2

R

= 4πb4d−1 2j+j |g(0)|4

Z

0

b j x/b)|2 dx,
b j x/b)|2 |ψ(2
|x|−4d |ψ(2

(A.12)

R

where we make a change of variable and the last equality is a consequence of Parseval’s theorem.
By (A.8) and (A.12), (A.6) holds. 

33

2 1/2

Theorem A.1 For a fixed set of octaves 0 < j1 < . . . < jm , let Ew(2j ) and w(2j ) be as in (A.1)
and (A.2), respectively. Then,

 

w(a2j1 )/a2d
Ew(a2j1 )/a2d


p
d




.
.
..
..
ν/a(ν) 
−
 → N (0, W ), ν → ∞,
Ew(a2jm )/a2d
w(a2jm )/a2d
where
Wii0 =

ji +ji0
4πb4d−1
|g(0)|4
ji ji0 2

Z

b j x/bj j 0 )|2 |ψ(2
b i x/bj j 0 ))|2 dx,
x−4d |ψ(2
i i
i i

R

and bji ji0 = gcd(2ji , 2ji0 ), i, i0 = 1, . . . , m.
Proof: The proof can be written as a simple adaptation of the proof of Theorem 3.1. 

B

Proofs and auxiliary results: Section 3

As typical in the asymptotic study of averages, we need investigate the asymptotic covariance of
the sample wavelet transforms W (2j ).
Recall that for a zero mean, Gaussian random vector Z ∈ Rm , the Isserlis theorem (e.g.,
Vignat (2012)) yields
XY
E(Z1 . . . Z2k ) =
E(Zi Zj ), E(Z1 . . . Z2k+1 ) = 0, k = 1, . . . , bm/2c.
(B.1)
PQ
The notation
stands for adding over all possible k-fold products of pairs E(Zi Zj ), where
the indices partition the set 1, . . . , 2k. Proposition B.1 below describes the asymptotic covariance
matrix for the wavelet transform of the mixed fractional process Y at fixed octaves.
Proposition B.1 Suppose Y = {Y (t)}t∈R satisfies the assumptions (A1 – 3). As ν → ∞, for
every pair of octaves j, j 0 ,
(i)
K

Kj
j0
p p
1 1 XX
0
0
ED(2j , k)D(2j , k 0 )∗ ⊗ ED(2j , k)D(2j , k 0 )∗
Kj Kj 0
Kj Kj 0
0
k=1 k =1

→2

(j+j 0 )/2

j

j0

gcd(2 , 2 )

∞
X

Φz gcd(2j ,2j 0 ) ⊗ Φz gcd(2j ,2j 0 ) ,

(B.2)

z=−∞

where

Z
Φz :=

b j x)ψ(2
b j 0 x)e−izx |x|−D G(x)|x|−D∗ dx;
ψ(2

(B.3)

R

(ii) there is a matrix Gjj 0 ∈ M (n(n + 1)/2, R), not necessarily symmetric, such that
p p
0
Kj Kj 0 Cov(vecS W (2j ), vecS W (2j )) → Gjj 0 ,

(B.4)

where the entries of Gjj 0 can be retrieved from (B.2) by means of (B.1) (see (2.3) on the
notation vecS ).

34

Proof: The statement (ii) is a direct consequence of (i), so we only prove the latter. We
proceed as in the proof of Proposition 3.3 (i) in Abry and Didier (2017). p
It suffices
to consider
p
0
0
the subsequence ν = 2j+j ν∗ , ν∗ → ∞. Then, Kj = 2j ν∗ , Kj 0 = 2j ν∗ , and Kj Kj 0 Kj−1 Kj−1
=
0
0 )/2
−(j+j
2
/ν∗ . The covariance between wavelet coefficients can be expressed as
Z Z
0
0
0
0
ED(2j , k)D(2j , k 0 )∗ = 2(j+j )/2 E
ψ(t)ψ(t0 )Y (2j t + 2j k)Y (2j t0 + 2j k 0 )dtdt0
R
0

= 2(j+j )/2

Z Z

Z

ψ(t)ψ(t0 )ei(2

dx
R

R

(j+j 0 )/2

R
j t+2j k)x

∗

|x|−D G(x)|x|−D ei(2j

0 0
t +2j 0 k0 )x

dtdt0

R

Z

b j x)ψ(2
b j 0 x)ei(2j k−2j
ψ(2

=2

0 0
k )x

∗

|x|−D G(x)|x|−D dx,

R

=: Φ2j k−2j 0 k0 .
Let Ξ2j k−2j 0 k0 = Φ2j k−2j 0 k0 ⊗ Φ2j k−2j 0 k0 . By Theorem 1.8 in Jones and Jones (1998), p.10, the
0
0
range of indices spanned by 2j k − 2j k 0 is Zgcd(2j , 2j ). Thus, we would like to show that
∞
X

kΞzgcd(2j ,2j 0 ) k < ∞.

(B.5)

z=−∞

Note that kΞ2j k−2j 0 k0 kl1 = kvec(Φzgcd(2j ,2j 0 ) )vec(Φzgcd(2j ,2j 0 ) )∗ kl1 ≤ kΦzgcd(2j ,2j 0 ) k2l1 . Thus, if
P∞
2
z=−∞ kΦz k < ∞, the expression (B.2) is now a consequence of Lemma E.4 below. In fact,
2

(j+j 0 )/2

kΦz k = 2

Z

b j x)ψ(2
b j 0 x)diag(|x|−2d1 |g ∗ (x)|2 , . . . , |x|−2dn |g ∗ (x)|2 )dxP ∗
eizx ψ(2
n
1

P

2

R

Z

4

≤ CkP k max

1≤i≤n

2

b j x)ψ(2
b j 0 x)|x|−2di |g ∗ (x)|2 dx .
eizx ψ(2
i

R

b j x)ψ(2
b j 0 x)|x|−2di |gi (x)|2 ∈ L2 (R). Thus, by Parseval’s theorem,
For any 1 ≤ i ≤ n, ψ(2
∞
X

Z

2

0

b j x)ψ(2
b j x)|x|−2di |gi (x)|2 dx
eizx ψ(2

R

z=−∞

Z
= 2π

2

b j x)ψ(2
b j 0 x)|x|−2di |gi (x)|2 dx < ∞,
ψ(2

R

this proves

P∞

z=−∞ kΦz

k2

< ∞, as claimed. 

Proof of Theorem 3.1: For notational simplicity, we will restrict ourselves to the bivariate
context (n = 2). The argument for general n can be worked out by a simple adaptation.
The proof is by means of Cramér-Wold device. Form the vector of wavelet coefficients
Vν = (ξ1 (2j1 , 1), ξ2 (2j1 , 1), . . . , ξ1 (2j1 , Kj1 ), ξ2 (2j1 , Kj1 ); . . . ;
ξ1(2jm , 1), ξ2 (2jm , 1), . . . , ξ1 (2jm , Kjm ), ξ2 (2jm , Kjm ))T ∈ RΥ(ν) ,
Pm
where Υ(ν) = 2 jj=j
Kj . Notice that m, j1 , . . . , jm are fixed, but each Kj goes to infinity with
1
ν. Let
α = (αj1 . . . , αjm )T ∈ R3m
35

where
αj = (αj,1 , αj,12 , αj,3 )T ∈ R3 ,

j = j1 , . . . , j m .

Now form the block-diagonal matrix
r
r
r
r


1
1
1
1
1
1
1
1
Dν = diag
Ωj , . . . ,
Ωj ; . . . ;
Ωj , . . . ,
Ωj
,
Kj
2j1 1
Kj1 2j1 1
Kj
2jm m
Kjm 2jm m
| m
| 1
{z
}
{z
}
Kj1

where


Ωj =

Kjm

αj,1
αj,12 /2
αj,12 /2
αj,2


,

j = j1 , . . . , jm .

Let Γ(ν) be the covariance matrix of Vν .
√
d
We would like to show ν(Vν∗ Dν Vν − EVν∗ Dν Vν ) → N (0, σ 2 ) for some σ 2 < ∞. By Lemma
E.1, we only need to prove that
√
(1) σ 2 := limν→∞ Var( νVν∗ Dν Vν ) < ∞;
√
(2) limν→∞ ρ( νDν )ρ(Γ(ν)) = 0,
where ρ(·) is the spectral radius of a matrix.
Statement (1) is a consequence of Proposition B.1, i.e.,

jm X
jm
X
√
αjT
Var( νV ∗ DV ) =
j=j1 j 0 =j1
jm X
jm
X

r

ν
2j

r


ν
j
j0
Cov(vec
W
(2
),
vec
W
(2
))
αj 0 →
S
S
2j 0

αjT Gjj 0 αj 0 < ∞,

ν → ∞.

j=j1 j 0 =j1

To show statement (2), note that, by Lemma E.2,
ρ(Γ(ν)) ≤ ρ(Γ1 ) + . . . + ρ(Γm ),
where Γi is the covariance matrix of Vi := (ξ1 (2ji , 1), ξ2 (2ji , 1), . . . , ξ1 (2ji , Kji ), ξ2 (2ji , Kji ))T , i =
1, . . . , m. Let Ti be the permutation matrix such that
Ti Vi = (ξ1 (2ji , 1), ξ1 (2ji , 2), . . . , ξ1 (2ji , Kji ); ξ2 (2ji , 1), ξ2 (2ji , 2), . . . , ξ2 (2ji , Kji ))T =: Vei ,
and let
e i = EVei Vei∗
Γ

(B.6)

be the covariance matrix of Vei . Then,
e i Ti .
Γi = EVi ViT = E(Ti−1 Vei VeiT Ti ) =: Ti−1 Γ
Since a similarity transformation of a matrix does not change its eigenvalues, we have ρ(Γi ) =
e i ). Let Ys be the s-th entry of Y , and {ξs (2j , k)}k∈Z be the wavelet transform of Ys at octave
ρ(Γ
j and shift k. By Lemma E.2 again, for the matrix in (B.6),
e i ) ≤ ρ(Γi1 ) + ρ(Γi2 ),
ρ(Γ
36

where Γis is the covariance matrix of
Vis := (ξs (2ji , 1), ξs (2ji , 2), . . . , ξs (2ji , Kji ))T ,

i = 1, . . . , m,

s = 1, 2.

On the other hand, note that the covariance between ξs (2j , k) and ξs (2j , k 0 ) is given by
j

0

j

Eξs (2 , k)ξs (2 , k ) =

2
X

p2sl 22jdl

Z
e

i(k−k0 )y

2 −2dl ∗
b
gl
|ψ(y)|
y

R

l=1



y
2j



2

dy.

Thus, {ξs (2j , k)}k∈Z is a stationary sequence for a fixed octave j and its the spectral density can
be expressed as
fj,s (y) =

2
X

p2sl 22jdl

l=1

≤ Cj

∞
X



2
−2dl ∗ y + 2wπ
b
|ψ(y + 2wπ)| |y + 2wπ|
gl
2j
w=−∞

2
∞
X
X

b + 2wπ)|2 |y + 2wπ|−2dl ,
|ψ(y

2

−π < y < π.

(B.7)

l=1 w=−∞

Fix l = 1, 2. The summation in (B.7) is bounded on (−π, π) by using
P (2.20) for w = 0,
and the decay of ψb given by (W 3) for bounding the remaining terms w6=0 . By Lemma E.3
below, ρ(Γis ) < ∞, i = 1, . . . , m, s = 1, 2. Thus, for some C > 0 that does not depend on ν,
√
√
ρ(Γν ) ≤ C < ∞. Since ρ( νDν ) = O( √1ν ), then limν→∞ ρ( νDν )ρ(Γν ) = 0. 
Proof of Theorem 3.2: We first show (i). Let C0 , C1 , R and Λ be as in (3.14) and (3.15). We
now show that, under (3.18), any solution B produced by the EJD algorithm is in the set MEJD .
In view of (3.8), consider the polar decomposition
R = PO,

P is positive definite,

O ∈ O(n).

(B.8)

The decomposition (B.8) always exists for nonsingular, real matrices, and is unique. Thus,
C0 = RR∗ = POO∗ P ∗ = P 2 ,
Since square roots are unique, Step 1 yields
W = P −1 .

(B.9)

W C1 W ∗ = W (POΛO∗ P ∗ )W ∗ = OΛO∗ .

(B.10)

Step 2 and (3.14) imply that

By (2.11), we can assume that the eigenvalues of Λ (see (3.15)) are ordered from smallest to
largest, in which case the column vector o·i in O is associated with the eigenvalue θi , where


 2 Z
2
Z
y
y
2di (J2 −J1 )
2
−2d
∗
2
−2d
∗
b
b
θi = 2
|ψ(y)|
|y| i gi
dy
|ψ(y)|
|y| i gi
dy,
(B.11)
2J2
2J1
R
R
i = 1, . . . , n. However, in the spectral decomposition in Step 2, each orthogonal eigenvector is
determined up to multiplication by −1. Thus, for I as in (3.16), Q∗ ∈ OI, and any demixing
matrix B produced by the EJD algorithm has the form
B = QW ∈ I(PO)−1 = IR−1 = I(P E(2J1 )
37

1/2

diag(2J1 d1 , . . . , 2J1 dn ))−1 .

(B.12)

In other words, B ∈ MEJD . Conversely, it is clear that any matrix in MEJD can be attained as
a solution to the EJD algorithm under (3.18). This establishes (i).
b0 = W (2J1 ), C
b1 = W (2J2 ) (we
To show (ii), consider the EJD algorithm with input matrices C
J
bk to avoid confusion with their deterministic counterparts Ck = EW (2 k ), k = 1, 2). By
write C
replacing all matrices in the proof of (i) with their sample counterparts and following the same
argument, the set of solutions to the EJD algorithm is made up of matrices of the form
b∗ C
b −1/2 ,
bν = Πν O
B
0

where Πν ∈ I,

b −1/2 C
b1 C
b −1/2 = O
bΛ
bO
b∗ ,
C
0
0

P
b and diagonal Λ.
b Note that C
b0 →
C0 , by
for some spectral decomposition with orthogonal O
Theorem 3.1. Since the square root is unique and C0 is invertible, then Theorem E.1 implies that,
b −1/2 exists. Thus, by Theorem 3.1, and
with probability going to 1, the inverse square root C
0
P
−1 C P −1 . However, P −1 C P −1 is a symmetric positive
b −1/2 C
b1 C
b −1/2 →
Slutsky’s theorem, C
P
1
1
0
0
definite matrix that admits the spectral decomposition OΛO∗ with pairwise distinct eigenvalues
b −1/2 with probability going to 1.
b −1/2 C
b1 C
(see (B.9) and (B.10)). Then, by Theorem E.1, so is C
0
0
b −1/2 C
b1 C
b −1/2 whose
Therefore, Theorem E.1 implies that there is a spectral decomposition of C
0
0
P
P
b and Λ,
b respectively, satisfy O
b →
b →
bν =
eigenvector and eigenvalue matrices O
O, Λ
Λ. So, B
−1/2 P
−1/2
∗
∗
−1
b
b
b
Πν O C 0
→ ΠO C0
= Π(PO) for some Π ∈ I, i.e., the sequence Bν satisfies (3.19).
√
d
b0 − C0 ), vecS (C
b1 − C1 ))T →
We now show (iii). From Theorem 3.1, ν(vecS (C
N (0, F ), where
F ∈ S+ (n(n + 1), R). Therefore, we can write

b1 = C1 + √1 Z2,ν .
C
ν

(B.13)


T
√
d
ν vecS (Z1,ν ), vecS (Z2,ν )
→ N (0, F ).

(B.14)

b0 = C0 + √1 Z1,ν ,
C
ν
for two random matrices Z1,ν and Z2,ν such that

Since

c=C
b −1/2 = (C0 + √1 Z1,ν )−1/2 = C −1/2 (In + √1 C −1 Z1,ν )−1/2
W
0
0
ν
ν 0

 
 
1 1 −3/2
1
1 1
1
−1/2
= C0
In − √ C0−1 Z1,ν + OP
= W − √ C0 Z1,ν + OP
,
2 ν
ν
2 ν
ν

(B.15)

where the fourth equality is the Taylor expansion of a matrix function (namely, the function
(1 + x)−1/2 , where we replace 1 and x with I and a matrix A, respectively; see Golub and Van
Loan (2012), p. 565). Then, we arrive at


√
c − W ) = − 1 C −3/2 Z1,ν + OP √1 .
ν(W
(B.16)
2 0
ν
On the other hand, by (B.13) and (B.15),


 
 
1 1 −3/2
1
1
1
−3/2 ∗
∗
∗ 1 1
∗
c
b
c
W C1 W = W − √ C0 Z1,ν +OP
C1 + √ Z2,ν
W − √ Z1,ν (C0 ) +OP
2 ν
ν
2 ν
ν
ν


 
1
1
1
1 −3/2
−3/2 ∗
∗
∗
∗
∗
+ OP
= W C1 W + √ W Z2,ν W − W C1 Z1,ν (C0 ) − C0 Z1,ν C1 W
,
2
2
ν
ν
38

thus,
√



1
1
1 −3/2
−3/2 ∗
∗
∗
∗
∗
∗
c
b
c
ν(W C1 W − W C1 W ) = W Z2,ν W − W C1 Z1,ν (C0 ) − C0 Z1,ν C1 W + OP √ .
2
2
ν
(B.17)
As a consequence, there are matrices

n(n + 1) 
A1 = A1 (C0 ) ∈ M n2 ,
,R
2

and
A2 = A2 (C0 , C1 ) ∈ M
such that




1 −3/2
vec − C0 Z1,ν
2

 n(n + 1)
2

T

, n + n2 , R



= A1 (vecS (Z1,ν ))T ,

(B.18)



T
1
1 −3/2
−3/2
∗
vecS W Z2,ν W ∗ − W C1 Z1,ν
(C0 )∗ − C0 Z1,ν C1 W ∗
2
2
= A2 (vecS (Z1,ν ), vecS (Z2,ν ))T .

(B.19)

By (B.14)–(B.19),
√

d

c − W ), vecS (W
cC
b1 W
c ∗ − W C1 W ∗ ))T → N (0, Σ1 ),
ν(vec(W

where

Σ1 =

e1
A
A2




Σ

e1
A
A2

∗

∈ M (n2 + n(n + 1)/2, R),

(B.20)

(B.21)

e1 = (A1 , 0n2 ×n(n+1)/2 ). In Step 2 of the EJD algorithm, write out the spectral decompoand A
cC
b1 W
c∗ = Q
b∗ D
b 1 Q.
b Recall that we
sition W C1 W ∗ = Q∗ D1 Q and also its estimated counterpart W
need to show the asymptotic normality of the random vector
bW
c ).
vec(Q

(B.22)

From the ordering of eigenvalues in (B.10) and expression (3.15), W C1 W ∗ has pairwise distinct
eigenvalues θ1 < . . . < θn , where θi is defined by (B.11). So, by the Delta method,
√

d

∗
c − W ), vec(Q
b ∗ − Q∗ ))T → N (0, JQ Σ1 JQ
ν(vec(W
).

(B.23)

In (B.23),
JQ = diag(In2 , Jq ) ∈ M (2n2 , n2 + n(n + 1)/2),

(B.24)

and Jq is given by
 
q1· ⊗ (θ1 In − W C1 W ∗ )+ D


..
2
Jq = 
 ∈ M (n , n(n + 1)/2, R)
.

qn· ⊗ (θn In − W C1 W ∗ )+ D


(cf. expression (6.4)), where the vector qi· denotes the i-th row of Q ∈ O(n). In view of
b instead of Q
b ∗ . So, let
(B.22), we need to establish the asymptotic behavior of the matrix Q,
39

T = (ti1 i2 )i1 ,i2 =1,...,n2 be the permutation operator defined by the transformation T (vec(R∗ ))T =
(vec(R))T , R ∈ M (n, R), i.e.,

1, i1 = (k − 1)n + p, i2 = (p − 1)n + k, k, p = 1, . . . , n;
t i1 i2 =
0, otherwise.
Thus,

√

d
c − W ), vec(Q
b − Q))T →
ν(vec(W
N (0, Σ2 ),

where
∗
Σ2 = diag(In2 , T )JQ Σ2 JQ
diag(In2 , T )∗ .

(B.25)

We arrive at the relations
c = W + √1 Z3,ν ,
W
ν

b = Q + √1 Z4,ν ,
Q
ν

d

where (vec(Z3,ν ), vec(Z4,ν ))T → N (0, Σ2 ). Therefore,
√



bW
c − QW ) = QZ3,ν + Z4,ν W + OP √1 .
ν(Q
ν

Therefore, for some matrix
A3 = A3 (O0 , Λ0 , C1 ) ∈ M (n2 , 2n2 ),
we can write
(vec(QZ3,ν + Z4,ν W ))T = A3 (vec(Z3,ν ), vec(Z4,ν ))T .
Hence,

√

bν ) − vec(B))T =
ν(vec(B

√

(B.26)

d

bW
c ) − vec(QW ))T → N (0, A3 Σ2 A∗ ),
ν(vec(Q
3

as claimed. 
The next proposition gives the asymptotic distribution of the main diagonal entries of the
b In its proof, we make use of the following
sample wavelet variance of the demixed process X.
lemma.
Lemma B.1 For a fixed Π ∈ I, let
bν (Π diag(2−J1 d1 , . . . , 2−J1 dn )E(2J1 )−1/2 P −1 )−1 ,
Ibν = B

(B.27)

bν is post-multiplied by the inverse of the limiting matrix on the right-hand side of (3.19).
i.e., B
Then,
√
d
ν(vec(Ibν − I))T → N (0, Σ(J1 , J2 )), ν → ∞,
(B.28)
for some positive semidefinite matrix Σ(J1 , J2 ).

40

Proof: There exists a matrix TP ∈ M (n2 , R) such that
bν − Π diag(2−J1 d1 , . . . , 2−J1 dn )E(2J1 )−1/2 P −1 ))T ,
(vec(Ibν − I))T = TP (vec(B
Then, by (3.20) and the Delta method, the limit in distribution (B.28) holds for Σ(J1 , J2 ) =
TP ΣF (J1 , J2 )TP∗ . 
bν be the demixing matrix described in (3.19). For Ibν as in (B.27), let
So, let B
D := Π diag(2−J1 d1 , . . . , 2−J1 dn )E(2J1 )−1/2 ,

(B.29)

b (see (3.22)) can be reexpressed as
which is a diagonal matrix. Then, the demixed process X
b := B
bν P D−1 DX = Ibν DX.
X
b be the demixed process defined by (3.22), let
Proposition B.2 For j = j1 , . . . , jm , let X
j
b and let EWX (a(ν)2j ) be the wavelet variance
WXb (a(ν)2 ) be the sample wavelet variance of X,
of the hidden process X. Then,

p

−2d1

ν/a(ν)diag(a(ν)

−2dn

, . . . , a(ν)

 T
) vecD (WXb (a(ν)2 ) − DEWX (a(ν)2 )D)


j

j

j=j1 ,...,jm

d

→ N (0, KWK∗ ),

(B.30)

as ν → ∞ (see (2.3) on the notation vecD ). In (B.30),
K = diag(D2 , . . . , D2 ),
| {z }

(B.31)

m

D is given by (B.29), and

W(k1 , k2 ) =

wl,v,i , k1 = (l − 1)n + i, k2 = (v − 1)n + i;
0,
otherwise,

where
4
i −1
wl,v,i = 4πb4d
jl jv |gi (0)|

Z

(B.32)

b jl x/bj j )|2 |ψ(2
b jv x/bj j ))|2 dx,
|x|−4di |ψ(2
l v
l v

R

and bjl jv =

gcd(2jl , 2jv ),

for l, v = 1, . . . , m, i = 1, . . . , n.

bν yields a
Remark B.1 Intuitively, Proposition B.2 says that the demixing matrix estimator B
b
demixed process X that is close to the hidden X up to a non-identifiability factor D. In fact, the
limiting distribution of


p
−2d1
−2dn
j
j
T
ν/a(ν)diag(a(ν)
, . . . , a(ν)
)(vecD (WDX (a(ν)2 ) − EWDX (a(ν)2 )))
j=j1 ,...,jm

is also N (0, KWK∗ ). In particular, the main diagonal entries of the sample wavelet variance
b are asymptotically independent.
WXb (a(ν)2j ) of the demixed process X

41

b = Ibν DX, then,
Proof of Proposition B.2: Since X
WXb (a(ν)2j ) = (Ibν )DWX (a(ν)2j )D(Ibν )∗ .
Thus,
WXb (a(ν)2j ) − DEWX (a(ν)2j )D = (Ibν )DWX (a(ν)2j )D(Ibν )∗ − DEWX (a(ν)2j ))D


j
−1
j
−1 ∗
b
b
b
= (Iν ) DWX (a(ν)2 )D − (Iν ) DEWX (a(ν)2 )D((Iν ) ) (Ibν )∗

 

j
j
−1
j
b
b
= (Iν ) D(WX (a(ν)2 ) − EWX (a(ν)2 ))D − ((Iν ) − I)DEWX (a(ν)2 )D



− DEWX (a(ν)2j )D ((Ibν )−1 )∗ − I



−1
j
−1 ∗
b
b
− ((Iν ) − I)DEWX (a(ν)2 )D ((Iν ) ) − I
(Ibν )∗ .

(B.33)

Recall that the operator vecD (WX (a(ν)2j )) picks out the main diagonal entries of the matrix
WX (a(ν)2j ), which are independent. Therefore, by Theorem A.1 for univariate processes,


p
−2d1
−2dn
j
j
T
ν/a(ν)diag(a(ν)
, . . . , a(ν)
)(vecD (D(WX (a(ν)2 ) − EWX (a(ν)2 ))D))
j=j1 ,...,jm



p
j
j
T
−2d1
−2dn
)(vecD (WX (a(ν)2 ) − EWX (a(ν)2 )))
=K
ν/a(ν)diag(a(ν)
, . . . , a(ν)
j=j1 ,...,jm
d

→ N (0, KWK∗ ),

ν → ∞.

(B.34)

In (B.34), the matrices W and K are defined by (B.32) and (B.31), respectively. By (B.28) and
the Delta method,
√
d
( νvec((Ibν )−1 − I))T → N (0, Σ(J1 , J2 )).
Since EWX (a(ν)2j ) = diag(EWX (a(ν)2j )11 , . . . , EWX (a(ν)2j )nn ), then
(vecD (((Ibν )−1 − I)DEWX (a(ν)2j ))D)T =
D2 diag(EWX (a(ν)2j )11 , . . . , EWX (a(ν)2j )nn )(vecD ((Ibν )−1 − I))T .
Therefore,
p
ν/a(ν)diag(a(ν)−2d1 , . . . , a(ν)−2dn )(vecD (((Ibν )−1 − I)DEWX (a(ν)2j )D))T


2
−2d1
−2dn
j
j
= D diag(a(ν)
, . . . , a(ν)
)diag(EWX (a(ν)2 )11 , . . . , EWX (a(ν)2 )nn )




√
1
1
−1
T
b
= OP p
· ( νvecD ((Iν ) − I))
·p
.
a(ν)
a(ν)
Similarly,
p
ν/a(ν)diag(a(ν)−2d1 , . . . , a(ν)−2dn )(vecD (DEWX (a(ν)2j )D((Ibν )−1 − I)))T
42

(B.35)


= OP


1
p
,
a(ν)

(B.36)

and
p
ν/a(ν)diag(a(ν)−2d1 , . . . , a(ν)−2dn )


1
−1
j
−1 ∗
T
b
b
·(vecD ((Iν ) − I)DEWX (a(ν)2 )D(((Iν ) ) − I)) = OP √ .
ν

(B.37)

Consequently, by (B.33)-(B.37) and Slutsky’s theorem, the limiting distribution of
r

ν
−2d1
−2dn
j
j
T
diag(a(ν)
, . . . , a(ν)
)(vecD (WXb (a(ν)2 ) − DEWX (a(ν)2 )D))
a(ν)
j=j1 ,...,jm
is equal to the limiting distribution of

r
ν
−2d1
−2dn
j
j
T
K
diag(a(ν)
, . . . , a(ν)
)(vecD (WX (a(ν)2 ) − EWX (a(ν)2 )))
,
a(ν)
j=j1 ,...,jm
as claimed. 
The next proposition provides a bound on the difference between the wavelet variance of the
entrywise process Xi and the scaling factor 2j2di |gi (0)|2 K(di ), i = 1, . . . , n. This bound is useful
because of the general absence of exact self-similarity in (2.9), and it is applied in the proof of
Theorem 3.3.
Proposition B.3 For i = 1, . . . , n, let EWX (·)ii0 be defined by (3.24). Then,
|EWX (2j )ii − 2j2di |gi (0)|2 K(di )| ≤ C2j(2di −β) ,
where K(h) =

R

R |ψ(x)|

b

(B.38)

2 |x|−2d dx.

Proof: In fact, for i = 1, . . . , n,
j

|EWX (2 )ii −

2j2di gi2 (0)K(d)|

j

Z

b j x)|2 (||g ∗ (x)|2 − |gi (0)|2 |)|x|−2di dx
|ψ(2
i

=2

R

≤2

j

Z
|x|<δ

b j x)|2 ||g ∗ (x)|2 − |gi (0)|2 ||x|−2di dx + 2j
|ψ(2
i

Z
|x|≥δ

b j x)|2 ||g ∗ (x)|2 − |gi (0)|2 ||x|−2di dx.
|ψ(2
i
(B.39)

By (2.13), the first sum term on the right-hand side of (B.39) is bounded by
Z
Z
j
j
2
β
−2di
j(2di −β)
2
b
b
C2
|ψ(2 x)| |x| |x|
dx = C2
|ψ(x)|
|x|−2di +β dx
|x|<2j δ

|x|<δ

≤ C2

j(2di −β)

Z

2
b
|ψ(x)|
|x|−2di +β dx ≤ C2j(2di −β)

R

Z

2
b
|ψ(x)|
|x|−2di +β dx.

(B.40)

R

By (2.20), the integrand in (B.40) behaves like |x|2Nψ −2di +β around the origin. Also, by (2.19),
the integrand is bounded
by |x|β−2α−2di as |x| → ∞, where β − 2α − 2di < −1 as a consequence
R
2
b
of (2.14). Thus, R |ψ(x)|
|x|−2di +β dx < ∞ and
Z
j
b j x)|2 ||g ∗ (x)|2 − |gi (0)|2 ||x|−2di dx ≤ C2j(2di −β) .
2
|ψ(2
i
|x|<δ

43

Moreover, by (2.19) and the fact that gi∗ (x) is bounded, the second sum term on the right-hand
side of (B.39) is bounded by
Z
−2jα j
|x|−(2di +2α) dx ≤ C2j(2di −β) .
C2
2
|x|≥δ

The last inequality holds because

R
|x|>π

|x|−(2di +2α) dx < ∞ and −2α < 2di − β − 1. Consequently,

|EWX (2j )ii − 2j2di |gi (0)|2 K(d)| < C2j(2di −β) ,
as claimed. 
The proof of Theorem 3.3, presented next, is similar to that of Proposition 3 in Moulines et
al. (2007a).
Proof of Theorem 3.3: Recast (B.30) as

 
a(ν)−2d1 WXb (a(ν)2j1 )11
a(ν)−2d1 EWX (a(ν)2j1 )11 D(1, 1)2

 
..
..

 
.
.

 
−2d
j
−2d
m
1
1



WXb (a(ν)2 )11   a(ν)
EWX (a(ν)2jm )11 D(1, 1)2
  a(ν)
r
ν



..
..

−
.
.


a(ν) 
 a(ν)−2dn W b (a(ν)2j1 )nn   a(ν)−2dn EWX (a(ν)2j1 )nn D(n, n)2

 
X

 
..
..

 
.
.
a(ν)−2dn WXb (a(ν)2jm )nn
a(ν)−2dn EWX (a(ν)2jm )nn D(n, n)2







 d
 → N (0, G),






(B.41)
where WXb (·)ii and EWX (·)ii are defined by (3.24). The limiting covariance matrix is block
diagonal and can be written as G = diag(G1 , . . . , Gn ). For i = 1, . . . , n, Gi is a m × m matrix
whose (k1 , k2 )-th entry is given by
Z
4di −1 jk1 +jk2
4
2
b jk2 x/bj ,j ))|2 dx,
b jk1 x/bj ,j )|2 |ψ(2
Gi (k1 , k2 ) = 4πbjk ,jk 2
|gi (0)| D(i, i)
x−4di |ψ(2
k1 k2
k1 k2
1

2

R

where bjk1 ,jk2 = gcd(2jk1 , 2jk2 ) for i = 1 . . . , n, k1 , k2 = 1 . . . , m. However, under condition (3.21),
relation (B.38) implies that
p
ν/aa−2di |EWX (a(ν)2j )ii − |gi (0)|2 K(di )(a2j )2di |
p
p
≤ C ν/aa−2hi a2di −β 2j(2di −β) ≤ C ν/aa−β → 0, ν → ∞,
(B.42)
for i = 1, . . . , n. As a consequence of (B.41) and (B.42),

 
|g1 (0)|2 K(d1 )22j1 d1 D(1, 1)2
a(ν)−2d1 WXb (a(ν)2j111 )



..
.
..

 
.

 
2
−2d
j



m
1
WXb (a(ν)2 )11   |g1 (0)| K(d1 )22jm d1 D(1, 1)2
  a(ν)
r
ν

 
..
..

−
.
.
 
a(ν) 
 a(ν)−2dn W b (a(ν)2j1 )nn   |gn (0)|2 K(dn )22j1 dn D(n, n)2
 

X

 
..
..
 

.
.
|gn (0)|2 K(dn )22jm dn D(n, n)2
a(ν)−2dn WXb (a(ν)2jm )nn

44







 d
 → N (0, G).





(B.43)

Define

X
m

f (x) =

wl1 log(x1l ), . . . ,

l=1

m
X

wln log(xnl )

T
,

l=1

i
for x = (x11 , . . . , x1m ; . . . ; xn1 , . . . , xnm )T ∈ Rnm
+ and w as in (3.25), i = 1, . . . , n. Let yν and y0
be the left and right vectors in the difference between parentheses on the left-hand side of (B.43).
Then, f (yν ) = (db1 , . . . , dbn ) and f (y0 ) = (d1 , . . . , dn ). By (B.43) and the Delta method,



p

ν/a(ν) 

 
db1
..  − 
.  
dbn


d1

d
..  →
N (0, ∇f (y0 )G∇f (y0 )T ),
. 
dn

where
∇f (y0 ) = diag(A1 , . . . , An ),
and

Ai =


i
w1i
wm
,...,
,
|gi (0)|2 K(di )22j1 di D(i, i)2
|gi (0)|2 K(di )22jm di D(i, i)2

i = 1, . . . , n.

This establishes (3.28). 
Proof of Corollary 3.1:
Note that Proposition 3.1 and Theorem 3.1 also hold under
0
assumptions (A1 ), (A2) and (A30 ). In addition, condition (3.12) follows from (A30 ), so, by the
same arguments for the proofs of Theorems 3.1, 3.1 and 3.3, the claim holds. 

C

Proofs and auxiliary results: Section 4

As a consequence of applying (4.2) and doing a direct computation, the integral representation
of the wavelet covariance in discrete time is provided in the following proposition.
Proposition C.1 Let {Y (k)}k∈Z be the sequence (4.1). For all j, j 0 ≥ 0 and k, k 0 ∈ Z,
Z
j
j0 0
∗
e j , k), D(2
e j 0 , k 0 )) =
Cov(D(2
Hj (x)Hj 0 (x)eix(2 k−2 k ) |x|−D G(x)|x|−D dx,
R

where
−j/2

Hj (x) = 2

Z X

ψ(2−j s)ϕ(s + l)e−ixl ds.

(C.1)

R l∈Z

P
−j/2
e j
Proof: Let Yet = ∞
l=−∞ Yl ϕ(t − l). Then, D(2 , k) = 2

R

−j
e
R Yt ψ(2 t

− k)dt. Therefore,

e j , k), D(2
e j 0 , k 0 ))
Cov(D(2
−j

=2

Z Z
E
R

−j

=2

Z Z
E
R

∞
X

∞
X

0

ψ(2−j t − k)ψ(2−j t0 − k 0 )ϕ(t − l)ϕ(t0 − l0 )Yl Yl∗0 dtdt0

R l=−∞ l0 =−∞
∞
X

∞
X

0

ψ(2−j t)ψ(2−j t0 )ϕ(t + l)ϕ(t0 + l0 )Y2j k−l Y2∗j 0 k0 −l0 dtdt0 .

R l=−∞ l0 =−∞

45

(C.2)

By (2.9), (2.10) and (4.2), we can reexpress (C.2) as
2−j

Z Z Z
R

R

∞
X

∞
X

0

ψ(2−j t)ψ(2−j t0 )ϕ(t + l)ϕ(t0 + l0 )

R l=−∞ l0 =−∞
j0 0

j

0

∗

ei(2 k−l)x e−i(2 k −l )x |x|−D G(x)|x|−D dtdt0 dx,
Z
j
j0 0
∗
=
Hj (x)Hj 0 (x)eix(2 k−2 k ) |x|−D G(x)|x|−D dx,

(C.3)

R

where G(x) and Hj (x) are defined by (3.6) and (C.1), respectively. Note that, by Proposition 3
in Moulines et al. (2007b),
|Hj (x)| = O(|x|Nψ ), x → 0,
(C.4)
and
|Hj (x)| ≤ C,

x ∈ R,

(C.5)

so the integral on the right-hand side of (C.3) is finite. 
The next result is the discrete time analogue of Proposition B.1.
Proposition C.2 Let {Y (k)}k∈Z be the sequence (4.1). For every pair of octaves j, j 0 ≥ 0,
(i)
K

p

Kj
j0
p
1 1 XX e j
e j 0 , k 0 )∗ ⊗ ED(2
e j , k)D(2
e j 0 , k 0 )∗
ED(2 , k)D(2
Kj Kj 0
Kj Kj 0
0
k=1 k =1

→2

(j+j 0 )/2

j

j0

gcd(2 , 2 )

∞
X

e
e
Φ
z gcd(2j ,2j 0 ) ⊗ Φz gcd(2j ,2j 0 ) ,

ν → ∞,

(C.6)

z=−∞

where

Z
e z :=
Φ

∗

Hj 0 (x)Hj (x)e−izx |x|−D G(x)|x|−D dx;

(C.7)

R

e jj 0 ∈ M (n(n + 1)/2, R), not necessarily symmetric, such that
(ii) there is a matrix G
p p
f (2j ), vecS W
f (2j 0 )) → G
e jj 0 ,
Kj Kj 0 Cov(vecS W

ν → ∞,

(C.8)

e jj 0 can be retrieved from (C.6) by means of (B.1) (see (2.3) on the
where the entries of G
notation vecS ).
Proof: Following the same argument as in the proof of Proposition B.1, we only need to show
e z k2 is summable, where
that kΦ
Z
∗
e
Φz :=
Hj (x)Hj 0 (x)eixz |x|−D G(x)|x|−D dx.
R

Since
Z
ezk = P
kΦ

2

eizx Hj 0 (x)Hj (x)diag(|x|−2d1 |g1∗ (x)|2 , . . . , |x|−2dn |gn∗ (x)|2 )dxP ∗

R

46

Z

4

1≤i≤n

2

0

eizx Hj 0 (x)Hj (x)(2j x)|x|−2di dx .

≤ CkP k max

R

Moreover, for any 1 ≤ i ≤ n, by (C.4) and (C.5), Hj 0 (x)Hj (x)x−2di ∈ L2 (R). Thus, by Parseval’s
theorem,
∞
X
z=−∞

Z

2

eizx Hj 0 (x)Hj (x)|x|−2di dx

2

Z

Hj 0 (x)Hj (x)|x|−2di dx < ∞.

=
R

R

e z k2 is summable. Therefore, (C.6) and (C.8) hold. 
Hence, kΦ
eν P Λ
e 1/2 Π and
Define the matrices Ieν = B
J1
e = ΠΛ
e −1/2 .
D
J1

(C.9)

Then, we can reexpress the demixed process (C.10) as
e
e := Ieν DX.
X

(C.10)

The following proposition gives the asymptotic distribution of the main diagonal entries of the
f e of the demixed process X.
e Note that there is a distinction between
sample wavelet variance W
X
f e and W
fX in the proof: the latter denotes the sample wavelet variance of X.
W
X
e be the demixed process (C.10), let W
f e (a(ν)2j ) be
Proposition C.3 For j = j1 , . . . , jm , let X
X
e EW
fX (a(ν)2j ) be the wavelet variance of the hidden process
the sample wavelet variance for X,
X. Then,
r

ν
e W
e T
f e (a(ν)2j ) − DE
fX (a(ν)2j )D))
diag(a(ν)−2d1 , . . . , a(ν)−2dn )(vecD (W
X
a(ν)
j=j1 ,...,jm
d
fK
eW
e ∗ ),
→ N (0, K

(C.11)

e is given by
e 2, . . . , D
e 2 ) and D
e = diag(D
as ν → ∞ (see (2.3) on the notation vecD ). In (C.11), K
| {z }
m

(C.9). The (k1 , k2 )-th entry of the limiting covariance matrix is given by

w
el,v,i , k1 = (l − 1)n + i, k2 = (v − 1)n + i;
f
W(k1 , k2 ) =
0,
otherwise,
Rπ
where w
el,v,i = 4π|gi (0)|4 24di max(jl ,jv )+min(jl ,jv ) −π |D|jl −jv | (x; di )|2 dx, for l, v = 1, . . . , m, i =
1, . . . , n,
X
b + 2kπ)ψ(2
b −u (x + 2kπ))
(C.12)
Du (x, d) =
|x + 2kπ|−2d eu (x + 2kπ)ψ(x
k∈Z

and, for all u ≥ 0,
−u x

eu (x) = 2−u/2 (1, ei2

u −1)2−u x

, . . . , e−i(2

47

)T ,

x ∈ R.

b with X.
e Then, the limiting
Proof: In the argument for proving Proposition B.2, replace X
distribution of
r

ν
−2h1
−2hn
j
j e T
e
f
f
diag(a(ν)
, . . . , a(ν)
)(vecD (WXe (a(ν)2 ) − DEWX (a(ν)2 )D))
a(ν)
j=j1 ,...,jm
is equal to the limiting distribution of
r

ν
−2h1
−2hn
j
j
e
f
f
K
diag(a(ν)
, . . . , a(ν)
)vecD (WX (a(ν)2 ) − EWX (a(ν)2 ))
,
a(ν)
j=j1 ,...,jm
which only involves main diagonal entries. So, fix i = 1, . . . , n. By (2.9), the generalized spectral
density (Yaglom (1958)) of the i-th component of X is
∞
X

fi (x) = |eix − 1|2

|x + 2lπ|−2di −2 |gi (x + 2lπ)|2 ,

x ∈ [−π, π)

l=−∞

for −1/2 ≤ di < 1/2, and
fi (x) =

∞
X

|x + 2lπ|−2di |gi (x + 2lπ)|2 ,

x ∈ [−π, π)

l=−∞

for di ≥ 1/2. Reexpress fi as

fi (x) = |1 − e−ix |−2di fi∗ (x),

where
fi∗ (x)

2 sin(x/2)
=
x

2di +2

X

|gi (x)|2 + |2 sin(x/2)|2di +2

|x + 2lπ|−2di −2 |gi (x + 2lπ)|2 ,

(C.13)

l6=0

for −1/2 ≤ di < 1/2, and
fi∗ (x)

2di

2 sin(x/2)
=
x

|gi (x)|2 + |2 sin(x/2)|2di

X

|x + 2lπ|−2di |gi (x + 2lπ)|2 ,

(C.14)

l6=0

for di ≥ 1/2. Then, fi∗ (0) = |gi (0)|2 , and when −1/2 ≤ di < 1/2
|fi∗ (x) − fi∗ (0)|
≤ |gi (x)|2

2 sin(x/2)
x

2di +2

−1 + |gi (x)|2 −|gi (0)|2 + 2 sin(x/2)2di +2

X

|x+2lπ|−2di −2 |gi (x+2lπ)|2

l6=0

(C.15)
2

β

2di +2

= O(|x| ) + O(|x| ) + O(|x|

),

x → 0.

Similarly, when di ≥ 1/2,
|fi∗ (x) − fi∗ (0)| = O(|x|2 ) + O(|x|β ) + O(|x|2di ),
So, |fi∗ (x) − fi∗ (0)| < C|x|β∗ for x ∈ [−π, π), where
β∗ = min{β, 2d1 + 2},
48

d1 < 1/2,

x → 0.

and
β∗ = min{β, 2d1 },

d1 ≥ 1/2.

Thus, by Theorem 2 in Moulines et al. (2007a),
r

ν
d
2di f
j
j
f
a(ν) (WX (a(ν)2 )ii − EWX (a(ν)2 )ii )
→ N (0, W (di )),
a(ν)
j=j1 ,...,jm
The (l, l0 )-th entry of the limiting covariance matrix is given by
Z π
|D|jl −jl0 | (x, di )|2 dx,
Wl,l0 (di ) = 4π|gi (0)|4 24di max(jl ,jl0 )+min(jl ,jl0 )

i = 1, . . . , n.

l, l0 = 1, . . . , m,

−π

where D|jl −jl0 | (x; di ) is defined in (C.12). Moreover, the entries Xi , i = 1, . . . , n, of X are independent, thus (C.11) holds. 
The following proposition justifies the claim made in Remark 4.2.
e j be defined in (4.7). Then, for large enough J1 and J2 , J1 < J2 the
Proposition C.4 Let Λ
−1
eJ Λ
e has pairwise distinct diagonal entries.
matrix Λ
2

J1

eJ Λ
e −1
Proof of Proposition C.4: Reexpressing Λ
2 J1
2(J2 −J1 )d1
eJ Λ
e −1
Λ
, . . . , 22(J2 −J1 )dn )
2 J1 = diag(2


Z
Z
−J2
J2 2
−2d1 ∗
J2 2
−J1
·diag 2
|HJ2 (x/2 )| |x|
|g1 (x/2 )| dx 2
HJ1 (x/2J1 )|2 |x|−2d1 |g1∗ (x/2J1 )|2 dx, . . . ,
R
−J2

R

Z

J2

2

−2dn

|HJ2 (x/2 )| |x|

2

|gn∗ (x/2J2 )|2 dx



−J1

Z

2

HJ1 (x/2 )| |x|

2

R

J1

−2dn

|gn∗ (x/2J1 )|2 dx

R

By Theorem 1 (a) of Moulines et al. (2007b),
Z
Z
−j
j 2
−2di ∗
j 2
2
b
2
|Hj (x/2 )| |x|
|gi (x/2 )| dx →
|x|−2di |ψ(x)|
|gi (0)|2 dx,
R

j → ∞.

R

Thus, for i = 1, . . . , n,
Z
Z
|HJ2 (x/2J2 )|2 |x|−2di |gi∗ (x/2J2 )|2 dx
|HJ1 (x/2J1 )|2 |x|−2di |gi∗ (x/2J1 )|2 dx → 1,
R


.

J1 , J2 → ∞.

R

The claim holds as a consequence of condition (2.11).
Proof of Theorem 4.3: The proof can be written as a direct adaptation of Theorem 3.3 by
using Theorem 1 in Moulines et al. (2007b) as the counterpart of Proposition B.3. 

D

Proofs and auxiliary results: Section 6

Proof of Theorem 6.1: For any matrix S ∈ S+ (n, R), define the vector-valued function
f : vecS (S) → (ξ1 , . . . , ξn , vec(O))

(D.1)

such that S = Odiag(ξ1 , . . . , ξn )O∗ , O ∈ O(n), ξ1 < . . . < ξn , is the spectral decomposition of
S, and O = (oi1 i2 )i1 ,i2 =1,...,n satisfies o1i ≥ 0, i = 1, . . . , n (cf. (6.2)). Since EW (2j ) has pairwise
49

distinct eigenvalues, Theorem E.1 implies that f is infinitely differentiable on a neighborhood
of EW (2j ). Moreover, the Jacobian matrix Jj of f at the point EW (2j ) is given by (6.4) with
S = EW (2j ). So, let J = diag(J1 , . . . , Jm ). Recall the notation (2.1) for block-diagonal matrices.
The Delta method and Theorem 3.1, imply that
p
p
bj − Oj )))T
( Kj (vecD (Lj − Λj )), Kj (vec(O
j=j1 ,...,jm
p
d
= ( Kj (f (vecS (W (2j )) − f (vecS (EW (2j ))))Tj=j1 ,...,jm → Nmn(n+1) (0, JF J ∗ ),

(D.2)

as claimed. 

E

Useful results

Lemma E.1 (Moulines et al. (2007a), Lemma 4) Let {ξν , ν ≥ 1} be a sequence of centered Gaussian vectors and let Γν be the covariance matrix of ξν . Let (Aν )ν≥1 be a sequence of deterministic
matrices with adapted dimensions such that
lim Var(ξνT Aν ξν ) = σ 2 ∈ [0, ∞].

ν→∞

Assume that
lim ρ(Aν )ρ(Γν ) = 0,

ν→∞

where ρ(·) denotes the spectral radius. Then
L

ξνT Aν ξν − E(ξνT Aν ξν ) → N (0, σ 2 ).
Lemma E.2 (Moulines et al. (2007a), Lemma 6) Let m ≥ 2 be an integer and Γ be a m × m
covariance matrix. Let p be an integer between 1 and m − 1. let Γ1 be the top left submatrix with
size p × p and Γ2 the bottom right submatrix with size (m − p) × (m − p). Then
ρ(Γ) ≤ ρ(Γ1 ) + ρ(Γ2 ).
Lemma E.3 (Moulines et al. (2007a), Lemma 5) Let {ξk , k ∈ Z} be a stationary process with
spectral density function f and let Γν be the covariance matrix of (ξ1 , . . . , ξν ). Then, ρ(Γν ) ≤
2π k f k∞ .
The following theorem provides the partial derivatives of the eigenvalues and eigenvectors of
a symmetric matrix with respect to the latter.
Theorem E.1 (Magnus (1985), Theorem 1) Let S0 ∈ S(n, R), and let u0 be a normalized eigenvector associated with a simple eigenvalue λ0 of S0 . Then, we can define a real-valued and a vector
function λ and u, respectively, for all symmetric matrix S in some neighborhood N (S0 ) ∈ S(n, R)
of S0 , where
λ(S0 ) = λ0 ,
u(S0 ) = u0 ,
and
Su = λu,

uT u = 1,

50

S ∈ S(n, R).

Moreover, the functions λ and u are infinitely differentiable on N (S0 ), and their differentials at
S0 are given by
∂λ
= (uT0 ⊗ uT0 )D,
∂[vecS (S)]

∂u
= [uT0 ⊗ (λ0 In − S0 )+ ]D.
∂[vecS (S)]

(E.1)

In (E.1), the symbol ⊗ and the superscript + denote the Kronecker product and the Moore-Penrose
inverse, respectively, and D is the duplication matrix defined by (6.3).
Lemma
E.4 (Abry and Didier (2017), Lemma B.3) Let {φ. } ∈ R be a sequence such that
P∞
|φ
z=−∞ zgcd(aj ,aj 0 ) | < ∞. Then,
a ν

aj ν
j0
∞
X
1XX
φaj k−aj 0 k0 → gcd(aj , aj 0 )
φzgcd(aj ,aj 0 ) ,
ν
0
z=−∞

ν → ∞.

k=1 k =1

F

Repeated eigenvalues

Following up on the discussion in Remark 6.2, the next proposition describes the limiting distribution for the eigenvalues of W (2j ) for a special case where EW (2j ) has one repeated eigenvalue.
In its statement, we use the multivariate gamma function Γq (·), which is defined by
Γq (t) = π q(q−1)/4

q 
Y
i=1


1
t − (i − 1) .
2

Moreover, we replace (A1) with the following assumption.
Assumption (A10 ): the observed process has the mixed form Y = P X, where P is nonsingular,
X is defined in (1.3) and satisfy
d1 = d2 = . . . = dn =: d,

d > 1/2,

(F.1)

and the high frequency functions gi (x), i = 1, . . . , n are constants, i.e.,
g1 (x) = g1 , . . . , gn (x) = gn .
Proposition F.1 Suppose the assumptions (A10 –A2) hold. Let
EW (2j ) = OΛO∗ ,

b O
b∗ ,
W (2j ) = OL

(F.2)

be the matrix spectral decompositions of the wavelet and sample wavelet variance matrices, respectively. Assume the diagonal matrix Λ has the form


Λ1
0
Λ=
(F.3)
0 λ∗ Iq
for some 1 < q < n, where the main diagonal entries of the matrix Λ1 are pairwise distinct and
less than λ∗ . Let
L = diag(l1 , . . . , ln ), Λ1 = diag(λ1 , . . . , λn−q ).
(F.4)
Then, as ν → ∞,
p
T d
Kj (l1 − λ1 , . . . , ln−q − λn−q ), (ln−q+1 − λ∗ , . . . , ln − λ∗ ) → (LT1 , LT2 )T ,
51

(F.5)

where L1 and L2 are independent random vectors. Moreover,
L1 ∼ N (0, 2b diag(λ21 , . . . , λ2n−q )),
where
b :=

∞ nZ
X

b j x)|2 e−i2
|ψ(2

j zx

|x|−2d dx

.Z

b j x)|2 |x|−2d dx
|ψ(2

o2

,

(F.7)

R

R

z=−∞

(F.6)

and L2 has density
− 12 q

2

 
n
√
1
−1 q
( bλ∗ π)q(q−1)/4 Γq 2
exp − √
2
2 bλ∗

n
X

a2i

i=n−q+1

oY
(ai − al ),

(F.8)

l<i

where
ai = li − λ∗ ,

i = n − q + 1, . . . , n.

(F.9)

b be as in expression (F.2), and define
Proof: Let O and O
p
T = O∗ W (2j )O, U = Kj (T − Λ),

(F.10)

where O is the orthogonal matrix in the expression (F.2). Then, we can write
T = Y LY ∗ ,

b ∈ O(n),
Y = O∗ O

(F.11)

and thus
U = O∗

p
Kj (W (2j ) − EW (2j ))O.

(F.12)

Let d be as in (F.1). From (3.8), we obtain
j

Λ=2 O

∗

P diag(g12 , . . . , gn2 )P ∗ O

Z

b j x)|2 |x|−2d dx.
|ψ(2

R

For z ∈ Z, let Φz be as in (B.3) (for j = j 0 ). Under the condition (F.1),
Z
∗
∗
2
2
∗
b j x)|2 eizx |x|−2d dx
O Φz O = O P diag(g1 , . . . , gn )P O
|ψ(2
R

= 2−j Λ

nZ

b j x)|2 e−izx |x|−2d dx
|ψ(2

.Z

R

o
b j x)|2 |x|−2d dx .
|ψ(2

R

By (B.2) (which also holds under (F.1)),
Kj Kj
p p
1 1 XX ∗
Kj Kj
O ED(2j , k)D(2j , k 0 )∗ O ⊗ O∗ ED(2j , k)D(2j , k 0 )∗ O
Kj Kj
0
k=1 k =1

→ 22j

∞
X

O∗ Φz2j O ⊗ O∗ Φz2j O = b(Λ ⊗ Λ),

ν → ∞,

(F.13)

z=−∞

where the scalar b is given by (F.7). Thus, from (F.12),
d

U → U = {ui1 i2 }i1 ,i2 =1,...,n ,

52

(F.14)

where (vecS (U))T ∼ N n(n+1) (0, Ω) and Ω can be retrieved from (F.13) by means of (B.1). In
2

particular, all entries of (vecS (U))T are independent. Moreover, for λ• as in (F.4),
Var(ui1 i1 ) = 2b λ2i1 ,

Var(ui1 i2 ) = b λi1 λi2 ,

Var(ui1 i1 ) = 2b λ2∗ ,

Var(ui1 i2 ) = b λ2∗ ,

1 ≤ i1 , i2 ≤ n − q,

(F.15)

n − q + 1 ≤ i1 , i2 ≤ n

(F.16)

(the remaining entries of U will not play a role in the ensuing development). It now suffices to
follow the same arguments as in Sections 13.5.1 and 13.5.2 of Anderson (2003). For the reader’s
convenience, we lay out the main steps. Recast the random matrices T , Y , U and L in (F.10)
and (F.11) as






T11 T12
Y11 Y12
U11 U12
T =
, Y =
, U=
, L = diag(L1 , L2 ), (F.17)
T21 T22
Y21 Y22
U21 U22
where T11 , Y11 , U11 , L1 ∈ M (n − q, R), and let
p
A = Kj (L − Λ) = diag(A1 , A2 ).
Define
Y22 = EJF,

C2 = EF ∈ O(q),

(F.18)

where the first relation is a singular value decomposition, J is diagonal and E, F ∈ O(q) are
orthogonal. Also let
p
p
p
p
W11 = Kj (Y11 − I), W12 = Kj Y12 , W21 = Kj Y21 , W22 = Kj (Y22 − C2 ). (F.19)
Based on (F.17) and (F.19), we can reexpress the system of equalities T = Λ + √1 U = Y LY ∗
Kj
as



 



1
1
Λ1
U11 U12
In−q
W11 W12
T =
+p
=
+p
λ∗ Iq
U21 U22
C2
W21 W22
Kj
Kj

Λ1
·
=



 


∗
∗
1
1
A1
In−q
W11
W21
p
+p
·
+
∗
∗
W12
W22
A2
C2∗
Kj
Kj



 
1
A1
W11 Λ1 λ∗ W12 C2∗
+p
+
λ∗ Iq
W21 Λ1 λ∗ W22 C2∗
C2 A2 C2∗
Kj


 
∗
∗
1
Λ1 W11
Λ1 W21
+
.
(F.20)
+ OP
∗
∗
λ∗ C2 W12 λ∗ C2 W22
Kj


λ∗ Iq

Λ1

On the other hand, I = Y Y ∗ and the relations (F.19) yield



 

 
∗
∗
1
1
In−q
W11 W12 C2∗
W11
W21
In =
+p
+
+ OP
. (F.21)
∗
∗
∗
Iq
W21 W22 C2
C2 W12 C2 W22
Kj
Kj
From (F.20) and (F.21), we obtain the system of equations




1
1
∗
∗
U11 = W11 Λ1 + A1 + Λ1 W11 + OP p
, 0 = W11 + W11 + OP p
,
Kj
Kj


1
∗
U22 = C2 A2 C2 + OP p
.
Kj
53

(F.22)

(F.23)

Recall that the limiting joint distribution of (U11 , U22 ) is given by U11 := {ui1 i2 }i1 ,i2 =1,...,n−q and
U22 := {ui1 i2 }i1 ,i2 =n−q+1,...,n from expression (F.14), where
U11 and U22 are independent.

(F.24)

By following the same argument as on pp. 546 and 547 in Anderson (2003), expressions (F.22)
can be used to show that the limiting distribution of the diagonal entries of D1 is (F.6). Next
note that A2 and Y22 are functions of U depending on ν (see (F.10) and (F.11)), and C2 , in turn,
is a function of Y22 depending on ν (see (F.18)). Therefore, by the same argument as in Anderson
(2003), p. 549, the limiting distribution of A2 and C2 is the distribution of A2 and Y22 defined by
the expression
∗
U22 = Y22 A2 Y22
.
In particular, the limiting distribution of the diagonal entries of A2 is (F.8). In view of (F.24),
the established limiting distributions for the diagonal entries of A1 and A2 yield (F.5).
Example F.1 For n = 3, consider the OFBM for which d1 = d2 = d3 =:R d, P ∈ O(3), and
2 |y|−2d dy <
b
0 < g1 < g2 = g3 . Then, by (3.8), the eigenvalues of E(2j ) are λ1 = 22jd g12 R |ψ(y)|
R
2 |y|−2d dy, where the latter has multiplicity 2. Now let l ≤ l ≤ l be the
b
λ∗ = 22jd g22 R |ψ(y)|
1
2
3
ordered eigenvalues of the sample wavelet variance W (2j ) (cf. (6.2)). Then, by Proposition F.1,
p
T d
Kj l1 − λ1 , l2 − λ∗ , l3 − λ∗ → (LT1 , LT2 )T , ν → ∞.
(F.25)
In (F.25), L1 is independent of L2 , L1 ∼ N (0, 2b λ21 ) and L2 has density
n
o
1
1 √
−1
( bλ∗ π)1/2 Γ2 2 (1) exp − √
(a22 + a23 ) (a3 − a2 ),
2
2 bλ∗
where ai = li − λ∗ , i = 2, 3, and b is given by (F.7).

References
Abry, P. & Didier, G. (2017), ‘Wavelet estimation for operator fractional Brownian motion’, To
appear in Bernoulli pp. 1–30.
Achard, S. & Gannaz, I. (2016), ‘Multivariate wavelet Whittle estimation in long-range dependence’, Journal of Time Series Analysis 37(4), 476–512.
Amblard, P., Coeurjolly, J.-F., Lavancier, F. & Philippe, A. (2012), ‘Basic properties of the
multivariate fractional Brownian motion’, Bulletin de la Société Mathématique de France,
Séminaires et Congrés 28, 65–87.
Amblard, P.-O. & Coeurjolly, J.-F. (2011), ‘Identification of the multivariate fractional Brownian
motion’, IEEE Transactions on Signal Processing 59(11), 5152–5168.
Anderson, T. (2003), An Introduction to Multivariate Statistical Analysis, 3 edn, Wiley.
Bardet, J.-M. (2002), ‘Statistical study of the wavelet analysis of fractional Brownian motion’,
IEEE Transactions on Information Theory 48(4), 991–999.
Bardet, J.-M. & Tudor, C. (2010), ‘A wavelet analysis of the Rosenblatt process: chaos expansion
and estimation of the self-similarity parameter’, Stochastic Processes and their Applications
120(12), 2331–2362.
54

Bardet, J.-M. & Tudor, C. (2014), ‘Asymptotic behavior of the Whittle estimator for the increments of a Rosenblatt process’, Journal of Multivariate Analysis 131, 1–16.
Bardet, J.-M., Lang, G., Oppenheim, G., Philippe, A., Stoev, S. & Taqqu, M. S. (2003), Semiparametric estimation of the long-range dependence parameter: a survey, in ‘Theory and
Applications of Long-Range Dependence’, pp. 557–577.
Becker-Kern, P. & Pap, G. (2008), ‘Parameter estimation of selfsimilarity exponents’, Journal of
Multivariate Analysis (99), 117–140.
Belouchrani, A., Abed-Meraim, K., Cardoso, J.-F. & Moulines, E. (1997), ‘A blind source separation technique using second-order statistics’, IEEE Transactions on Signal Processing
45(2), 434–444.
Beran, J. (1994), Statistics for Long-memory Processes, Vol. 61, CRC press.
Beran, J., Feng, Y., Ghosh, S. & Kulik, R. (2013), Long memory processesprobabilistic properties
and statistical models, Springer, Heidelberg.
Biermé, H., Meerschaert, M. M. & Scheffler, H.-P. (2007), ‘Operator scaling stable random fields’,
Stochastic Processes and their Applications 117(3), 312–332.
Cardoso, J.-F. (1998), ‘Blind signal separation: statistical principles.’, Proceedings of the IEEE
86(10), 2009–2025.
Chan, K. & Tsai, K. (2010), ‘Inference of bivariate long-memory aggregate time series’, available
at http://public.econ.duke.edu/ brossi/NBERNSF/Tsai.pdf pp. 1–17.
Cheridito, P., Kawaguchi, H. & Maejima, M. (2003), ‘Fractional Ornstein-Uhlenbeck processes’,
Electronic Journal of Probability 8(3), 1–14.
Choi, S., Cichocki, A., Park, H.-M. & Lee, S.-Y. (2005), ‘Blind source separation and independent
component analysis: A review’, Neural Information Processing-Letters and Reviews 6(1), 1–
57.
Ciuciu, P., Abry, P. & He, B. (2014), ‘Interplay between functional connectivity and scale-free
dynamics in intrinsic fmri networks’, Neuroimage 95, 248–263.
Clausel, M. & Vedel, B. (2011), ‘Explicit construction of operator scaling Gaussian random fields’,
Fractals 19(01), 101–111.
Clausel, M. & Vedel, B. (2013), An optimality result about sample path properties of operator
scaling Gaussian random fields, in ‘Proceedings of the “XIème Colloque Franco-Roumain de
Mathématiques Appliquées” 4 (LXII)’, pp. 375–409.
Clausel, M., Roueff, F., Taqqu, M. S. & Tudor, C. (2014a), ‘Asymptotic behavior of the quadratic
variation of the sum of two Hermite processes of consecutive orders’, Stochastic Processes
and their Applications 124(7), 2517–2541.
Clausel, M., Roueff, F., Taqqu, M. S. & Tudor, C. (2014b), ‘Wavelet estimation of the long
memory parameter for Hermite polynomial of Gaussian processes’, ESAIM: Probability and
Statistics 18, 42–76.

55

Coeurjolly, J.-F., Amblard, P.-O. & Achard, S. (2013), ‘Wavelet analysis of the multivariate
fractional Brownian motion’, ESAIM: Probability and Statistics 17, 592–604.
Comon, P. & Jutten, C. (2010), Handbook of Blind Source Separation: Independent Component
Analysis and Applications, Academic Press.
Cryer, J. & Chan, K. (2008), Time Series Analysis with Application in R, Springer.
Dahlhaus, R. (1989), ‘Efficient parameter estimation for self-similar processes’, Annals of Statistics
pp. 1749–1766.
Daubechies, I. (1992), Ten Lectures on Wavelets, Vol. 61, Society for Industrial and Applied
Mathematics, Philadelphia-PA.
Didier, G. & Pipiras, V. (2011), ‘Integral representations and properties of operator fractional
Brownian motions’, Bernoulli 17(1), 1–33.
Didier, G. & Pipiras, V. (2012), ‘Exponents, symmetry groups and classification of operator
fractional Brownian motions’, Journal of Theoretical Probability 25, 353–395.
Didier, G., Helgason, H. & Abry, P. (2015), Demixing multivariate-operator selfsimilar processes,
in ‘IEEE International Conference on Acoustics, Speech and Signal, Processing (ICASSP),
Brisbane, Australia’, pp. 1–5.
Didier, G., Meerschaert, M. M. & Pipiras, V. (2017a), ‘Domain and range symmetries of operator
fractional brownian fields’, To appear in Stochastic Processes and their Applications.
Didier, G., Meerschaert, M. M. & Pipiras, V. (2017b), ‘Exponents of operator self-similar random
fields’, Journal of Mathematical Analysis and Applications 448(2), 1450–1466.
Dobrushin, R. & Major, P. (1979), ‘Non-central limit theorems for non-linear functional of Gaussian fields’, Probability Theory and Related Fields 50(1), 27–52.
Dogan, M., Van Dam, R., Liu, G., Meerschaert, M. M., Butler, J. J., Bohling, G., Benson, D. A.
& Hyndman, D. W. (2014), ‘Predicting flow and transport in highly heterogeneous alluvial
aquifers’, Geophysical Research Letters 41(21), 7560–7565.
Embrechts, P. & Maejima, M. (2002), Selfsimilar Processes, Princeton Series in Applied Mathematics, Princeton University Press, Princeton, NJ.
Fevotte, C. & Godsill, S.-J. (2006), ‘A bayesian approach for blind separation of sparse sources’,
IEEE Transactions on Audio, Speech, and Language Processing 14(6), 2174–2188.
Flandrin, P. (1992), ‘Wavelet analysis and synthesis of fractional brownian motion’, IEEE Transactions on Information Theory 38, 910 – 917.
Fontugne, R., Abry, P., Fukuda, K., Veitch, D., Cho, K., Borgnat, P. & Wendt, H. (2017),
‘Scaling in Internet traffic: a 14 year and 3 day longitudinal study, with multiscale analyses
and random projections’, IEEE/ACM Transactions on Networking 25(4), 1–14.
Foufoula-Georgiou, E. & Kumar, P. (2014), Wavelets in Geophysics, Vol. 4, Academic Press.
Fox, R. & Taqqu, M. (1986), ‘Large-sample properties of parameter estimates for strongly dependent stationary Gaussian time series’, Annals of Statistics pp. 517–532.
56

Frecon, J., Didier, G., Pustelnik, N. & Abry, P. (2016), ‘Non-linear wavelet regression and branch
and bound minimization for the full identification of bivariate operator fractional Brownian
motion’, IEEE Transactions on Signal Processing 64(15), 4040–4049.
Golub, G. & Van Loan, C. (2012), Matrix Computations, 4 edn, JHU Press.
Granger, C. & Joyeux, R. (1980), ‘An introduction to long-memory time series models and fractional differencing’, Journal of Time Series Analysis 1(1), 15–29.
Guo, H., Lim, C. Y. & Meerschaert, M. M. (2009), ‘Local Whittle estimator for anisotropic
random fields’, Journal of Multivariate Analysis 100(5), 993–1028.
Horn, R. & Johnson, C. (1985), Matrix Analysis, Cambridge University Press.
Hosking, J. R. M. (1981), ‘Fractional differencing’, Biometrika pp. 165–176.
Hosoya, Y. (1996), ‘The quasi-likelihood approach to statistical inference on multiple time-series
with long-range dependence’, Journal of Econometrics 73, 217–236.
Hosoya, Y. (1997), ‘A limit theory for long-range dependence and statistical inference on related
models’, The Annals of Statistics pp. 105–137.
Hualde, J. & Robinson, P. (2010), ‘Semiparametric inference in multivariate fractionally cointegrated systems’, Journal of Econometrics 157(2), 492–511.
Hudson, W. & Mason, J. (1982), ‘Operator-self-similar processes in a finite-dimensional space’,
Transactions of the American Mathematical Society 273(1), 281–297.
Ivanov, P., Nunes Amaral, L., Goldberger, A., Havlin, S., Rosenblum, M., Struzik, Z. & Stanley,
H. (1999), ‘Multifractality in human heartbeat dynamics’, Nature 399(6735), 461–465.
Jones, G. A. & Jones, J. M. (1998), Elementary Number Theory, Berlin: Springer-Verlag.
Kechagias, S. & Pipiras, V. (2015a), ‘Definitions and representations of multivariate long-range
dependent time series’, Journal of Time Series Analysis 36(1), 1–25.
Kechagias, S. & Pipiras, V. (2015b), ‘Identification, estimation and applications of a bivariate
long-range dependent times series model with general phase’, Preprint pp. 1–25.
Laha, R. G. & Rohatgi, V. K. (1981), ‘Operator self-similar stochastic processes in Rd ’, Stochastic
Processes and their Applications 12(1), 73–84.
Li, Y. & Xiao, Y. (2011), ‘Multivariate operator-self-similar random fields’, Stochastic Processes
and their Applications 121(6), 1178–1200.
Li, Y., Adali, T., Wang, W. & Calhoun, V. (2009), ‘Joint blind source separation by multiset
canonical correlation analysis’, IEEE Transactions on Signal Processing 57(10), 3918–3929.
Lobato, I. (1997), ‘Consistency of the averaged cross-periodogram in long memory series’, Journal
of Time Series Analysis 18(2), 137–155.
Maejima, M. & Mason, J. (1994), ‘Operator-self-similar stable processes’, Stochastic Processes
and their Applications 54, 139–163.

57

Magnus, J. (1985), ‘On differentiating eigenvalues and eigenvectors’, Econometric Theory
1(2), 179–191.
Magnus, J. & Neudecker, H. (1980), ‘The elimination matrix: some lemmas and applications’,
SIAM Journal on Algebraic Discrete Methods 1(4), 422–449.
Mallat, S. (1999), A Wavelet Tour of Signal Processing, Academic Press.
Mandelbrot, B. (1974), ‘Intermittent turbulence in self-similar cascades: divergence of high moments and dimension of the carrier’, J. Fluid Mech. 62, 331–358.
Mandelbrot, B. & Van Ness, J. (1968), ‘Fractional brownian motions, fractional noises and applications’, SIAM Review 10(4), 422–437.
Marinucci, D. & Robinson, P. M. (2000), ‘Weak convergence of multivariate fractional processes’,
Stochastic Processes and their Applications 86(1), 103–120.
Marinucci, D. & Robinson, P. M. (2001), ‘Semiparametric fractional cointegration analysis’, Journal of Econometrics 105(1), 225–247.
Mason, J. & Xiao, Y. (2002), ‘Sample path properties of operator-self-similiar Gaussian random
fields’, Theory of Probability and Its Applications 46(1), 58–78.
Masry, E. (1993), ‘The wavelet transform of stochastic processes with stationary increments and
its application to fractional Brownian motion’, IEEE Transactions on Information Theory
39(1), 260–264.
Moreau, E. (2001), ‘A generalization of joint-diagonalization criteria for source separation’, IEEE
Transactions on Signal Processing 49(3), 530–541.
Moulines, E., Roueff, F. & Taqqu, M. (2007a), ‘Central limit theorem for the log-regression
wavelet estimation of the memory parameter in the gaussian semi-parametric context’, Fractals 15(4), 301–313.
Moulines, E., Roueff, F. & Taqqu, M. (2007b), ‘On the spectral density of the wavelet coefficients
of long-memory time series with application to the log-regression estimation of the memory
parameter’, Journal of Time Series Analysis 28(2), 155–187.
Moulines, E., Roueff, F. & Taqqu, M. (2008), ‘A wavelet whittle estimator of the memory parameter of a nonstationary Gaussian time series’, Annals of Statistics pp. 1925–1956.
Nielsen, F. S. (2011), ‘Local Whittle estimation of multi-variate fractionally integrated processes’,
Journal of Time Series Analysis 32(3), 317–335.
Nielsen, M. & Frederiksen, P. (2011), ‘Fully modified narrow-band least squares estimation of
weak fractional cointegration’, The Econometrics Journal 14(1), 77–120.
O’Grady, P., Pearlmutter, B. & Rickard, S. (2005), ‘Survey of sparse and non-sparse methods in
source separation’, International Journal of Imaging Systems and Technology 15(1), 18–33.
Parra, L. & Sajda, P. (2003), ‘Blind source separation via generalized eigenvalue decomposition’,
Journal of Machine Learning Research 4, 1261–1269.

58

Pham, D.-T. & Cardoso, J.-F. (2001), ‘Blind separation of instantaneous mixtures of nonstationary sources’, IEEE Transactions on Signal Processing 49(9), 1837–1848.
Pipiras, V. & Taqqu, M. S. (2017), Long-Range Dependence and Self-Similarity, Cambridge Series on Statistical and Probabilistic Mathematics, Cambridge University Press, Cambridge,
United Kingdom.
Prakasa Rao, B. L. S. (2010), Statistical Inference for Fractional Diffusion Processes, Wiley Series
in Probability and Statistics.
Puplinskaitė, D. & Surgailis, D. (2015), ‘Scaling transition for long-range dependent Gaussian
random fields’, Stochastic Processes and their Applications 125(6), 2256–2271.
Robinson, P. (1995a), ‘Gaussian semiparametric estimation of long range dependence’, Annals of
Statistics 23(5), 1630–1661.
Robinson, P. (1995b), ‘Log-periodogram regression of time series with long range dependence’,
Annals of Statistics pp. 1048–1072.
Robinson, P. (2008), ‘Multivariate Local Whittle estimation in stationary systems’, Annals of
Statistics 36(5), 2508–2530.
Robinson, P. M. & Yajima, Y. (2002), ‘Determination of cointegrating rank in fractional systems’,
Journal of Econometrics 106(2), 217–241.
Sela, R. & Hurvich, C. (2012), ‘The averaged periodogram estimator for a power law in coherency’,
Journal of Time Series Analysis 33(2), 340–363.
Shimotsu, K. (2007), ‘Gaussian semiparametric estimation of multivariate fractionally integrated
processes’, Journal of Econometrics 137, 277–310.
Shimotsu, K. (2012), ‘Exact local Whittle estimation of fractionally cointegrated systems’, Journal
of Econometrics 169(2), 266–278.
Stone, J. (2004), Independent Component Analysis: a Tutorial Introduction, MIT Press.
Taqqu, M. S. (1975), ‘Weak convergence to fractional Brownian motion and to the Rosenblatt
process’, Probability Theory and Related Fields 31(4), 287–302.
Taqqu, M. S. (1979), ‘Convergence of integrated processes of arbitrary Hermite rank’, Probability
Theory and Related Fields 50(1), 53–83.
Taqqu, M. S. (2003), Fractional Brownian motion and long range dependence, in ‘Theory and
Applications of Long-Range Dependence (P. Doukhan, G. Oppenheim and M. S. Taqqu,
eds.)’, Birkhäuser, Boston, pp. 5–38.
Taqqu, M., Willinger, W. & Sherman, R. (1997), ‘Proof of a fundamental result in self-similar
traffic modeling’, ACM SIGCOMM Computer Communication Review 27(2), 5–23.
Tsai, H. & Chan, K. (2005), ‘Quasi-maximum likelihood estimation for a class of continuous-time
long-memory processes’, Journal of Time Series Analysis 26(5), 691–713.
Tsai, H., Rachinger, H. & Chan, K.-S. (2017), ‘Inference of bivariate long-memory aggregate time
series’, To appear in Statistica Sinica.
59

Veitch, D. & Abry, P. (1999), ‘A wavelet-based joint estimator of the parameters of long–range
dependence’, IEEE Transactions on Information Theory 45(3), 878–897.
Vignat, C. (2012), ‘A generalized Isserlis theorem for location mixtures of Gaussian random
vectors’, Statistics and Probability Letters 82(1), 67–71.
Wendt, H., Didier, G., Combrexelle, S. & Abry, P. (2017), ‘Multivariate Hadamard self-similarity:
testing fractal connectivity’, available at https://arxiv.org/pdf/1701.04366v1.pdf pp. 1–51.
Wendt, H., Scherrer, A., Abry, P. & Achard, S. (2009), Testing fractal connectivity in multivariate
long memory processes , in ‘IEEE International Conference on Acoustics, Speech and Signal
Processing, Taipei, Taiwan, April 19–24, 2009’, pp. 2913–2916.
Whitcher, B., Guttorp, P. & Percival, D. B. (2000), ‘Wavelet analysis of covariance with application to atmospheric time series’, Journal of Geophysical Research 105(D11), 14941–14962.
Wornell, G. & Oppenheim, A. (1992), ‘Estimation of fractal signals from noisy measurements
using wavelets’, IEEE Transactions on Signal Processing 40(3), 611–623.
Xiao, Y. (2009), Sample path properties of anisotropic Gaussian random fields, in ‘A Minicourse
on Stochastic Partial Differential Equations’, Springer, pp. 145–212.
Yaglom, A. (1958), ‘Correlation theory of processes with random stationary nth increments’,
American Mathematical Society Translations 8, 87–141.
Yeredor, A. (2002), ‘Non-orthogonal joint diagonalization in the least-squares sense with application in blind source separation’, IEEE Transactions on signal processing 50(7), 1545–1553.
Ziehe, A., Laskov, P., Nolte, G. & Müller, K.-R. (2004), ‘A fast algorithm for joint diagonalization
with non-orthogonal transformations and its application to blind source separation’, Journal
of Machine Learning Research 5, 777–800.
Patrice Abry
Univ Lyon, ENS de Lyon,
Univ Claude Bernard, CNRS,
Laboratoire de Physique,
F-69342 Lyon,
France
patrice.abry@ens-lyon.fr

Gustavo Didier and Hui Li
Mathematics Department
Tulane University
6823 St. Charles Avenue
New Orleans, LA 70118
USA
gdidier@tulane.edu
hli15@tulane.edu

60

