arXiv:1608.07839v1 [math.PR] 28 Aug 2016

Non-Linear Wavelet Regression and Branch & Bound
Optimization for the Full Identification of Bivariate
Operator Fractional Brownian Motion
Jordan Frecon∗, Gustavo Didier†, Nelly Pustelnik∗ , and Patrice Abry∗‡

Abstract
Self-similarity is widely considered the reference framework for modeling the scaling
properties of real-world data. However, most theoretical studies and their practical
use have remained univariate. Operator Fractional Brownian Motion (OfBm) was
recently proposed as a multivariate model for self-similarity. Yet it has remained seldom
used in applications because of serious issues that appear in the joint estimation of
its numerous parameters. While the univariate fractional Brownian motion requires
the estimation of two parameters only, its mere bivariate extension already involves
7 parameters which are very different in nature. The present contribution proposes
a method for the full identification of bivariate OfBm (i.e., the joint estimation of all
parameters) through an original formulation as a non-linear wavelet regression coupled
with a custom-made Branch & Bound numerical scheme. The estimation performance
(consistency and asymptotic normality) is mathematically established and numerically
assessed by means of Monte Carlo experiments. The impact of the parameters defining
OfBm on the estimation performance as well as the associated computational costs are
also thoroughly investigated.

∗ Jordan Frecon (Corresponding author), Nelly Pustelnik and Patrice Abry are with Univ Lyon,
Ens de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique, F-69342 Lyon, France (e-mail:
firstname.lastname@ens-lyon.fr).
† Gustavo Didier is with Mathematics Department, Tulane University, New Orleans, LA 70118, USA,
(e-mail: gdidier@tulane.edu).
‡ Work supported by French ANR grant AMATIS #112432, 2010-2014 and by the prime award no.
W911NF-14-1-0475 from the Biomathematics subdivision of the Army Research Office, USA. G.D. gratefully
acknowledges the support of Ens de Lyon.

1

1

Introduction

Scale invariance and self-similarity. Scale invariance, or scaling, is now recognized as
an ubiquitous property in a variety of real-world applications which are very different in
nature (cf. e.g., [25] and references therein for reviews). The so-named scale invariance
paradigm is based on the assumption that temporal dynamics in data are not driven by
one, or a few, representative time scales, but by a large continuum of them. Self-similar
stochastic processes provide the basal mathematical framework for the modeling of scaling
phenomena. In essence, self-similarity states that a signal X cannot be distinguished from
any of its dilated copies (cf. e.g., [22]):
f dd

{X(t)}t∈R = {aH X(t/a)}t∈R , ∀a > 0,

(1)

f dd

where = stands for the equality of finite dimensional distributions. The key information
on scale-free dynamics is summed up under a single parameter 0 < H < 1, called the Hurst
exponent, whose estimation is the main goal in scaling analysis. Amongst the numerous
estimators of H proposed in the literature (cf. e.g., [5] for a review), one popular methodology draws upon the computation of the sample variance of a set of multiscale quantities
(e.g., wavelet coefficients) TX (a, t) that behave like a power law with respect to the scale a:
X
2
TX
(a, t) ' a2H+1 .
(2)
t

In view of the relation (2), H can be estimated by means of a linear regression in log–log
coordinates (cf. e.g., [24]).
Fractional Brownian motion (fBm) BH (t) – i.e., the only Gaussian, self-similar, stationaryincrement process – has massively been used as a reference process in the modeling of scaling
properties in univariate real-world signals.
Multivariate scaling. Notwithstanding its theoretical and practical importance, fBm falls
short of providing an encompassing modeling framework for scaling because most modern
contexts of application involve the recording of multivariate time series that hence need to be
jointly analyzed. The construction of a comprehensive multivariate estimation paradigm is
still an open problem in the literature. The so-named Operator fractional Brownian motion
(OfBm), henceforth denoted by B W,H (t), is a natural extension of fBm. It was recently been
defined and studied in [3, 10, 9, 7] as the only Gaussian, multivariate self-similar process
with stationary increments. Multivariate self-similarity translates into the relation:
f dd

{B W,H (t)}t∈R = {aH B W,H (t/a)}t∈R , ∀a > 0,

(3)

where the scaling exponent consists of a Hurst matrix H = W diagHW −1 . In the latter
expression, W represents a P × PPinvertible matrix, and H is a P -dimensional vector of
+∞
k
k
Hurst eigenvalues, where aH :=
k=0 log (a)H /k!. The full parametrization of OfBm
further requires a P × P point-covariance matrix Σ. OfBm remains so far rarely used in
applications, mostly because its actual use requires, in a general setting, the estimation of
P + P 2 + P (P − 1)/2 parameters which are very different in nature (cf. Section 2). In
particular, Eq. (2) above results in a mixture of power laws (cf. Eqs. (13)-(15) in Section 2).
However, the identification of OfBm has been thoroughly studied in the entry-wise scaling
case (corresponding to a diagonal mixing matrix W ) [3] and often used in applications (cf.
2

e.g., [2, 6]). Identification has also been recently achieved under a non-diagonal mixing
matrix W , yet with more restrictive assumptions on Σ [8]. Even more recently, [1] proposed a general estimator for the vector of Hurst eigenvalues H in the bivariate setting,
yet requiring additional assumptions for the estimation of the extra parameters W and Σ.
The full identification of OfBm without parametric assumptions has remained an open issue.
Goals, contributions and outline. In this work, our contribution is two-fold. First,
the full identification of bivariate (P = 2) OfBm (Biv-OfBm) is formulated as a non-linear
wavelet domain regression. Second, an algorithmic solution for the associated optimization
problem is devised by means of a Branch & Bound procedure, which is essential in view
of the highly non-convex nature of the latter. To this end, definitions and properties of
Biv-OfBm are recapped in Section 2. A parsimonious parametrization of the process is
also proposed, which prevents the potential parametric under-determination of Biv-OfBm
[10]. In Section 3, the properties of the wavelet coefficients and of the wavelet spectrum of
Biv-OfBm are explicitly laid out and computed. This provides a mathematical framework
for the proposed estimation method. The full identification of Biv-OfBm is formulated as
a minimization problem whose solution is developed based on a Branch & Bound strategy
(cf. Section 4). The consistency and asymptotic normality of the proposed estimator is
mathematically established in the general multivariate setting P ≥ 2 (cf. Section 5), and
numerically assessed in the bivariate setting P = 2 by means of Monte Carlo experiments
conducted on large numbers of synthetic Biv-OfBm paths (cf. Section 6). Comparisons
with the Hurst eigenvalues estimators proposed in [1] are also reported. The routines for
the identification and synthesis of OfBm will be made publicly available at the time of
publication.

2
2.1
2.1.1

Bivariate Operator fractional Brownian motion
Definitions
Preamble

The most general definitions of OfBm were formulated in [3, 10, 9, 7] as the only multivariate
Gaussian, self-similar (i.e., satisfying Eq. (3)) process with stationary increments. Targeting
real-world data and applications, the present contribution is restricted to the (slightly)
narrower class of time reversible OfBm (cf. [3, 10, 9, 7]) whose scaling exponent matrix H
can be diagonalized as H = W diagHW −1 , where W is an invertible matrix. The definitions
and properties of OfBm are stated only in the bivariate setting.
2.1.2

Entry-wise scaling OfBm

The entry-wise scaling, time-reversible OfBm {X(t)}t∈R ≡ B Id,H (t) is defined by the condition W = Id. Hence, the Hurst exponent has the form H = diag(h1 , h2 ), 0 < h1 ≤ h2 < 1.
Let ΣX ≡ EX(1)X(1)∗ denote the point covariance matrix of X with entries σxm σxn ρxm ,xn ,
where σx2m is the variance of component m and ρxm ,xn is the correlation between components
m and n. It was shown in [4], [9] that the process X is well-defined (i.e., that its covariance

3

matrix EX(t)X(s)∗ is always positive definite) if and only if (with ρx ≡ ρx1 ,x2 ):
g(h1 , h2 , ρx ) ≡ Γ(2h1 + 1)Γ(2h2 + 1) sin(πh1 ) sin(πh2 )

− ρ2x Γ(h1 + h2 + 1)2 sin2 (π(h1 + h2 )/2) > 0. (4)

For entry-wise scaling OfBm, self-similarity in Eq. (3) simplifies to:
f dd

{X1 (at), X2 (at)}t∈R = {ah1 X1 (t), ah2 X2 (t)}t∈R , ∀a > 0.

(5)

The estimation of the parameters (h1 , h2 , ρx , σx1 , σx2 ), which fully characterize the process,
can thus be conducted by following univariate-type strategies, i.e., by making use of extensions of Eq. (2) to all auto- and cross-components (cf. [3, 7] for a theoretical study estimation
performance, or [6] for wavelet-based estimation on real-world data).
2.1.3

Mixing

Let W denote a 2 × 2 invertible matrix, hereinafter called the mixing matrix. OfBm is
defined as {B W,H (t) ≡ Y (t)}t∈R = {W B Id,H (t) ≡ W X(t)}t∈R . Following [3, 9], it is
straightforward to show that Y is self-similar as in Eq. (3), with H = W diagHW −1 . When
W is not diagonal, OfBm is no longer entry-wise scaling. Instead, the entry-wise scaling
behavior of OfBm consists of mixtures of univariate power laws (cf. Eqs. (13)-(15)). For
this reason, the construction of estimators in the bivariate setting cannot rely on a direct
extension of a univariate procedure.

2.2
2.2.1

Properties
Under-determination

Because W is invertible, one can show that, for Σy (t) ≡ EY (t)Y (t)∗ ,
Σy (t) = W EX(t)X(t)∗ W ∗ ≡ W Σx (t)W ∗ ,

(6)

which reveals three forms of under-determination in the parametrization of OfBm:
∗
i) Writing TX = diag(σx1 , σx2 ) and ΣX = TX CX TX
, where CX ≡ {1 ρx ; ρx 1} is the
correlation matrix of X, one cannot discriminate between Y = W X and Y = W 0 X 0 , where
−1
W 0 = W TX and X 0 = TX
X.
ii) Let Π denote a 2 × 2 permutation matrix, i.e., there is only one non-zero entry (equal
to 1) per column or line. Then, Y = W X = W 0 X 0 , where W 0 = W Π and X 0 = ΠT X.
iii) Let S be a diagonal matrix with entries ±1 and X 0 = SX, then Y = W X = W 0 X 0 ,
where W 0 ≡ W S.
2.2.2

Parametrization

To fix the parametric under-determination of OfBm, we adopt the following conventions:
i) the columns of W are normalized to 1; ii) h1 ≤ h2 ; iii) the diagonal entries of W are
positive. This leads us to propose the following generic 7-dimensional parametrization Θ =
(h1 , h2 , ρx , σx1 , σx2 , β, γ) of Biv-OfBm {Y (t)}t∈R :




√1 2 √β 2
σx21
σx1 σx2 ρx
1+γ
1+β


, ΣX =
.
(7)
W =
√−γ 2 √ 1 2
σx1 σx2 ρx
σx22
1+γ

1+β

4

3

Wavelet Analysis of OfBm

3.1

Multivariate discrete wavelet transform (DWT)

R
Let ψ0 be a mother wavelet, namely, ψ0 ∈ L2 (R) and R tk ψ0 (t)dt ≡ 0, k = 0, 1, . . . , Nψ − 1.
Let {ψj,k (t) = 2−j/2 ψ0 (2−j t−k)}(j,k)∈Z2 denote the collection of dilated and translated templates of ψ0 that forms an orthonormal basis of L2 (R). The multivariate DWT coefficients
of {Y (t)}t∈R are defined as (Dy (j, k)) ≡ (Dy1 (j, k), Dy2 (j, k)), where
Z
Dym (j, k) =
2−j/2 ψ0 (2−j t − k)Ym (t)dt, m = 1, 2.
(8)
R

For a detailed introduction to wavelet transforms, interested readers are referred to, e.g.,
[17].

3.2

Wavelet spectrum

The properties of the wavelet coefficients of OfBm in a P -variate setting were studied in
detail in [1]. Here, we only recall basic properties and expand on what is needed for actual
full identification (i.e., the estimation of all parameters entering its definition) of Biv-OfBm.
3.2.1

Mixture of power laws

From Eq. (5) and Y (t) = W X(t), it can be shown that the wavelet spectrum reads:
EDy (j, k)Dy (j, k)∗ = W 2j(H+Id/2) E0 2j(H

∗

+Id/2)

W∗

(9)

with E0 ≡ EDx (0, k)Dx (0, k)∗ =
σx21 ηh1
ρx σx1 σx2 η h1 +h2

ρx σx1 σx2 η h1 +h2
2
σx22 ηh2

2

!
(10)

Z
Z
1
|u|2h du ψ0 (v)ψ0 (v − u)∗ dv > 0.
(11)
2 R
R
The OfBm parametrization proposed above yields the following explicit form of the
wavelet spectrum:


E11 (2j , Θ) E12 (2j , Θ)
∗
j
EDy (j, k)Dy (j, k) ≡ E(2 , Θ) =
,
(12)
E12 (2j , Θ) E22 (2j , Θ)
and ηh = −

with E11 (2j , Θ) = (1 + γ 2 )−1 σx21 ηh1 2j(2h1 +1)
+ 2β(1 + β 2 )−1/2 (1 + γ 2 )−1/2 ρx σx1 σx2 η h1 +h2 2j(h1 +h2 +1)
2

+ β 2 (1 + β 2 )−1 σx22 ηh2 2j(2h2 +1) , (13)
E12 (2j , Θ) = −γ(1 + γ 2 )−1 σx21 ηh1 2j(2h1 +1)

+ (1 − βγ)(1 + β 2 )−1/2 (1 + γ 2 )−1/2 ρx σx1 σx2 η h1 +h2 2j(h1 +h2 +1)
2

+ β(1 + β 2 )−1 σx22 ηh2 2j(2h2 +1) , (14)
5

E22 (2j , Θ) = γ 2 (1 + γ 2 )−1 σx21 ηh1 2j(2h1 +1)
− 2γ(1 + β 2 )−1/2 (1 + γ 2 )−1/2 ρx σx1 σx2 η h1 +h2 2j(h1 +h2 +1)
2

+ (1 + β 2 )−1 σx22 ηh2 2j(2h2 +1) . (15)
3.2.2

Further under-determination

Eqs. (13), (14) and (15) reveal that E11 (2j , Θ), |E12 (2j , Θ)| and E22 (2j , Θ) are invariant
under the transformation (β, γ, ρx ) → −(β, γ, ρx ). Therefore, the definition of ρx can be
restricted to ρx ≥ 0.

3.3

Empirical wavelet spectrum

The goal is to estimate the Biv-OfBm parameters Θ = (h1 , h2 , ρx , σx1 , σx2 , β, γ) starting from
the wavelet spectrum EDy (j, ·)Dy (j, ·)∗ . The plug-in estimator of the ensemble variance
EDy (j, ·)Dy (j, ·)∗ is the sample variance
S(2j ) =

Kj
1 X
Dy (j, k)Dy (j, k)∗ ,
Kj

Kj =

k=1

N
,
2j

where N denotes the sample size. Fig. 1 illustrates the fact that S(2j ) is a satisfactory
estimator for EDy (j, ·)Dy (j, ·)∗ .

4
4.1

Non-linear regression based estimation and Branch
and Bound algorithm
Identification procedure as a minimization problem

The estimation of the parameter vector Θ of Biv-OfBm is challenging because its entry-wise
wavelet (or Fourier) spectrum is a mixture of power laws (cf. Eqs. (13)-(15)). This precludes
the direct extension of classical univariate techniques, based on the scalar relation Eq. (2)
[24]. For this reason, we formulate the full identification of Biv-OfBm (i.e., the estimation
of Θ) as a minimization problem1 :
Θ̂M
N = argmin CN (Θ), where

(16)

Θ∈Q0

CN (Θ) ≡

P
=2
X

j2
X

i1 ,i2 =1 j=j1
i1 ≤i2

2
log2 |Si1 ,i2 (2j )| − log2 |Ei1 ,i2 (2j , Θ)| .

(17)

The use of log2 ensures that the scales 2j , j = j1 , . . . , j2 contribute equally to CN . The search
space incorporates prior information in the shape of constraints: Sections 2.2.1 and 3.2.2
impose h1 ≤ h2 and ρx ∈ [0, 1] ; feasible solutions must satisfy the constraint g(h1 , h2 , ρx ) >
0 (cf. Eq. (4)) and (β, γ) ∈ [−1, 1]2 . For the sake of feasibility, we further restrict (σx1 , σx2 ) ∈
1 The superscript ·M has been added in order to refer to the M -estimator whose theoretical details are
given in Section 5.

6

3
30

3
30

3
30

20

2

20

2

20

2

10

1

10

1

10

1

0
1

0
5
10
j = log 2 scale

0
15

1

0
5
10
j = log 2 scale

0
15

1

5
10
j = log 2 scale

0
15

Figure 1: Wavelet Spectrum. Superimposition of log2 |Ep,p0 (2j , Θ)| (red ’+’) and
log2 |Sp,p0 (2j )| (solid black line) (with (p, p0 ) = (1, 1), (1, 2), (2, 2) from left to right) for
a single realization of Biv-ofBm with Θ = (h1 = 0.4, h2 = 0.8, ρx = 0.1, σx1 = 1, σx2 =
1, β = 0.5, γ = 0.5) (absolute difference between data and model is shown in dashed blue).
q
by21 + σ
by22 , where σ
by2m denotes the sample variance estimates of
[0, σmax ]2 , with σmax = σ
the increments of Ym . We arrive at the parameter space
n
Q0 = Θ = (h1 , h2 ,ρx , σx1 , σx2 , β, γ) ∈ R7 | Θ ∈ [0, 1]3 × [0, σmax ]2 × [−1, 1]2 ,
o
g(h1 , h2 , ρx ) > 0, h1 ≤ h2 .
(18)
The minimization of CN (Θ) is an intricate task for two reasons. First, because it involves disentangling a mixture of power laws, which yields a highly non-convex function.
Second, because the parameters to be estimated in Θ (scaling exponents, mixing coefficients,
variances and correlation) are very different in nature. The present contribution makes the
original proposition of searching the global minimum of Eq. (16) by means of a Branch &
Bound procedure detailed in the next section.

4.2

Global minimization via a Branch & Bound strategy

Branch & Bound algorithms consist of smart enumeration methods, which were shown to
solve a variety of constrained global non-convex optimization problems [14, 16, 18, 21]. In
the context of the estimation problem (16), it amounts to partitioning (branching) the search
space Q0 into smaller and smaller subregions, bounding the range of the objective function
CN in each subregion, and then identifying the region containing the global minimum. This
can be rephrased as 4 steps which are repeated until a stopping criterion is reached:
- Selecting: Choose any region R from the search space and relax it into a closed convex
set, i.e. an interval, as illustrated by the dashed line in Fig. 2 (left plot).
- Partitioning: Divide R into two smaller regions Ra and Rb .
- Bounding: Compute lower and upper bounds of CN on Ra and Rb . Upper bounds can be
obtained by evaluating CN anywhere in the region at hand. Lower bounds are computed by
resorting to interval arithmetic techniques (cf. Appendix A and [19, 15, 20]), which combine
elementary operations to produce rough lower bounds for the range of a given function, here
CN , on any interval.
- Pruning: Pruning is driven by three mechanisms: discard regions that do not satisfy
constraints (infeasibility) ; discard regions whose lower bound is larger than the smallest
upper bound as they cannot contain the global minimum (bound ) ; discard regions whose
size (for all parameters) has reached the targeted precision (size).

7

... ... ...
... ... ... ... ... ... C∆2
C2 C4 ...
... ...
C1 C3
... ...

Q0

Initialization
2
S0 = ∪∆
i=1 Ci

Initialization
Q0

Figure 2: Schematic view of the proposed approximation S0 ⊆ Q0 .

4.3
4.3.1

Branch & Bound procedure for Biv-OfBm identification
Convex relaxation

By nature, interval arithmetic techniques apply only to intervals, i.e., to convex sets. Therefore, in most Branch & Bound procedures, a convex relaxation of the search space Q0 is
required at initialization, as sketched in Fig. 2 (left). However, in the present case, a convex
relaxation of Q0 is not feasible because of the constraint g(h1 , h2 , ρx ) > 0.
Instead, we propose to approximate Q0 by an inner convexrelaxation S0 , consisting
2
of the union of ∆2 separable convex sets Ci , i.e., S0 = ∪∆
⊂ Q0 , as illustrated in
i=1 Ci

Fig. 2 (right). In practice, (h1 , h2 , ρx ) ∈ [0, 1]3 is approximated by a union of ∆2 nonoverlapping parallelepipedic sets {Pi }1≤i≤∆2 , denoted P ⊂ [0, 1]3 . They are obtained by
dividing (h1 , h2 ) ∈ [0, 1]2 into squares Ti with a discretization step ∆−1 and defining Pi =
Ti × [0, ρi ] where ρi is the largest value such that (∀(h1 , h2 ) ∈ Ti ), g(h1 , h2 , ρi ) > 0:
n
o
Ci = Θ = (h1 , h2 , ρx , σx1 , σx2 , β, γ) ∈ R7 | Θ ∈ Pi × [0, σmax ]2 × [−1, 1]2 , h1 ≤ h2 . (19)
In this inner convex relaxation strategy, the constraint g(h1 , h2 , ρx ) > 0 is necessarily
satisfied, a major practical benefit as infeasible regions need not be explored.
4.3.2

Algorithm

The full identification of Biv-OfBm is achieved via the following proposed sequence of operations:
Inputs.
- From data, compute the wavelet spectrum S(2j ), j = j1 , . . . , j2 ;
- Pick the Biv-OfBm model (i.e., Eqs. (13)-(15)) ;
- Set the precision δ for each parameter ;
- Set the inner convex relaxation parameter ∆ in order to approximate the set Q0 by
2
S0 = ∪∆
i=1 Ci ;
Initialization.
- Set Sb = ∅ and k = 0.
- Compute lower bounds li of CN on Ci (∀i = 1, . . . , ∆2 ) ;
- Compute upper bounds ui of CN on Ci (∀i = 1, . . . , ∆2 ) ;
- Set U = min(u1 , . . . , u∆2 ) and L = min(l1 , . . . , l∆2 ).
Iteration. Let Sk denote the partitioning at the step k.
- Selecting: Select region R ⊂ Sk with lowest lower bound L (best-first-search strategy) 2 .
2 At the step k = 0, this amounts in choosing one C . More generally, it consists in selecting one element
i
of the partitioning Sk .

8

- Cutting: Divide R into Ra and Rb such that R = Ra ∪ Rb and Ra ∩ Rb = ∅, along
its longest edge, in half, where the length of an edge is defined relatively to the maximum
accuracy δ prescribed by the practitioner.
- Lower bound: Compute lower bounds of Ra and Rb , using interval arithmetic.
- Upper bound: Compute upper bound of Ra (resp. Rb ), by evaluating CN (Θ) for Θ chosen
at the center of Ra (resp. Rb ).
- Branching: Update the partitioning Sk+1 = (Sk \R) ∪ Ra ∪ Rb , and update U and L on
this new partition.
- Pruning: Discard regions R∗ of Sk+1 either by bound, infeasibility or size, i.e. Sk+1 ←
Sk+1 \R∗ . Append to Sb the regions of Sk+1 discarded by size ; Discard regions in Sb for
bound.
- Set k ← k + 1
Stop and Output. Stop iterations when Sk is empty. Output Sb as the list of potential solutions at targeted precision δ. Output region in Sb with lowest upper bound and
corresponding best estimate at the targeted precision
Θ̂M,BB
.
N

5

(20)

Asymptotic theory in the multivariate setting

b M (see (16)) are studied theIn this section, the asymptotic properties of the exact solution Θ
N
oretically for a general multivariate OfBm [3, 10, 9]. In other words, the results encompass,
but are not restricted to, the bivariate framework of Section 2. Regarding the OfBm B W,H ,
it is assumed that: (OFBM1) B W,H is an RP -valued OfBm with Hurst parameter H, not
necessarily diagonalizable, where the eigenvalues of H satisfy 0 < <(hk ) < 1, k = 1, . . . , n ;
(OFBM2) EB W,H (t)B W,H (t)∗ , t 6= 0, is a full rank matrix (properness) ; (OFBM3) B W,H
is a time reversible stochastic process. Regarding ψ0 ∈ L1 (R), it is assumed that: (W1)
Nψ ≥ 2 ; (W2) supp(ψ0 ) is a compact interval ; (W3) supt∈R |ψ0 (t)|(1 + |t|)α < ∞ for α > 1.
The summation range in the objective function (17) is generalized to i1 , i2 = 1, . . . , P ≥ 2.
The proofs of the statements in this section can be found in Appendix B.

5.1

The asymptotic normality of the wavelet spectrum

b M draws upon the asymptotic normality of the
The asymptotic behavior of the estimator Θ
N
wavelet variance for fixed scales. Under the assumptions (OFBM1-3) and (W1-3), the latter
property can be established by an argument that is almost identical to that in [1, Theorem
3.1]. For the reader’s convenience, we reproduce the claim here.
Theorem 5.1. Suppose the assumptions (OFBM1-3) and (W1-3) hold. Let F ∈ S( P (P2+1) m, R)3
be the asymptotic covariance matrix described in [1, Proposition 3.3]. Then,

p
d
Kj (vecS S(2j ) − vecS E(2j , Θ))
→ Z,
(21)
j=j1 ,...,j2
d

as N → ∞, where j1 < . . . < j2 and m = j2 − j1 + 1 and Z = N P (P +1) ×m (0, F ), where vecS
2
defines the operator that vectorizes the upper triangular entries of a symmetric matrix:
vecS (S) = (s11 , . . . , s1P ; . . . ; sP −1,P −1 , sP −1,P ; sP,P )∗ .
3 S(n, R)

is the space of real symmetric matrices of size n × n.

9

5.2

Consistency of Θ̂M
N

Let Θ0 be the true parameter value. To prove the consistency of Θ̂M
N , the following additional
assumptions on the parametrization Θ are required:
i) The parameter space Ξ ⊆ Q0 (see (18)) is a finite-dimensional compact set and
Θ0 ∈ int Ξ;

(22)

ii) For some j ∗ = j1 , . . . , j2 ,
∗

∗

Θ 6= Θ0 ⇒ |Ei∗1 ,i∗2 (2j , Θ)| =
6 |Ei∗1 ,i∗2 (2j , Θ0 )|

(23)

for some matrix entry (i∗1 , i∗2 );
iii) ∀i1 , i2 = 1, . . . , P, j = j1 , . . . , j2 ,
Ei1 ,i2 (2j , Θ0 ) 6= 0 ;
iv) The mapping Θ 7→ E(2j , Θ)

(24)
(25)

is three times continuously differentiable on int Ξ.
Under (24), the functions log2 |Ei1 ,i2 (2j , Θ)|, i1 , i2 = 1, . . . , P , are well-defined. This
fact and Theorem 5.1 then imply that the functions log2 |Si1 ,i2 (2j )| are well-defined with
probability going to 1. In turn, condition (23) implies that the (entry-wise) absolute value
∗
of the target matrix E(2j , Θ) is (parametrically) identifiable, namely, there is an injective
j∗
function Ξ 3 Θ 7→ |E(2 , Θ)|.
The objective function CN (Θ) is a function of N , and so is S(2j ), j = j1 , . . . , j2 . Since
b N is attained (a.s.),
CN (·) is continuous and Ξ is compact, then for all N a minimum Θ
whence we can form one such sequence
bM
{Θ
N }N ∈N .

(26)

Any sequence (26) defines an M -estimator of Θ0 , e.g., [23, chapter 5]. The next theorem
shows that (26) is consistent.
Theorem 5.2. Under the assumptions of Theorem 5.1, suppose in addition that the conditions i) to iv) hold. Then, the sequence of minima (26) is consistent for Θ, namely,
bM
Θ
N → Θ0

in probability.

(27)

b M for a given N is not ensured by the conditions i) to
Remark 5.3. The uniqueness of Θ
N
iv), but it is not needed in Theorem 5.2.

5.3

Asymptotic normality of Θ̂M
N

By comparison to consistency, showing asymptotic normality will require an additional assumption, laid out next.
v) det

j2
X

X

Λi1 ,i2 (2j , Θ0 )Λi1 ,i2 (2j , Θ0 )∗



> 0,

j=j1 1≤i1 ≤i2 ≤n

(28)

where we define the score-like vector Λi1 ,i2 (2j , Θ)∗ = ∇Θ log2 |Ei1 ,i2 (2j , Θ)|.
10

Theorem 5.4. Under the assumptions of Theorem 5.1, suppose in addition that the condib M }N ∈N be a consistent sequence of minima of {CN }N ∈N , respectively.
tion v) holds. Let {Θ
N
Then,
√
d
bM
N (Θ
N → ∞, where
(29)
N − Θ0 ) → W,
d

W=

j2
X

X

j2
−1  X
2j/2
Λi1 ,i2 (2 , Θ)Λi1 ,i2 (2 , Θ)
log 2
j =1
j

j=j1 1≤i1 ≤i2 ≤P

X

∗

j

1

1≤i1 ≤i2 ≤P

Zi1 ,i2 (2j )

Λi1 ,i2 (2j , Θ) 
,
Ei1 ,i2 (2j , Θ0 )

j

where Z = (Zi1 ,i2 (2 ))j is a random vector whose distribution is obtained in the weak limit
(21).
The next result is a corollary to Theorems 5.2 and 5.4.
b M }N ∈N be a sequence of
Corollary 5.5. Under the assumptions of Theorem 5.4, let {Θ
N
M,BB
b
minima of the objective function (17). Also, let {Θ
}N ∈N be an estimator of the form
N
(20) which satisfies
C
bM − Θ
b M,BB k ≤
a.s.
(30)
kΘ
N
N
N 1/2+ε
for constants C, ε > 0. Then,
√

d
b M,BB − Θ0 ) →
N (Θ
W,
N

N → ∞,

where the random vector W is given in (29).
Remark 5.6. The condition (30) is easily satisfied in practice, since over a compact set
and at a low computational cost a Branch and Bound algorithm is guaranteed to yield a
solution which lies at a controllable distance of the true minimum.
Remark 5.7. The technical condition (28) should be satisfied in many cases of interest, as
discussed in Appendix D.

6
6.1

Estimation performance: empirical study
Numerical simulation setting

Monte Carlo experiments were performed to empirically quantify the finite size performance
of the estimator Θ̂M,BB
. To examine the influence of Θ on the estimation performance, 9
N
different values of Θ were used, obtained essentially by varying the strength of the correlation
amongst components (ρx = 0.1, 0.45 and 0.8) and of the mixing factor (no mixing, β = γ =
0 ; orthogonal mixing, β = γ = 0.5, and β = −γ = 0.5, referred to as anti-orthogonal).
Three different sample sizes (short, N = 210 , medium, N = 214 and large, N = 218 ) were
investigated. Results are reported here for (h1 , h2 ) = (0.4, 0.8). Equivalent conclusions
are obtained for other choices of (h1 , h2 ). For each set of parameters Θ, the estimation
performance were assessed by means of box plots computed from 100 independent copies.
The synthesis of OfBm is achieved by using the multivariate toolbox devised in [12, 13], cf.
www.hermir.org. The computational loads are also quantified as a percentage of the number
of iterations that would be required by a systematic greedy grid search. The wavelet analysis
was based on least asymmetric orthonormal Daubechies wavelets [17]. All available scales
were used to compute CN : j1 = 1 ≤ j ≤ j2 = log2 N − Nψ − 1. The results are reported
11

ρx = 0.45

ρx = 0.8

log2 N

log2 N

log2 N

β = −γ = 0.5

β = γ = 0.5

β=γ=0

ρx = 0.1

Figure 3: Estimation performance of h2 as function of log2 N .

for Nψ = 2; it was found that further increasing Nψ did not improve the performance. The
2
proposed Branch & Bound procedure was run with S0 = ∪∆
i=1 for δ = 50. The impact of
varying the requested precision δ was also investigated.
b M for simplicity, is compared against other
The performance of Θ̂M,BB
, re-labelled Θ
N
existing estimation procedures. The scaling exponents (h1 , h2 ) are estimated by means of
bU
the univariate wavelet based estimator for Hurst parameter, (b
hU
1 , h2 ), as described in [24]
and applied to each component independently: The univariate estimate of h1 (resp. h2 ) is
obtained by taking the minimum (resp. maximum) between the linear regression coefficients
2
of log2 S11 (2j ) and log2 S22 (2j ) versus j ∈ J fs (resp. j ∈ J cs ), with J fs = {j1 , . . . , b j1 +j
2 c}
j1 +j2
cs
(fine scales) and J = {b 2 c + 1, . . . , j2 } (coarse scales). The parameters (h1 , h2 , β) are
b W bW
also estimated using the multivariate semiparametric estimator (b
hW
1 , h2 , β ) proposed in
j
b W bW
b M , (b
[1], which relies on the multiscale eigenstructure of S(2 ). The statistics Θ
hW
1 , h2 , β )
bU
and (b
hU
1 , h2 ) are compared in Figs. 4 to 7 in yellow, blue and magenta colors, respectively.
Though the full parametrization of Biv-OfBm requires a 7-dimensional vector parameter
Θ = (h1 , h2 , ρx , σx1 , σx2 , β, γ),
for ease of exposition we focus only on the 5 most interesting parameters (h1 , h2 , ρx , β, γ).
This follows the univariate literature that focuses on the estimation of H for fBm, while
neglecting the less interesting parameter σ 2 .

6.2

Estimation performance

Estimation of the dominant scaling exponent h2 (Fig. 3). As expected, all procedures yield accurate estimates of the largest scaling exponent h2 . While all methods show
comparable performance for large sample sizes, it is interesting that b
hM
2 displays better perbU
formance with lower bias and dispersion for small sample sizes by comparison to b
hW
2 and h2 .
12

ρx = 0.45

ρx = 0.8

log2 N

log2 N

log2 N

β = −γ = 0.5

β = γ = 0.5

β=γ=0

ρx = 0.1

Figure 4: Estimation performance of h1 as function of log2 N .
The impact of the correlation ρx or the mixing parameters (β, γ) on the performance is weak.
Estimation of the non-dominant scaling exponent h1 (Fig. 4). Estimating the lowest scaling exponent h1 is intrinsically more difficult because the mixture of power laws
masks the non-dominant Hurst eigenvalue. As expected, univariate-type analysis fails to
bW
estimate correctly h1 (except when there is no mixing (β = γ = 0)). While b
hM
1 and h2
show essentially the same performance for large sample size, it is interesting to note that b
hM
1
displays a far superior performance with lower bias and dispersion for small sample sizes.
However, the bias of b
hM
1 for (β = −γ) and small correlation ρx = 0.1 are observed, showing
M
that b
h1 is more strongly affected by the conjunction of low correlation amongst components
and anti-orthogonal mixing than b
hM
2 .
Estimation of β (Fig. 5). A significant benefit of βbM consists of its being robust to small
sample sizes, when βbW is not. While the mixing parameters seem not to impact the performance of βbM , a low correlation value ρx hurts its performance. Other results not reported
here for reasons of space also show that the performance of βbM is robust to a decrease of
h2 − h1 , while that of βbW drastically deteriorates when h2 − h1 → 0.
Estimation of γ (Fig. 6). The performance of γ
bM is very satisfactory, yet it is observed
to be affected by low correlation.
Estimation of ρx (Fig. 7). The parameter ρx appears the most difficult to estimate, with
significant bias for low correlation and anti-orthogonal mixing, a result in consistency with
[3].

13

ρx = 0.45

ρx = 0.8

log2 N

log2 N

log2 N

β = −γ = 0.5

β = γ = 0.5

β=γ=0

ρx = 0.1

Figure 5: Estimation performance of β as function of log2 N .

6.3

Computational costs

The computational costs of the proposed identification algorithm described in Section 4.3.2
are reported in Fig. 8 (top plots) as a function of log2 N for each parameter setting. They
are significantly smaller than those required by a systematic greedy grid search. Fig. 8
also clearly shows that the stronger the correlation amongst components, the easier the
minimization of the functional CN , thus indicating that the cross terms play a significant
role in the identification of Biv-OfBm. Though surprising at first, the clear decrease of the
computational costs with the increase of the sample size may be interpreted as the fact that
it is obviously far easier to disentangle three different power laws when a large number of
scales 2j is available, which then requires large sample sizes. It is also worth noting that
the orthogonal mixing, which may intuitively be thought of as the easiest, appears to be the
most demanding in terms of iterations to minimize CN . Unsurprisingly, the computational
cost increases when the requested precision δ on the estimates is increased (δ → 0) Fig. 9
(left).

6.4

Sample size versus precision.

bU
Figs. 3 to 7 show that increasing the sample size N improves the performance of (b
hU
1 , h2 )
W
W
W
and (b
h1 , b
h2 , βb ), as both their median-bias and variance decrease with N . Fig. 9 (right)
b M . As long as the dispersions
indicates that the impact of N is slightly more involved for Θ
of the estimates remain above the desired precision δ (targeted independently of N ), one
observes an expected decrease in the bias and dispersion when N is increased. Because the
minimization is stopped when the prescribed accuracy δ is reached, increasing N without
decreasing δ does not bring about any performance improvement, showing that the precision
should be decreased when the sample size increases (empirically as N −1/4 ) to improve the
performance.

14

ρx = 0.45

ρx = 0.8

log2 N

log2 N

log2 N

β = −γ = 0.5

β = γ = 0.5

β=γ=0

ρx = 0.1

Figure 6: Estimation performance of γ as function of log2 N .

6.5

Asymptotic normality of Θ̂

For a thorough study of normality, Θ is restricted to (h1 , h2 ), while all other parameters are
fixed a priori and known. Averages are obtained over 1000 independent realizations of BivOfBm for each sample size. Fig. 10 (left) visually compares the empirical distributions of
b
bM
hM
1 , h2 to their best Gaussian fits. Fig. 10 (right) measures the Kullback-Leibler divergence
bM
between the empirical distributions of b
hM
1 , h2 and their best Gaussian fits, as a function of
bM
the sample size N . Fig. 10 confirms the asymptotic normality of the estimates b
hM
1 , h2 as
theoretically predicted in Section 5. It further shows that the higher the requested precision
δ, the faster the convergence to normality. Moreover, asymptotic normality is reached for
far smaller sample sizes for the largest Hurst exponent estimate b
hM
2 than for that of the
M
b
smallest one h1 .

7

Conclusion

To the best of our knowledge, this contribution proposes the first full identification procedure
for Biv-OfBm. Its originality is to formulate identification as a non-linear regression as well
as to propose a Branch & Bound procedure to provide efficient and elegant solutions to the
corresponding non-convex optimization problem.
Consistence and asymptotic normality of the estimates are shown theoretically in a
general multivariate setting. The estimation performance is assessed for finite sample sizes
by Monte Carlo simulations and found globally satisfactory for all parameters. Estimation of
parameters γ (mixing) and ρx (correlation amongst components) remain the most difficult
parameters to estimate, though no other estimation procedure has yet been proposed in
the literature. However, including estimation of γ and ρx into the non-linear regression
b W bW
formulation permits to outperform the state-of-the art method, (b
hW
1 , h2 , β ), for estimating
h1 , h2 and β, at the price of massively increasing computational costs. The proposed
Branch & Bound procedure is yet still shown to have a significantly lower computational
15

ρx = 0.45

ρx = 0.8

log2 N

log2 N

log2 N

β = −γ = 0.5

β = γ = 0.5

β=γ=0

ρx = 0.1

Figure 7: Estimation performance of ρx as function of log2 N .
cost compared to the infeasible greedy grid search strategy. The estimation performance is
satisfactory and controlled enough so that actual use of real world data can be investigated.
This estimation procedure, together with its performance assessments are paving the
road for hypothesis testing, where, e.g., testing the absence of mixing (i.e., W is diagonal) or
of correlation amongst components (i.e., ρx ≡ 0) are obviously interesting issues in practice.
Routines permitting both the identification and synthesis of OfBm will be made publicly
available at time of publication.
This work has been explored on real Internet traffic data [11].

A

Interval arithmetic

Interval arithmetic is classically used to compute lower bounds for an objective criterion CN
on a convex set R [19, 15, 20]. It relies on the explicit decomposition of CN into several
elementary functions such as sum, product, inverse, square, logarithm, exponential,. . . ,
referred to as the calculus tree. For the sake of readability, we do not detail the calculus
tree associated to the full CN , but only its first term (13), sketched in Fig. 11 (right). The
leaves of the tree, i.e., the bottom line, consist of occurrences of the variables involved in Θ.
Each node of the graph corresponds to an elementary function applied to its children. The
intervals are composed and propagated from bottom to top to obtain a bound for C within
R. The literature on interval arithmetic provides boundaries for most of the elementary
functions on intervals:
• [x, x] + [y, y] = [x + y, x + y]
• [x, x] − [y, y] = [x − y, x − y].
• [x, x] × [y, y] = [min{x × y, x × y, x × y, x × y}, max{x × y, x × y, x × y, x × y}]
• [x, x] ÷ [y, y] = [x, x] × [ y1 , y1 ], if 0 6∈ [y, y]
16

Execution time

β=γ=0
2000

β = γ = 0.5

β = −γ = 0.5

ρ x = 0.1
ρ x = 0.45
ρ x = 0.8

1500
1000
500

Percentage of iterations

0

0.2

0.1

0
10

14

18 10

14

log2 N

18 10

14

log2 N

18

log2 N

4

-4
N = 210
N = 216

2

-5

log 2 Std(b
hM
2 )

log 2 computational cost (second)

Figure 8: Computational costs as functions of log2 N . Top: Average execution time
(in seconds). Bottom: percentage of performed iterations (compared to a greedy grid search
algorithm with same accuracy).

0

-2

-6
-7
δ = 0.05
δ = 0.01
δ = 0.005
δ = 0.001

-8
-9
-10
10

-4
0

0.01

0.02

0.03

0.04

0.05

12

14

16

18

20

log 2 N

δ

Figure 9: Precision. Computational costs (left) and dispersion of the estimates, as functions of targeted precision δ

log 2 DKL ( b
hM
1)

P (b
hM
1 )

15
10
5
0
0.4

0.45

0.5

0.55

log 2 DKL ( b
hM
2)

P (b
hM
2 )

b
hM
1
20
10
0
0.75

0.8

0.85

b
hM
2

2
0
−2
−4
10

12

14
16
log 2 N

18
δ
δ
δ
δ

14
16
log 2 N

18

0

=
=
=
=

20
0.05
0.01
0.005
0.001

−2

−4
10

12

20

bM
Figure 10: Asymptotic Normality. Left: empirical distributions of b
hM
1 , h2 and best
20
−3
Gaussian fit, N = 2 , δ = 10 . Right: Kullback-Leibler divergence between the empirical
bM
distributions of b
hM
1 , h2 and their best Gaussian fits, as functions of the sample size N .
• log([x, x]) = [log(x), log(x)], if 0 6∈ [x, x]
• exp([x, x]) = [exp(x), exp(x)], if 0 6∈ [x, x]

17

0.08
0.07
0.06

ηh

0.05
0.04
0.03
0.02
0.01
0
0

0.2

0.4

0.6

0.8

1

h

Figure 11: Interval arithmetic and function ηh . Left: calculus tree associated with the
first term in Eq.(13). Right: function ηh .
The only non-elementary function involved in CN is h 7→ ηh , defined in (11) and illustrated
in Fig. 11 (right), for a least asymmetric orthogonal Daubechies wavelet ψ0 with Nψ = 2.
From the study of its monotonicity, we devise the following empirical bounding scheme:


[ηh , ηh ], if h < 0.3,
(31)
(∀[h, h] ⊆ [0, 1]) η[h,h] = [ηh , ηh ], if h ≥ 0.3,


[min(ηh , ηh ), 0.071], otherwise.

B

Proof of Theorem 5.2

Consider any sequence {ΘN }N ∈N ∈ Ξ (i.e., not necessarily composed of minima). We claim
that
P
P
CN (ΘN ) → 0 ⇒ ΘN → Θ0 .
(32)
By contradiction, assume that we can choose a subsequence {ΘN (r) }r∈N such that, with
positive probability, CN (r) (ΘN (r) ) < r−1 and kΘN (r) − Θ0 k ≥ C0 > 0. Then, the conditions (23) and (25) imply that there are indices j ∗ , i∗1 and i∗2 , a constant δ > 0 and a
∗
∗
sequence of sets Eδ,N (r) = {ω : | log2 |Ei∗1 ,i∗2 (2j , ΘN (r) )| − log2 |Ei∗1 ,i∗2 (2j , Θ0 )|| ≥ δ} such
that P (Eδ,N (r) ) ≥ C1 > 0 for some C1 > 0. So, choose ε ∈ (0, δ). By Theorem 5.1, for some
C2 ,
j2
X
1
>
r j=j

1

X
1≤i1 ≤i2 ≤P

{log2 |Si1 ,i2 (2j )|−log2 |Ei1 ,i2 (2j , Θ0 )|+log2 |Ei1 ,i2 (2j , Θ0 )|−log2 |Ei1 ,i2 (2j , ΘN (r) )|}2
n
o2
∗
∗
≥ | log2 |Ei∗1 ,i∗2 (2j , Θ0 )| − log2 |Ei∗1 ,i∗2 (2j , ΘN (r) )|| − ε ≥ C2 > 0

with non-vanishing probability (contradiction). Thus, (32) holds. Now consider the sequence
(26), and note that
P

bM
0 ≤ CN (Θ
N ) = inf CN (Θ) ≤ CN (Θ0 ) → 0,
Θ∈Ξ

by Theorem 5.1. Therefore, by (32), the limit (27) holds. 

18

(33)

C

Proof of Theorem 5.4

Rewrite
CN (Θ) =

j2
X

X

j=j1 1≤i1 ≤i2 ≤P

(fN )i1 ,i2 ,j (Θ) ∈ R,

where
(fN )i1 ,i2 ,j (Θ) = {log2 |Si1 ,i2 (2j )| − log2 |Ei1 ,i2 (2j , Θ)|}2 .

It is clear that for k ∈ N∗ ,

(k)
CN (Θ)

=

j2
X

X

(k)

(fN )i1 ,i2 ,j (Θ).

(34)

j=j1 1≤i1 ≤i2 ≤P

Fix a triple (i1 , i2 , j). By (22), (24) and (25), the first three derivatives of (fN )i1 ,i2 ,j (Θ)
with respect to Θ are well-defined in int Ξ. The first two derivatives at Θ can be expressed
as
0
(fN
)i1 ,i2 ,j (Θ) = ∇Θ (fN )i1 ,i2 ,j (Θ)∗

= −2{log2 |Si1 ,i2 (2j )| − log2 |Ei1 ,i2 (2j , Θ)|}Λi1 ,i2 (2j , Θ),

(35)
(36)

00
(fN
)i1 ,i2 ,j (Θ) = ∇Θ [∇Θ (fN )i1 ,i2 ,j (Θ)∗ ]
(37)
n
j
j
∗
j
j
= 2 Λi1 ,i2 (2 , Θ)Λi1 ,i2 (2 , Θ) − [log2 |Si1 ,i2 (2 )| − log2 |Ei1 ,i2 (2 , Θ)|]∇Θ Λi1 ,i2 (2j , Θ)}

(38)

000
)i1 ,i2 ,j (Θ) consists of sums and products of log2 |S(2j )| − log2 |E(2j , Θ)| and
Similarly, (fN
derivatives of Λi1 ,i2 (2j , Θ).
0
0
0
Rewrite CN
(Θ) = {(CN
)l (Θ)}l=1,...,dim Ξ . By a second order Taylor expansion of (CN
)l (Θ)
with Lagrange remainder,
0
b N ) − (C 0 )l (Θ0 ) = ∇Θ (C 0 )l (Θ0 )(Θ
b N − Θ0 )
R 3 (CN
)l (Θ
N
N
0
∗
e
b N − Θ0 ),
b N − Θ0 )∗ ∇Θ [∇Θ (CN )l ((ΘN )l ) ] (Θ
+ (Θ
2!

(39)

e N )l , l = 1, . . . , dim Ξ, is a parameter value lying in a segment between
where each entry (Θ
b
Θ0 and ΘN . Thus,
0 b
0
00
b N − Θ0 )
Rdim Ξ 3 CN
(ΘN ) − CN
(Θ0 ) = CN
(Θ0 )(Θ


0
∗ 
e
b N − Θ0 )∗ ∇Θ [∇Θ (CN )l ((ΘN )l ) ]
b N − Θ0 ),
+ (Θ
(Θ
(40)
2!
l=1,...,dim Ξ
b N − Θ0 )∗ ∇Θ [∇Θ (C 0 )l ((Θ
e N )l )∗ ], l = 1, . . . , dim Ξ, is a row vector. By
where each entry (Θ
N
0 b
b
(22), ΘN ∈ int Ξ for large N with probability going to 1. Thus, CN
(ΘN ) = 0. Solving (40)
b N − Θ0 yields
for Θ
√

n

o−1 √
0
∗ 
e
00
0
b N −Θ0 ) = − CN
b N −Θ0 )∗ ∇Θ [∇Θ CN ((ΘN )l ) ]
N (Θ
(Θ0 )+ (Θ
N CN
(Θ0 ).
2!
l=1,...,dim Ξ
(41)
19

b N for Θ0 and the expression for (f 000 )i ,i ,j (Θ), we have
Under (25), by the consistency of Θ
N 1 2
for l = 1, . . . , dim Ξ,
P
0
0
e N )l )∗ ] →
∇Θ [∇Θ CN
((Θ
∇Θ [∇Θ CN
(Θ0 )∗ ].
(42)
The invertibility of the matrix between braces in (41) is ensured with probability going to
1 by the condition (28) and Theorem 5.1, since the expression (37) entails
j2
X

P

00
CN
(Θ0 ) → 2

X

Λi1 ,i2 (2j , Θ0 )Λi1 ,i2 (2j , Θ0 )∗ .

(43)

j=j1 1≤i1 ≤i2 ≤P

Let k · kl1 be the entry-wise l1 matrix norm. By the relations (34) and (35), as well as a
first order Taylor
of log2 | · | around Ei1 ,i2 (2j , Θ0 ) under the condition (24), we
√ expansion
0
(Θ0 ) as
can reexpress N CN
−2

j2
X
j=j1

p X Si1 ,i2 (2j ) − Ei1 ,i2 (2j , Θ0 )
2j/2 Kj
Λi1 ,i2 (2j , Θ0 )
(log 2)Ei1 ,i2 (2j , Θ0 )
1≤i1 ≤i2 ≤P

j2

X
p
Kj kS(2j ) − E(2j , Θ0 )kl1 ∈ Rdim Ξ .
+o

(44)

j=j1

Then, the weak limit (29) is a consequence of (41), (42), (43), (44), Theorem 5.1 and the
Cramér-Wold device. 

D

Discussion of remark 5.7

The condition (30) amounts to requiring the full rank of a sum of rank 1 terms. To fix ideas,
consider a sum of the form vv ∗ + ww∗ , where v, w ∈ R2 \{0}. Then, the sum has deficient
rank 1 if and only if v and w are collinear. Indeed, assume that there is a vector u 6= 0 such
that u∗ {vv ∗ + ww∗ }u = 0. Thus, u∗ vv ∗ u = 0 = u∗ ww∗ u, whence collinearity follows.
So, for notational simplicity, we assume that we can rewrite the wavelet spectrum as
|Ei1 ,i2 (2j , Θ)| = ai1 ,i2 2j2h1 + bi1 ,i2 2j(h1 +h2 ) + ci1 ,i2 2j2h2 ,
where i1 , i2 = 1, . . . , P (c.f. [1], Lemma 4.2). Further assume that the only parameters to
be estimated are h1 < h2 . Then, for a fixed j and a pair of indices (i1 , i2 ),
∂
1 ai1 ,i2 log(22j )2j2h1 + bi1 ,i2 log(2j )2j(h1 +h2 )
log2 |Ei1 ,i2 (2j , Θ)| =
,
∂h1
log 2 ai1 ,i2 2j2h1 + bi1 ,i2 2j(h1 +h2 ) + ci1 ,i2 2j2h2
1 bi1 ,i2 log(22j )2j(h1 +h2 ) + ci1 ,i2 log(2j )2j2h2
∂
log2 |Ei1 ,i2 (2j , Θ)| =
.
∂h2
log 2 ai1 ,i2 2j2h1 + bi1 ,i2 2j(h1 +h2 ) + ci1 ,i2 2j2h2
This suggests that for at least two triplets (i1 , i2 , j), the vectors
 ∂
∗
∂
log2 |Ei1 ,i2 (2j , Θ)|,
log2 |Ei1 ,i2 (2j , Θ)|
∂h1
∂h2
will not be collinear for most parametrizations in practice.
20

References
[1] P. Abry and G. Didier. Wavelet estimation for operator fractional Brownian motion.
Bernoulli, page to appear, 2015.
[2] S. Achard, D. Bassett, A. Meyer-Lindenberg, and E. Bullmore. Fractal connectivity of
long-memory networks. Phys. Rev. E, 77(3):036104, 2008.
[3] P.-O. Amblard and J.-F. Coeurjolly. Identification of the Multivariate Fractional Brownian Motion. IEEE Trans. Signal Process., 59(11):5152–5168, Nov. 2011.
[4] P.-O. Amblard, J.-F. Coeurjolly, F. Lavancier, and A. Philippe. 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, 2012.
[5] J.-M. Bardet, G. Lang, G. Oppenheim, A. Philippe, S. Stoev, and M. Taqqu.
Semi-parametric estimation of the long-range dependence parameter: a survey. In
P. Doukhan, G. Oppenheim, and M. S. Taqqu, editors, Theory and applications of
Long-range dependence, pages 557–577, Boston, 2003. Birkhäuser.
[6] P. Ciuciu, P. Abry, and B. J. He. Interplay between functional connectivity and scalefree dynamics in intrinsic fMRI networks. NeuroImage, 95:248–263, 2014.
[7] J.-F. Coeurjolly, P.-O. Amblard, and S. Achard. Wavelet analysis of the multivariate
fractional Brownian motion. ESAIM: Probability and Statistics, 17:592–604, Aug. 2013.
[8] G. Didier, H. Helgason, and P. Abry. Demixing multivariate-operator self-similar processes. In Proc. Int. Conf. Acoust., Speech Signal Process., Brisbane, Australia, 2015.
[9] G. Didier and V. Pipiras. Integral representations and properties of operator fractional
Brownian motions. Bernoulli, 17(1):1–33, 2011.
[10] G. Didier and V. Pipiras. Exponents, symmetry groups and classification of operator
fractional Brownian motions. Journal of Theoretical Probability, 25(2):353–395, 2012.
[11] J. Frecon, R. Fontugne, G. Didier, N. Pustelnik, K. Fukuda, and P. Abry. Non-linear
regression for bivariate self-similarity identification - application to anomaly detection
in Internet traffic based on a joint scaling analysis of packet and byte counts. In Proc.
Int. Conf. Acoust., Speech Signal Process., Shanghai, China, March 2016.
[12] H. Helgason, V. Pipiras, and P. Abry. Fast and exact synthesis of stationary multivariate
Gaussian time series using circulant embedding. Signal Process., 91(5):1123 – 1133,
2011.
[13] H. Helgason, V. Pipiras, and P. Abry. Synthesis of multivariate stationary series with
prescribed marginal distributions and covariance using circulant matrix embedding.
Signal Process., 91(8):1741 – 1758, 2011.
[14] K. Ichida and Y. Fujii. An interval arithmetic method for global optimization. Computing, 23(1):85–97, 1979.
[15] L. Jaulin, M. Kieffer, and O. Didrit. Applied interval analysis : with examples in
parameter and state estimation, robust control and robotics. Springer, London, 2001.
21

[16] R. B. Kearfott. An interval branch and bound algorithm for bound constrained optimization problems. J. Global Optim., 2(3):259–280, 1992.
[17] S. Mallat. A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way.
Academic Press, 3rd edition, 2008.
[18] R. Moore, E. Hansen, and A. Leclerc. Recent advances in global optimization. chapter
Rigorous Methods for Global Optimization, pages 321–342. Princeton University Press,
Princeton, NJ, USA, 1992.
[19] R. E. Moore. Interval analysis. Prentice-Hall Inc., Englewood Cliffs, N.J., 1966.
[20] R. E. Moore, R. B. Kearfott, and M. J. Cloud. Introduction to Interval Analysis. Society
for Industrial and Applied Mathematics, 2009.
[21] H. Ratschek and J. Rokne. New Computer Methods for Global Optimization. Halsted
Press, New York, NY, USA, 1988.
[22] G. Samorodnitsky and M. Taqqu. Stable non-Gaussian random processes. Chapman
and Hall, New York, 1994.
[23] A. Van der Vaart. Asymptotic Statistics, volume 3. Cambridge University Press, 2000.
[24] D. Veitch and P. Abry. A wavelet-based joint estimator of the parameters of long-range
dependence. IEEE Trans. Inform. Theory, 45(3):878–897, 1999.
[25] H. Wendt, P. Abry, and S. Jaffard. Bootstrap for empirical multifractal analysis. IEEE
Signal Process. Mag., 24(4):38–48, 2007.

22

