Multivariate Hadamard self-similarity: testing fractal connectivity
Herwig Wendt, Gustavo Didier, Sébastien Combrexelle, Patrice Abry
MMULTIVARIATE HADAMARD SELF-SIMILARITY:TESTING FRACTAL CONNECTIVITY
HERWIG WENDT , GUSTAVO DIDIER , S´EBASTIEN COMBREXELLE , PATRICE ABRY Abstract.
While scale invariance is commonly observed in each component of real worldmultivariate signals, it is also often the case that the inter-component correlation structureis not fractally connected, i.e., its scaling behavior is not determined by that of the indi-vidual components. To model this situation in a versatile manner, we introduce a class ofmultivariate Gaussian stochastic processes called Hadamard fractional Brownian motion(HfBm). Its theoretical study sheds light on the issues raised by the joint requirementof entry-wise scaling and departures from fractal connectivity. An asymptotically normalwavelet-based estimator for its scaling parameter, called the Hurst matrix, is proposed, aswell as asymptotically valid confidence intervals. The latter are accompanied by originalfinite sample procedures for computing confidence intervals and testing fractal connectiv-ity from one single and finite size observation. Monte Carlo simulation studies are used toassess the estimation performance as a function of the (finite) sample size, and to quantifythe impact of omitting wavelet cross-correlation terms. The simulation studies are shownto validate the use of approximate confidence intervals, together with the significance leveland power of the fractal connectivity test. The test performance and properties are furtherstudied as functions of the HfBm parameters.
Keywords.
Multivariate self-similarity, Hadamard self-similarity, fractal connectivity,wavelet analysis, confidence intervals, hypothesis testing Introduction
Scale Invariance.
The relevance of the paradigm of scale invariance is evidencedby its successful use, over the last few decades, in the analysis of the dynamics in dataobtained from a rather diverse spectrum of real world applications. The latter range fromnatural phenomena – physics (hydrodynamic turbulence [40], out-of-equilibrium physics),geophysics (rainfalls), biology (body rhythms [37], heart rate [10, 31], neurosciences andgenomics [16, 17, 23, 24]) – to human activity – Internet traffic [1, 7], finance [39], urbangrowth and art investigation [6, 8, 28].In essence, scale invariance – also called scaling, or scale-free dynamics – implies that thephenomenical or phenomenological dynamics are driven by a large continuum of equallyimportant time scales, rather than by a small number of characteristic scales. Thus, the IRIT-ENSEEIHT, CNRS (UMR 5505), Universit´e de Toulouse, France, [email protected] Mathematics Department, Tulane University, New Orleans, USA, [email protected] Universit´e de Lyon, ENS de Lyon, Universit´e Claude Bernard, CNRS, Laboratoire dePhysique, F-69342 Lyon, France, [email protected]
P.A. and H.W. are supported in part by the French ANR grant MultiFracs. G.D. was supported in partby the ARO grant W911NF-14-1-0475. G.D. gratefully acknowledges the support of ENS de Lyon for hisvisits, and thanks Alexandre Belloni for the insightful mathematical discussions. a r X i v : . [ m a t h . S T ] J a n MULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY investigation’s focus is on identifying a relation amongst relevant scales rather than pickingout characteristic scales.Historically, self-similarity was one of the first proposed mathematical frameworks forthe modeling of scale invariance (e.g., [44]). A random system is called self-similar whendilated copies of a single signal X are statistically indistinguishable, namely,(1.1) { X ( t ) } t ∈ R fdd = { a H X ( t/a ) } t ∈ R , ∀ a > , where fdd = stands for the equality of finite-dimensional distributions. An example of a sto-chastic process that satisfies the property (1.1) is fractional Brownian motion (fBm). Indeed,the latter is the only self-similar, Gaussian, stationary increment process, and it is the mostwidely used scaling model for real-world signals [46].Starting from (1.1), the key parameter for quantifying scale-free dynamics is the scaling,or Hurst, exponent 0 < H <
1. The estimation of H is the central task in scaling analysis,and it has received considerable effort and attention in the last three decades (see [14] fora review). The present contribution is about wavelet-based estimation [22, 47]. It relies onthe key scaling property(1.2) 1 T (cid:88) t T X ( a, t ) (cid:39) Ca α , α := 2 H, where T X ( a, t ) is the wavelet coefficient of an underlying self-similar stochastic process and T is the number of available coefficients. In other words, the sample wavelet variance of thestochastic process behaves like a power law with respect to the scale a .1.2. Multivariate scaling.
In many modern fields of application such as Internet trafficand neurology, data is collected in the form of multivariate time series. Univariate-likeanalysis in the spirit of (1.2) – i.e., independently on each component – does not accountfor the information stemming from correlations across components. The classical fBmparametric family, for example, provides at best a model for component-wise scaling, andthus cannot be used as the foundation for a multivariate modeling paradigm.To model self-similarity in a multivariate setting, a natural extension of fBm, calledOperator fractional Brownian motion (OfBm), was recently defined and studied (see [11,18, 21]). An OfBm X satisfies the m -variate self-similarity relation(1.3) { X ( t ) } t ∈ R fdd = { a H X ( t/a ) } t ∈ R , ∀ a > , where the scaling exponent is a m × m matrix H , and a H stands for the matrix exponential (cid:80) ∞ k =0 ( H log a ) k /k !. Likewise, the wavelet spectrum of each individual component is not asingle power law as in (1.2); instead, it behaves like a mixture of distinct univariate powerlaws. In its most general form, OfBm remains scarcely used in applications; recent effortshave tackled many difficulties that arise in the identification of its parameters [3, 20].1.3. Entry-wise multivariate scaling.
We call an OfBm entry-wise scaling when theHurst parameter is simply a diagonal matrix H = diag( H , . . . , H m ). This instance ofOfBm has been used in many applications (e.g., [9, 16]) and its estimation is thoroughlystudied in [11]. Since H is diagonal, the relation (1.3) takes the form(1.4) { X ( t ) , . . . , X m ( t ) } t ∈ R fdd = { a H X ( t/a ) , . . . , a H m X m ( t/a ) } t ∈ R , ∀ a > , ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 3 which is reminiscent of the univariate case. This implies that the extension of (1.2) to allauto- and cross-components of m -variate data can be written as(1.5) 1 T (cid:88) t T X q ( a, t ) T X q ( a, t ) (cid:39) Ca α q q , α q q := H q + H q , q , q = 1 , . . . , m, where T is as in (1.2).1.4. Fractal connectivity.
Yet, entry-wise scaling OfBm is a restrictive model since thecross-scaling exponents α q q , q (cid:54) = q , are determined by the auto-scaling exponents α q q and α q q , i.e.,(1.6) α q q = H q + H q = ( α q q + α q q ) / . In this situation, called fractal connectivity [9, 16, 49], no additional scaling information canbe extracted from the analysis of cross-components . However, in real world applications,cross-components are expected to contain information on the dynamics underlying the data,e.g., cross-correlation functions [32, 33]. As an example, recent investigation of multivariatebrain dynamics in [16] produced evidence of departures from fractal connectivity, notablyfor subjects carrying out prescribed tasks. This means that there is a clear need for moreversatile models than entry-wise scaling OfBm (see also Remark 2.2 below). The covariancestructure of the new model should satisfy the following two requirements:(1) all auto- and cross-components are (approximately) self-similar;(2) departures from fractal connectivity are allowed, i.e., the exponents of the cross-components are not necessarily determined by the exponents of the correspondingauto-components.Hereinafter, a departure from fractal connectivity (1.6) on a given covariance structureentry ( q , q ) will be quantified by means of the parameter(1.7) δ q q = α q q + α q q − α q q ≥ , q , q = 1 , . . . , m, where nonnegativeness is a consequence of the Cauchy-Schwarz inequality (see (2.9) below).It is clear that δ q q = 0 when q = q .1.5. Goals, contributions and outline.
Our contribution comprises four main compo-nents. First, we propose a new class of multivariate Gaussian stochastic processes, called
Hadamard fractional Brownian motion (HfBm), that combines scale-free dynamics and po-tential departures from fractal connectivity. Moreover, we provide a precise discussion ofthe issues entailed by the presence of these two properties (Section 2). Second, we study themultivariate discrete wavelet transform (DWT) of HfBm, define wavelet-based estimatorsfor the scaling exponents α q q and the fractal connectivity parameter δ q q , mathematicallyestablish their asymptotic performance (i.e., asymptotic normality and covariance struc-ture), and computationally quantify finite sample size effects (Section 3). Third, startingfrom a single sample path, we construct approximate confidence intervals for the proposedestimators. The procedure is backed up by the mathematical identification of the approx-imation orders as a function of the sample path size and the limiting coarse scale. This isfurther investigated by means of Monte Carlo simulations, as well as by means of a study ofthe ubiquitous issue of the impact of (partially) neglecting the correlation amongst waveletcoefficients (Section 4). Beyond being of interest in multivariate modeling, the study shedsmore light on the same issue for the univariate case. Fourth, we devise an efficient test for MULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY the presence of fractal connectivity from a single sample path. In addition, we assess thefinite sample performance of the test in terms of targeted risk by means of Monte Carlosimulations (Section 5). Finally, routines for the synthesis of HfBm, as well as for estima-tion, computation of confidence intervals and testing will be made publicly available at timeof publication. All proofs can be found in the Appendix.2.
Hadamard fractional Brownian motion
For Hadamard fractional Brownian motion, defined next, the fractal connectivity relation(1.6) does not necessarily hold.
Definition 2.1.
A Hadamard fractional Brownian motion B H = { B H ( t ) } t ∈ R (HfBm) is aproper, Gaussian (stationary increment) process whose second moments can be written as(2.1) E [ B H ( s ) B H ( t ) ∗ ] = (cid:90) R (cid:16) e isx − ix (cid:17)(cid:16) e − itx − − ix (cid:17) f H ( x ) dx, s, t ∈ R , For 0 < h min ≤ h max <
1, the matrix exponent H = (cid:16) h q q (cid:17) q ,q =1 ,...,m satisfies the condi-tions(2.2) h q q ∈ [ h min , h max ] , q , q = 1 , . . . , m. The matrix-valued function f is a spectral density of the form(2.3) f H ( x ) q q = (cid:16) ρ q q σ q σ q | x | − h q q − / (cid:17) g q q ( x ) , q , q = 1 , . . . , m, i.e., the Hadamard scaling parameters are given by(2.4) α q q = 2 h q q , q , q = 1 , . . . , m, where ρ q q ∈ [ − , σ q , σ q ∈ R + . The real-valued functions g q q ∈ C ( R ) satisfy(2.5) max l =0 , , sup x ∈ R (cid:12)(cid:12)(cid:12) d l dx l g q q ( x ) (cid:12)(cid:12)(cid:12) ≤ C, (2.6) (cid:12)(cid:12)(cid:12) d l dx l ( g q q ( x ) − (cid:12)(cid:12)(cid:12) ≤ C (cid:48) | x | (cid:36) − l , x ∈ ( − ε , ε ) , l = 0 , , , for constants C, C (cid:48) , ε >
0, where(2.7) 2 h max < (cid:36) ≤ h min ) . Example 2.1.
An HfBm with parameters h q q = ( h q q + h q q ) / g q q ( x ) ≡ q , q = 1 , . . . , m , is an entry-wise scaling OfBm with diagonal Hurstmatrix H = diag( h , . . . , h m ) (see [11, 18, 21]).By the known properties of spectral densities, f H ( x ) = ( f H ( x ) q q ) q ,q =1 ,...,m ∈S > ( m, R ) (symmetric positive semidefinite) a.e. and satisfies(2.8) | f H ( x ) q q | ≤ (cid:113) f H ( x ) q q (cid:113) f H ( x ) q q dx -a.e.( [15], p.436). The relation (2.8) further implies that(2.9) h q q ≤ h q q + h q q , q , q = 1 , . . . , m. ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 5
Whenever convenient we also write h qq = h q , α qq = α q , q = 1 , . . . , m. The name “Hadamard” comes from Hadamard (entry-wise) matrix products. If onerewrites HfBm componentwise as B H ( t ) = ( B H, ( t ) , . . . , B H,m ( t )) ∗ , then the conditions(2.1), (2.2) and (2.6) yield the asymptotic equivalence E [ B H,q ( cs ) B H,q ( ct )] ∼ ς {| cs | h q q + | ct | h q q − | c ( s − t ) | h q q } , c → ∞ ,q , q = 1 , . . . , m , for some ς ∈ R . In other words, over large scales, the covariance betweeneach pair of entries of an HfBm approaches that of a univariate fBm, up to a change ofsign (see also Proposition 3.1, ( iii )). In this sense, for large c , an HfBm behaves like itsideal counterpart B H, ideal = { B H, ideal ( t ) } t ∈ R , defined as a generally non-existent stochasticprocess satisfying the also ideal Hadamard (entry-wise) self-similarity relation(2.10) E [ B H, ideal ( cs ) B H, ideal ( ct ) ∗ ] = c ◦ H ◦ E [ B H, ideal ( s ) B H, ideal ( t ) ∗ ] , c > , where ◦ denotes the Hadamard (entry-wise) matrix product and c ◦ H := (cid:16) c h q q (cid:17) q ,q =1 ,...,m .The process B H, ideal can be viewed as a heuristic tool for developing intuition on multivariateself-similarity. Mathematically, though, it can only exist in fractally connected instances,the reason being that distinct (spectral) power laws cross over. Indeed, we must have α q q ≤ α q q + α q q for x close to 0 and α q q ≥ α q q + α q q for large | x | , whence α q q ≡ α q q + α q q . This shows that HfBm is a perturbation of its virtual counterpart, where theregularization functions g q q ( x ) in (2.3) introduce high-frequency corrections. Example 2.2.
An illustrative subclass of HfBm is obtained by setting g qq ( x ) ≡ σ q , q = 1 , . . . , m , and g q q ( x ) = e − x , q (cid:54) = q . Note that g q q ( · ) satisfies (2.7) with (cid:36) = 2.In this case, the expression for the main diagonal spectral entry of an HfBm is identical tothat of an ideal-HfBM, and the difference lies on the off-diagonal entries:(2.11) ρ q q | x | − α q q e − x ≤ ρ q q ρ q q | x | − ( α q q + α q q ) dx -a.e.In this case, each individual entry { B H ( t ) q } t ∈ R , q = 1 , . . . , m , of HfBm B H is by itself afBm with Hurst parameter 0 < h q <
1. In particular, { B H ( ct ) q } t ∈ R fdd = { c h q B H ( t ) q } t ∈ R , c > , q = 1 , . . . , m. However, it is generally not true that { B H ( ct ) } t ∈ R fdd = { diag( c h , . . . , c h q , . . . , c h m ) B H ( t ) } t ∈ R , c > . Otherwise, B H would necessarily be fractally connected. Remark 2.1.
When simulating HfBm via Circulant Matrix Embedding, one verifies thatregularization is rarely necessary for ensuring the positive definiteness of the covariancematrix for finite sample sizes. In other words, ideal-HfBM is also a useful approximation inpractice.
Remark 2.2.
Let Y H ( t ) = B H ( t ) − B H ( t − { B H ( t ) } t ∈ R is an HfBm. Then, { Y H ( t ) } t ∈ Z is a (discrete time) stationary process with spectral density f Y H ( x ) = 2(1 − cos( x )) ∞ (cid:88) k = −∞ f H ( x + 2 πk ) | x + 2 πk | , x ∈ ( − π, π ] , MULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY where f H is the HfBm spectral density (2.1). Under (2.5), (2.6) and (2.7), we obtain theentry-wise limiting behavior f Y H ( x ) q q ( x ) ∼ ρ q q σ q σ q | x | − h q q − / , x → + , q , q = 1 , . . . , m. In particular, { Y H ( t ) } t ∈ Z does not fall within the scope of the usual definitions of multi-variate scaling behavior or long range dependence [30, 36, 42, 43], which are restricted tofractally connected processes.3. Wavelet-based analysis of HfBm
In this section, we construct the wavelet analysis and estimation for HfBm. Due to themathematical convenience of the notion of Hadamard (approximate) self-similarity, most ofthe properties of wavelet-based constructs resemble their univariate analogues.3.1.
Multivariate discrete wavelet transform.
Throughout the rest of the paper, wewill make the following assumptions on the underlying wavelet basis. Such assumptions willbe omitted in the statements.
Assumption ( W ψ ∈ L ( R ) is a wavelet function, namely,(3.1) (cid:90) R ψ ( t ) dt = 1 , (cid:90) R t q ψ ( t ) dt = 0 , q = 0 , , . . . , N ψ − , N ψ ≥ . Assumption ( W ) : the functions ϕ (a bounded scaling function) and ψ correspond to(3.2) a MRA of L ( R ), and supp( ϕ ) and supp( ψ ) are compact intervals. Assumption ( W β > (cid:98) ψ ∈ C ( R )and(3.3) sup x ∈ R | (cid:98) ψ ( x ) | (1 + | x | ) β < ∞ . Under (3.1)–(3.3), ψ is continuous, (cid:98) ψ ( x ) is everywhere differentiable and(3.4) (cid:98) ψ ( l ) (0) = 0 , l = 0 , . . . , N ψ − Definition 3.1.
Let B H = { B H ( t ) } t ∈ R ∈ R m be an HfBm. For a scale parameter j ∈ N and a shift parameter k ∈ Z , its ( L -normalized) wavelet transform is defined by(3.5) D (2 j , k ) = 2 − j/ (cid:90) R − j/ ψ (2 − j t − k ) B H ( t ) dt =: (cid:16) d q ( j, k ) (cid:17) q =1 ,...,m . Under (3.1)-(3.3) and the continuity of the covariance function (2.1), the wavelet trans-form (3.5) is well-defined in the mean-square sense and E D (2 j , k ) = 0, k ∈ Z , j ∈ N (see [19],p. 86). ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 7
Multivariate wavelet spectrum.
Fix j ≤ j , j , j ∈ N . Because of the approxi-mate nature of Hadamard self-similarity, analysis and estimation must be considered in thecoarse scale limit, by means of a sequence of dyadic numbers { a ( n ) } n ∈ N satisfying(3.6) 1 ≤ na ( n )2 j ≤ na ( n )2 j ≤ n, a ( n ) h max − h min )+1 n + na ( n ) (cid:36) → , n → ∞ . Example 3.1.
An example of a scaling sequence satisfying (3.6) for large enough n is a ( n ) := 2 (cid:98) n η h max − h min)+1 (cid:99) , h max − h min ) + 11 + 2 (cid:36) < η < . In other words, a wide parameter range [ h min , h max ] or a low regularity parameter value (cid:36) implies that a ( n ) must grow slowly by comparison to n . Remark 3.1.
For a fixed octave range { j , . . . , j } associated with an initial scaling fac-tor value a ( n ) = 1 and sample size n , define the scale range { j ( n ) , . . . , j ( n ) } = { a ( n )2 j , . . . , a ( n )2 j } for a general sample size n (where a ( n ) is assumed dyadic). Then,under (3.6), for every n the range of useful octaves is constant and given by j ( n ) − j ( n ) = j − j , where the new octaves are j l ( n ) = log a ( n ) + j l ∈ N , l = 1 , Definition 3.2.
Let B H = { B H ( t ) } t ∈ R be an HfBm. Let { D ( a ( n )2 j , k ) } k =1 ,...,n j , j = j ,...,j be its wavelet coefficients, where n a,j = na ( n )2 j is the number of wavelet coefficients D ( a ( n )2 j , · ) at scale a ( n )2 j , of a total of n (wavelet)data points. The wavelet variance and the sample wavelet variance, respectively, at scale a ( n )2 j are defined by(3.7) E (cid:2) S n ( a ( n )2 j ) (cid:3) , S n ( a ( n )2 j ) = n a,j (cid:88) k =1 D ( a ( n )2 j , k ) D ( a ( n )2 j , k ) ∗ n a,j , j = j , . . . , j . Let(3.8) (cid:37) ( q q ) ( a ( n )2 j ) = E (cid:2) d q ( a ( n )2 j , d q ( a ( n )2 j , (cid:3) . The standardized counterparts of (3.7) are(3.9) E (cid:2) W n ( a ( n )2 j ) (cid:3) , W n ( a ( n )2 j ) = (cid:37) ( a ( n )2 j ) ◦− ◦ n a,j (cid:88) k =1 D ( a ( n )2 j , k ) D ( a ( n )2 j , k ) ∗ n a,j ,j = j , . . . , j , where (cid:37) ( a ( n )2 j ) = (cid:16) (cid:37) ( q q ) ( a ( n )2 j ) (cid:17) q ,q =1 ,...,m . Entry-wise, for 1 ≤ q , q ≤ m ,(3.10) S ( q q ) n ( a ( n )2 j ) = n a,j (cid:88) k =1 d q ( a ( n )2 j , k ) d q ( a ( n )2 j , k ) n a,j , W ( q q ) n ( a ( n )2 j ) = S ( q q ) n ( a ( n )2 j ) (cid:37) ( q q ) ( a ( n )2 j ) . Proposition 3.1, stated next, provides basic results on the moments and asymptotic dis-tribution, in the coarse scale limit a ( n )2 j → ∞ , of the wavelet transform (3.5) and variance(3.9). In particular, in regard to the limits in distribution, the vector of random ma-trices { S n ( a ( n )2 j ) } j = j ,...,j can be intuitively interpreted as an asymptotically unbiased MULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY and Gaussian estimator of its population counterpart { E (cid:2) S n ( a ( n )2 j , k ) (cid:3) } j = j ,...,j . The cel-ebrated decorrelation property of the wavelet transform (e.g., [13], Proposition II.2) liesbehind the Gaussian limits in the proposition, as well as of the fact that the random ma-trices S n ( a ( n )2 j ) have weak inter-component, intra-scale and inter-scale correlations. Also,each entry S ( q q ) n ( a ( n )2 j ), q , q = 1 , . . . , m , displays asymptotic power law scaling. InSection 3.3 below, these properties are used to define estimators of the scaling exponentsand to analytically establish their asymptotic performance.In the statement of Proposition 3.1, we make use of the operator(3.11) vec S S = ( s , s , . . . , s m ; s , . . . , s m ; . . . ; s m − ,m − , s m − ,m ; s m,m ) ∗ . In other words, vec S · vectorizes the upper triangular entries of the symmetric matrix S . Proposition 3.1.
Let B H = { B H ( t ) } t ∈ R be an HfBm. Consider j, j (cid:48) ∈ N and k, k (cid:48) ∈ Z .Then, ( i ) for fixed j , j (cid:48) , k , k (cid:48) , we can write (3.12) E (cid:104) D ( a ( n )2 j , k ) D ( a ( n )2 j (cid:48) , k (cid:48) ) ∗ (cid:105) = (cid:16) Ξ jj (cid:48) ,aq q ( a ( n )(2 j k − j (cid:48) k (cid:48) )) (cid:17) q ,q =1 ,...,m for some matrix-valued function Ξ jj (cid:48) ,a ( · ) that depends on a ( n ) . In particular, forfixed scales j , j (cid:48) , { D ( a ( n )2 j , k ) } k ∈ Z is a stationary stochastic process; ( ii ) for q , q = 1 , . . . , m , (3.13) (cid:37) ( q q ) ( a ( n )2 j ) = Ξ jj,aq q (0) (cid:54) = 0 , j ∈ N , i.e., (3.9) is well-defined; ( iii ) for Ξ jj (cid:48) ,a ( · ) as in (3.12) , q , q = 1 , . . . , m , and z ∈ Z , as n → ∞ , (3.14) Ξ jj (cid:48) ,aq q ( a ( n ) z ) a ( n ) h q q → Φ jj (cid:48) q q ( z ) = ρ q q σ q σ q (cid:90) R e izx | x | − (2 h q q +1) (cid:98) ψ (2 j x ) (cid:98) ψ (2 j (cid:48) x ) dx. In particular, the wavelet spectrum E S n ( a ( n )2 j ) (see (3.7) ) can be approximated bythat of an ideal-HfBM (see (2.10) ) in the sense that (3.15) a ( n ) ◦ − H ◦ E (cid:2) S n ( a ( n )2 j ) (cid:3) → E (cid:2) S ideal (2 j ) (cid:3) = (cid:16) Φ jjq q (0) (cid:17) q ,q =1 ,...,m ;( iv ) for Ξ jj (cid:48) as in (3.12) and q , q = 1 , . . . , m , as n → ∞ , √ n a,j √ n a,j (cid:48) n a,j (cid:88) k =1 n a,j (cid:48) (cid:88) k (cid:48) =1 Ξ jj (cid:48) ,aq q ( a ( n )(2 j k − j (cid:48) k (cid:48) )) a ( n ) h q q (3.16) → − ( j + j (cid:48) )2 gcd(2 j , j (cid:48) ) ∞ (cid:88) z = −∞ Φ jj (cid:48) q q ( z gcd(2 j , j (cid:48) )); where Φ jj (cid:48) q q ( · ) is given by (3.14) ; (v) for ≤ q ≤ q ≤ m , ≤ q ≤ q ≤ m , as n → ∞ , (3.17) √ n a,j a ( n ) δ q q √ n a,j (cid:48) a ( n ) δ q q Cov (cid:104) W ( q q ) n ( a ( n )2 j ) , W ( q q ) n ( a ( n )2 j (cid:48) ) (cid:105) → G jj (cid:48) ( q , q , q , q ) , ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 9 for δ ·· as in (1.7) , where Φ jjq q (0)Φ j (cid:48) j (cid:48) q q (0)2 − ( j + j (cid:48) )2 G jj (cid:48) ( q , q , q , q )= φ jj (cid:48) ( q , q , q , q ) , if δ q q = 0 = δ q q and ( δ q q > δ q q > φ jj (cid:48) ( q , q , q , q ) , if ( δ q q > δ q q > δ q q = 0 = δ q q ; φ jj (cid:48) ( q , q , q , q ) + φ jj (cid:48) ( q , q , q , q ) , if δ q q = δ q q = δ q q = δ q q = 0;0 , otherwise , and (3.18) φ jj (cid:48) ( q , q , q , q ) := gcd(2 j , j (cid:48) ) ∞ (cid:88) z = −∞ Φ jj (cid:48) q q ( z gcd(2 j , j (cid:48) ))Φ jj (cid:48) q q ( z gcd(2 j , j (cid:48) ));(vi) as n → ∞ , (3.19) (cid:16) vec S (cid:104)(cid:16) √ n a,j a ( n ) δ q q (cid:17) q ,q =1 ,...,m ◦ ( W n ( a ( n )2 j ) − ) (cid:105)(cid:17) j = j ,...,j d → N m ( m +1)2 × J (0 , G ) , where (3.20) J = j − j + 1 , and is a vector of ones. Each entry of the asymptotic covariance matrix G in (3.19) is given by the terms G jj (cid:48) ( q , q , q , q ) in (3.17) for appropriate values of j , j (cid:48) , q , q , q and q . Remark 3.2.
Note that, up to a change of sign, each entry Φ jj (cid:48) q q ( z ) on the right-hand sideof (3.14) corresponds to the covariance between the wavelet coefficients at octaves j and j (cid:48) of a fBm with parameter h q q . Remark 3.3.
In Proposition 3.1, ( v ), the asymptotic covariance between same-entrywavelet variance terms is always nontrivial, irrespective of the values of fractal connectivityparameters. For example, when (1 ,
2) =: ( q , q ) = ( q , q ), it is given by G jj (cid:48) (1 , , ,
2) = 2 − ( j + j (cid:48) )2 Φ jj (0)Φ j (cid:48) j (cid:48) (0) · (cid:26) φ jj (cid:48) (1 , , , , δ > φ jj (cid:48) (1 , , ,
2) + φ jj (cid:48) (1 , , , , δ = 0 . This is not true for terms associated with different pairs of indices. For example, if ( q , q ) =(1 , (cid:54) = ( q , q ) = (2 , G jj (cid:48) (1 , , ,
2) = 2 − ( j + j (cid:48) )2 Φ jj (0)Φ j (cid:48) j (cid:48) (0) · (cid:26) φ jj (cid:48) (1 , , , , δ = 0;0 , δ > . In other words, the phenomenon of the asymptotic decorrelation of wavelet variance termsis only observed for instances involving departures from fractal connectivity.3.3.
Estimation of the scaling exponents.
Definition of the estimators.
As in the univariate case, the fact that the waveletvariance (3.7) satisfies the Hadamard scaling relation in the coarse scale limit (see (3.15))points to the development of a log-regression regression method based on the sample waveletvariance across scales.Estimators can be defined in a standard way by means of the log-regression relationsˆ α q = j (cid:88) j = j w j log S ( qq ) n ( a ( n )2 j ) , ≤ q ≤ m, (3.21) ˆ α q q = j (cid:88) j = j w j log | S ( q q ) n ( a ( n )2 j ) | , ≤ q < q ≤ m, which are well-defined with probability 1, where w j , j = j , . . . , j , are linear regressionweights satisfying(3.22) j (cid:88) j = j w j = 0 , j (cid:88) j = j jw j = 1(e.g., [5, 45]). The derived estimator(3.23) ˆ δ q q = ˆ α q q + ˆ α q q − ˆ α q q , q (cid:54) = q , for the parameter (1.7) will be used in Section 5 in the construction of a hypothesis test forfractal connectivity (see also [11, 22, 47]). With the standard Gaussian asymptotics for thesample wavelet variances { S n ( a ( n )2 j ) } j = j ,...,j established in Proposition 3.1, the system ofequations (3.21) and (3.23) is expected to yield efficient estimators. In fact, this is provedin the next section.3.3.2. Asymptotic distribution.
In the next result, Theorem 3.1, we build upon Proposition3.1 to show that the vector ( (cid:98) α q q ) q ,q =1 ,...,m is an asymptotically unbiased and Gaussianestimator of the scaling exponents ( α q q ) q ,q =1 ,...,m (whence the analogous statement holdsfor the estimator (3.23)). As in the aforementioned proposition, the decorrelation propertyof the wavelet transform contributes to the Gaussianity of the asymptotic distribution. Theorem 3.1.
Let w j , . . . , w j be the weight terms in (3.21) and let G be as in (3.19) .Suppose (3.24) ρ q q (cid:54) = 0 , q , q = 1 , . . . , m, q (cid:54) = q (see (2.3) ). Then, as n → ∞ , (3.25) vec S (cid:104)(cid:16) √ na ( n ) δ q q +1 / ( (cid:98) α q q − α q q ) (cid:17) q ,q =1 ,...,m (cid:105) d → N m ( m +1)2 (0 , M GM ∗ ) . Remark 3.4.
Condition (3.24) is needed to ensure the consistency of the estimator. Fur-thermore, it is clear that ρ qq > q = 1 , . . . , m , since otherwise the process is not properor its spectral density is not positive semidefinite a.e. Remark 3.5.
The convergence rate in (3.25) depends on the unknown fractal connectivityparameters δ ·· . In practice, the latter can be replaced by their estimates or, in some cases,ignored, since they exponentiate the slow growth term a ( n ). ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 11
Numerical simulation setting.
We conducted Monte Carlo experiments over 1000 in-dependent realizations of bivariate HfBm. Though several parameter settings were tested,results are reported only for two representative cases: fractally connected ideal-HfBm( α , α , ρ ) = (0 . , . , . α , α , δ , ρ ) = (0 . , . , . , . . For thewavelet analysis, least asymmetric Daubechies wavelet with N Ψ = 2 vanishing momentswere used [38]. Estimation by means of weighted linear regression was performed on theoctave range ( j , j ) = (3 , log n − N Ψ ). This choice of regression range (with fine scale j fixed and j ∼ log n ) is used here in order to, first, obtain results that are consistent andcomparable with those reported in the literature for the estimation of univariate scalingexponents, cf., e.g., [5] and references therein, and, second, avoid the issue of fixing h min , h max and (cid:36) in (3.6).Monte Carlo parameter estimates ( (cid:98) α , (cid:98) α , (cid:98) α , (cid:98) δ ) are reported in Figures 1 and 2 interms of bias, standard deviation, skewness and kurtosis, as functions of (the log of) samplesizes n = { , , , , , } .3.3.4. Finite sample performance.
Figures 1 and 2 show that for all four estimates( ˆ α , ˆ α , ˆ α , ˆ δ ), the bias becomes negligible as the sample size grows. It can also beseen that the bias is slightly larger for non-fractally connected data, notably for the param-eter δ .The figures further show that standard deviations for all estimates decrease as n − / (the latter trend being plotted as superimposed dashed lines), with no significant differencebetween fractally and non-fractally connected instances. For ( ˆ α , ˆ α , ˆ α ), there is nodetectable dependence of the outcomes on the actual values of the parameters themselves(in accordance with theoretical calculations reported in Table 3). Accordingly, one observesthat the variances for ˆ α and ˆ α are identical and depend on the sample size, but noton any parameter of the stochastic process. For ˆ α , the variance does not depend on α ,but results not shown demonstrate that it does vary with ρ , as confirmed by approximatecalculations (see (4.9)). For ˆ δ , further simulations not displayed indicate that Var ˆ δ dependson δ roughly proportionally to 1 + δ . Departures from fractal connectivity tend to imply adecrease in the variance of ˆ δ ; this is a counterintuitive phenomenon that has yet to be fullyexplained.In regard to convergence in distribution, Figures 1 and 2 show that the Monte Carloskewness and kurtosis estimates for ( α , α ) are very close to 0 even at very small samplesizes, both with and without fractal connectivity. This characterizes a fast convergenceto limiting normal distributions, and lies in agreement with the asymptotics presented inTheorem 3.1. However, the weak convergence is observed to be much slower in practicefor the cross-exponents ( α , δ ), yet with no noticeable difference between fractally andnon-fractally connected instances, i.e., different values of the parameter δ .4. Confidence intervals
To complement the asymptotic and finite sample results in Sections 3.3.2 and 3.3.4, wenow construct confidence intervals for the estimators (3.21) and (3.23). In particular, weinvestigate the effect of omitting the covariance among sample wavelet variance terms, aswell as the dependence of the asymptotic variance in (3.25) and confidence intervals on the
12 MULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY log Sample Size10 12 14 16 18 20 B i a s f o r α -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( α X ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 K u r t o s i s ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 B i a s f o r α -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( α Y ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 K u r t o s i s ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 B i a s f o r α -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( α XY ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 K u r t o s i s ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 B i a s f o r δ -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( δ , ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( δ ) -1-0.500.51 log Sample Size
10 12 14 16 18 20 K u r t o s i s ( δ ) -1-0.500.51 Figure 1.
Estimation performance and asymptotic normality.
Study based on 1000 Monte Carlo realizations of (fractally connected) ideal-HfBm with parameters ( ↵ , ↵ , , ⇢ ) = (0 . , . , , . Asymptotic covariance.
In Theorem 4.1 below, we provide a finite sample approxi-mation for the covariance between the logarithms of sample wavelet variances, coupled withthe convergence rate of the covariance approximation to its asymptotic expression in termsof ideal-HfBm covariances (3.14). As in Theorem 3.1, in part due to the decorrelation prop-erty of the wavelet transform, the estimators’ second order properties are comparable tothose observed in the i.i.d. case (e.g., [35], p. 430). In particular, the finite sums provided inTheorem 4.1 mathematically justify heuristic expressions obtained from Taylor expansions(c.f. (4.5) below), and the theorem also accurately quantifies the approximation error as afunction of the sample size and scale.In the theorem, we use indicator functions to regularize log | S ( q q ) n ( a ( n )2 j ) | around theorigin. The truncation converges fast to zero and is of no impact in practice. Indeed, Figure 1.
Estimation performance and asymptotic normality.
Study based on 1000 Monte Carlo realizations of (fractally connected) ideal-HfBm with parameters ( α , α , δ , ρ ) = (0 . , . , , . Asymptotic covariance.
In Theorem 4.1 below, we provide a finite sample approxi-mation for the covariance between the logarithms of sample wavelet variances, coupled withthe convergence rate of the covariance approximation to its asymptotic expression in termsof ideal-HfBm covariances (3.14). As in Theorem 3.1, in part due to the decorrelation prop-erty of the wavelet transform, the estimators’ second order properties are comparable tothose observed in the i.i.d. case (e.g., [35], p. 430). In particular, the finite sums provided inTheorem 4.1 mathematically justify heuristic expressions obtained from Taylor expansions(c.f. (4.5) below), and the theorem also accurately quantifies the approximation error as afunction of the sample size and scale.In the theorem, we use indicator functions to regularize log | S ( q q ) n ( a ( n )2 j ) | around theorigin. The truncation converges fast to zero and is of no impact in practice. Indeed,one can redefine the estimator (3.21) with the indicators and obtain the same asymptoticdistribution (3.25). ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 13
MULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 13 log Sample Size10 12 14 16 18 20 B i a s f o r α -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( α X ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 K u r t o s i s ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 B i a s f o r α -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( α Y ) -11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 K u r t o s i s ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 B i a s f o r α -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( α XY ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 K u r t o s i s ( α ) -1-0.500.51 log Sample Size10 12 14 16 18 20 B i a s f o r δ -0.1-0.08-0.06-0.04-0.0200.020.04 log Sample Size
10 12 14 16 18 20 l og s t d ( δ , ) -12-11-10-9-8-7-6-5-4-3-2-1 log Sample Size10 12 14 16 18 20 S k e w ne ss ( δ ) -1-0.500.51 log Sample Size
10 12 14 16 18 20 K u r t o s i s ( δ ) -1-0.500.51 Figure 2.
Estimation performance and asymptotic normality.
Study based on 1000 Monte Carlo realizations of (non-fractally connected)HfBm with parameters ( ↵ , ↵ , , ⇢ ) = (0 . , . , . , . Theorem 4.1.
Let q , q , q , q = 1 , . . . , m be generic indices, and suppose condition (3.24) holds. Consider S n ( a ( n )2 j ) = { S ( q q ) n ( a ( n )2 j ) } , W n ( a ( n )2 j ) = { W ( q q ) n ( a ( n )2 j ) } as in (3.10) and jj q q ( · ) as in (3.14) . For notational simplicity, denote (4.1) n ⇤ = na ( n )2 j + j and also recall the definition of $ in (2.7) . Consider a sequence { r n } n N ✓ (0 , / satisfying (4.2) 1 r n = O ⇣ exp n C ⇣ n ⇤ a ( n ) h max h min ⌘ ⇠ o⇣ a ( n ) h max h min n ⇤ ⌘ / ⌘ Figure 2.
Estimation performance and asymptotic normality.
Study based on 1000 Monte Carlo realizations of (non-fractally connected)HfBm with parameters ( α , α , δ , ρ ) = (0 . , . , . , . Theorem 4.1.
Let q , q , q , q = 1 , . . . , m be generic indices, and suppose condition (3.24) holds. Consider S n ( a ( n )2 j ) = { S ( q q ) n ( a ( n )2 j ) } , W n ( a ( n )2 j ) = { W ( q q ) n ( a ( n )2 j ) } as in (3.10) and Φ jj (cid:48) q q ( · ) as in (3.14) . For notational simplicity, denote (4.1) n ∗ = na ( n )2 j + j (cid:48) and also recall the definition of (cid:36) in (2.7) . Consider a sequence { r n } n ∈ N ⊆ (0 , / satisfying (4.2) 1 r n = O (cid:16) exp (cid:110) C (cid:16) n ∗ a ( n ) h max − h min (cid:17) − ξ (cid:111)(cid:16) a ( n ) h max − h min n ∗ (cid:17) / (cid:17) for any C > and any < ξ < . Then, Cov (cid:104) log | S ( q q ) n ( a ( n )2 j ) | {| W ( q q n ( a ( n )2 j ) | >r n } , log | S ( q q ) n ( a ( n )2 j (cid:48) ) | {| W ( q q n ( a ( n )2 j (cid:48) ) | >r n } (cid:105) = (log e ) − ( j + j (cid:48) ) Φ jjq q (0)Φ j (cid:48) j (cid:48) q q (0) { O ( a ( n ) − (cid:36) ) } (cid:110) a ( n ) h q q + h q q ) − h q q + h q q ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) ) (cid:105) + a ( n ) h q q + h q q ) − h q q + h q q ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) ) (cid:105) (4.3) + o (cid:16) a { h q q + h q q ,h q q + h q q }− h q q + h q q ) n ∗ (cid:17)(cid:111) + O (cid:16)(cid:16) a ( n ) h max − h min n ∗ (cid:17) (cid:17) . Remark 4.1.
Note that the finite summations on the right-hand side of (4.3) converge as n → ∞ . For example,1 n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) ) → φ jj (cid:48) ( q , q , q , q ) , where φ jj (cid:48) ( · ) is defined by (3.18) (this can be shown by the same argument for establishing(A.16) and (A.25) in the proof of Proposition 3.1, ( iv )). Remark 4.2.
In (4.3), the main residual term satisfies(4.4) O (cid:16)(cid:16) a ( n ) h max − h min n ∗ (cid:17) (cid:17) = o (1)under the condition (3.6). In practice, the modeling of instances involving extreme de-viations from fractal connectivity ( δ ·· close to 1) may require greater regularity from thefunctions g ·· (see (2.3)) or the wavelet basis. Otherwise, for conservative choices of theregularity parameters (e.g., (cid:36) = 2 and β = 1 + ε for small ε > Remark 4.3.
In (3.5), we assume that an HfBm continuous time sample path is available.However, it is well known that, under mild assumptions, the availability of discrete timeobservations does not generally alter the nature of the asymptotic distribution for waveletestimators (see, for instance, [13], Section III; [45], Section 3.2; [2], Section C). Establishingthis for the case of HfBm is a topic for future work.
Remark 4.4.
When q = q = q , for the univariate case the expansion (4.3) appearsimplicitly in [41], expression (86). Also, when computing the moments of W ( qq ) n ( a ( n )2 j ),the truncation is unnecessary (see Remark C.2 below).4.2. Variances and covariances of ˆ α q q and ˆ δ q q . Closed-form approximations.
Turning back to expression (4.3), by ignoring the scal-ing factor ( a ( n ) ≡ (cid:104) log | S ( q q ) n (2 j ) | , log | S ( q q ) n (2 j (cid:48) ) | (cid:105) ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 15 ≈ (log e ) − ( j + j (cid:48) ) Φ jjq q (0)Φ j (cid:48) j (cid:48) q q (0) 1 n j (cid:110) n j (cid:48) n j (cid:88) k =1 n j (cid:48) (cid:88) k (cid:48) =1 Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) )(4.5) + 1 n j (cid:48) n j (cid:88) k =1 n j (cid:48) (cid:88) k (cid:48) =1 Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) q q (2 j (cid:48) k (cid:48) − j k ) (cid:111) , where n j := n j . Note that the approximation (4.5) is a finite sample one, and ignores the potential asymp-totic decorrelation effect stemming from the shifting scaling factor a ( n ) (see Proposition3.1, ( v ), and (4.3)). Consider the normalization(4.6) r q q ( j, k ; j (cid:48) , k (cid:48) ) := Φ jj (cid:48) q q (2 j k − j (cid:48) k (cid:48) ) (cid:113) Φ jjq q (0)Φ j (cid:48) j (cid:48) q q (0) ∈ [ − , . Expressions (3.21), (4.5) and (4.6) yield the covariance approximation(4.7) Cov [ ˆ α q q , ˆ α q q ] ≈ (log e) j (cid:88) j,j (cid:48) = j w j w j (cid:48) n j n j (cid:48) n j − (cid:88) k =0 n j (cid:48) − (cid:88) k (cid:48) =0 r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j, k ) r q q ( j (cid:48) , k (cid:48) ; j (cid:48) , k (cid:48) )+ r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j, k ) r q q ( j (cid:48) , k (cid:48) ; j (cid:48) , k (cid:48) ) . In particular, (4.7) further allows us to compute(4.8) Var (cid:104) ˆ δ q q (cid:105) = 14 (Var [ ˆ α q q ] + Var [ ˆ α q q ]) + Var [ ˆ α q q ]+ 12 Cov [ ˆ α q q , ˆ α q q ] − Cov [ ˆ α q q , ˆ α q q ] − Cov [ ˆ α q q , ˆ α q q ] . The variances (4.8) will be used in Section 5 in the construction of a test for fractal con-nectivity, i.e., for the hypothesis H : δ q q ≡
0. Table 1 summarizes the closed-form ap-proximations for Var [ ˆ α q q ], Var [ ˆ α q q ], Cov [ ˆ α q q , ˆ α q q ] and Cov [ ˆ α q q , ˆ α q q ] establishedin (4.7). Var[ˆ α qq ](log e) (cid:80) j j,j (cid:48) = j w j w j (cid:48) n j n j (cid:48) (cid:80) n j ,n (cid:48) j k,k (cid:48) =1 r qq ( j, k ; j (cid:48) , k (cid:48) ) Var [ ˆ α q q ] (log e) (cid:80) j j,j (cid:48) = j w j w j (cid:48) n j n j (cid:48) (cid:80) n j ,n (cid:48) j k,k (cid:48) =1 r q q ( j,k ; j (cid:48) ,k (cid:48) ) r q q ( j (cid:48) ,k (cid:48) ; j,k )+ r q q ( j,k ; j (cid:48) ,k (cid:48) ) r q q ( j,k ; j (cid:48) ,k (cid:48) ) r q q ( j,k ; j,k ) r q q ( j (cid:48) ,k (cid:48) ; j (cid:48) ,k (cid:48) )Cov [ ˆ α q q , ˆ α q q ] (log e) (cid:80) j j,j (cid:48) = j w j w j (cid:48) n j n j (cid:48) (cid:80) n j ,n (cid:48) j k,k (cid:48) =1 r q q ( j,k ; j (cid:48) ,k (cid:48) ) r q q ( j,k ; j,k ) r q q ( j (cid:48) ,k (cid:48) ; j (cid:48) ,k (cid:48) )Cov [ ˆ α q q , ˆ α q q ] (log e) (cid:80) j j,j (cid:48) = j w j w j (cid:48) n j n j (cid:48) (cid:80) n j ,n (cid:48) j k,k (cid:48) =1 r q q ( j,k ; j (cid:48) ,k (cid:48) ) r q q ( j,k ; j (cid:48) ,k (cid:48) ) r q q ( j,k ; j,k ) r q q ( j (cid:48) ,k (cid:48) ; j (cid:48) ,k (cid:48) ) Table 1.
Closed-form expression for
Cov [ ˆ α q q , ˆ α q q ] . Impact of inter- and intra-scale correlations.
The expression of Cov [ ˆ α q q , ˆ α q q ] canfurther be split into three terms, namely,(4.9) Cov [ ˆ α q q , ˆ α q q ] ≈ (log e) j (cid:88) j = j (cid:34) w j n j (cid:18) r q q ( j, j, r q q ( j, j, (cid:19) + w j n j (cid:88) k (cid:88) k (cid:48) (cid:54) = k r q q ( j, k ; j, k (cid:48) ) r q q ( j, k ; j, k (cid:48) ) + r q q ( j, k ; j, k (cid:48) ) r q q ( j, k ; j, k (cid:48) ) r q q ( j, k ; j, k ) r q q ( j, k (cid:48) ; j, k (cid:48) )+ (cid:88) j (cid:48) (cid:54) = j w j w j (cid:48) n j n j (cid:48) (cid:88) k (cid:88) k (cid:48) r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j (cid:48) , k (cid:48) ) + r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j (cid:48) , k (cid:48) ) r q q ( j, k ; j, k ) r q q ( j (cid:48) , k (cid:48) ; j (cid:48) , k (cid:48) ) (cid:35) , where the first term in the sum over j reflects the variance only of wavelet coefficients, thesecond term the covariance of wavelet coefficients at a given scale, and the third term, thecovariance of wavelet coefficients at different scales. In other words, if wavelet coefficientswere independent, the second and third terms would equal zero.The relative contributions of the three terms to the final variances are quantified bymeans of Monte Carlo simulations conducted following the same protocol and settings asthose described in Section 3.3.3. Table 2, reporting the relative contributions of each of thethree terms for various sample sizes under fractal connectivity (i.e., δ q q ≡ δ q q > n Var [ (cid:98) α ] Var [ (cid:98) α ] Var [ (cid:98) α ] Var (cid:104)(cid:98) δ (cid:105) var × (cid:98) α , (cid:98) α ] Cov [ (cid:98) α , (cid:98) α ] Cov [ (cid:98) α , (cid:98) α ]co × Table 2.
Relative contributions of the three terms in (4.9) toVar (cid:2) (cid:98) α ( · ) (cid:3) , Cov [ ˆ α q q , ˆ α q q ] and Var (cid:104)(cid:98) δ ( · ) (cid:105) for various sample sizes.(( α , α , δ , ρ ) = (0 . , . , , . j = 2 and j = { , , , } , n = { , , , } ).4.2.3. First order approximations for the variances and covariances of ˆ α q q . It is of interestto further examine the leading order approximations for the variances and covariances ofˆ α q q and δ q q , corresponding to neglecting all intra- and inter-scale correlations amongst ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 17 wavelet coefficients. The first order approximations neglecting all inter- and intra- scale cor-relations amongst wavelet coefficients of Cov [ ˆ α q q , ˆ α q q ] and of Var (cid:104) ˆ δ q q (cid:105) are summarizedin Table 3. Var[ˆ α qq ](log e) (cid:80) j j = j w j n j Var [ ˆ α q q ] (log e) (cid:80) j j = j w j n j (1 + r q q ( j, j, ) Cov [ ˆ α q q , ˆ α q q ] (log e) (cid:80) j j = j w j n j r q q ( j, j, Cov [ ˆ α q q , ˆ α q q ] (log e) (cid:80) j j = j w j n j Var [ ˆ δ q q ] (log e) (cid:80) j j = j w j n j (cid:16) ( r q q ( j, j,
0) + r q q ( j, j, ) − (cid:17) Table 3.
First-order approximations in (4.7) for the variances andcovariances of ˆ α q q and δ q q neglecting all inter- and intra-scale correla-tions amongst wavelet coefficients.The results show that Var [ ˆ α q q ] and Var [ ˆ α q q ], q (cid:54) = q , do not depend on the ac-tual values of the scaling exponents ( α q q , α q q , α q q ), which corroborates the numericalperformance reported in Section 3.3.4.Note that for an ideal-HfBm with fractal connectivity, it is straightforward to show that r qq ( j, j, ≡ r q q ( j, j, ≡ ρ q q when q (cid:54) = q . While Var [ ˆ α q q ] does not dependon correlations ρ q q , as expected, Var [ ˆ α q q ] , q (cid:54) = q does vary with ρ q q according to1 /ρ q q , showing that Var [ ˆ α q q ] → + ∞ when ρ q q →
0. This can be interpreted as the factthat when ρ q q →
0, the scaling exponent α q q , q (cid:54) = q , loses its meaning. Furthermore,Cov [ ˆ α q q , ˆ α q q ] depends on ρ as ρ q q , not surprisingly indicating that when ρ q q → α q q , ˆ α q q ] → (cid:104) ˆ δ q q (cid:105) is observed not to depend on theactual value of δ q q , while Var (cid:104) ˆ δ q q (cid:105) clearly depends on δ q q based on the numericalsimulations reported in Section 3.3.4. This can be interpreted as the fact that, for ˆ δ q q , thefirst order approximation (neglecting all intra- and inter-scale correlations amongst waveletcoefficients) is not sufficient to approximate well Var (cid:104) ˆ δ q q (cid:105) , as opposed to what is observedfor Var [ ˆ α q q ] and Var [ ˆ α q q ] , q (cid:54) = q .To finish with, Var (cid:104) ˆ δ q q (cid:105) varies with ρ q q as ρ q q + 1 /ρ q q −
2. This shows again thatwhen ρ q q →
0, parameter δ q q becomes irrelevant. Moreover, it also shows that when ρ q q → ±
1, Var [ δ q q ] →
0. This can be understood as the fact that when ρ q q → ±
1, adeparture from fractal connectivity is no longer permitted, as indicated by (2.11). Thus, ρ q q → ± δ q q →
0, which is then no longer a random variable.4.3.
Practical computation of the variances and covariances of ˆ α q q and δ q q . Computation of E [ d q ( j, k ) d q ( j (cid:48) , k (cid:48) )] and r q q ( j, k ; j (cid:48) , k (cid:48) ) . Evaluation of (4.7) re-quires knowledge of the covariance between wavelet coefficients (4.6), which will be de-veloped here explicitly for the discrete and the dyadic wavelet transforms. Let h ( k ) and g ( k ), k = 1 , · · · , L , be the coefficients of the high pass and low pass filters of the dis-crete wavelet transform, respectively, and let ↑ [ · ] and ↓ [ · ] be the dyadic upsamplingand decimation operators. The wavelet transform of the discrete-time process compo-nent X q ( k ) yields, at each scale j = 1 , . . . , J , sequences of approximation coefficients a q ( j, k ) and detail coefficients d q ( j, k ), q = 1 , . . . , m . The corresponding dyadic coeffi-cients are given by ˜ a q ( j, k ) = a q ( j, j k ) and ˜ d q ( j, k ) = d q ( j, j k ), respectively. Pick theinitialization a q (0 , k ) = X q ( k ) (see [4] for a discussion of the initialization of the discretewavelet transform). At scale j = 1, a q (1 , · ) = h ∗ a q (0 , · ) and d q (1 , · ) = g ∗ a q (0 , · ), atscale j = 2, a q (2 , · ) = ↑ [ h ] ∗ a q (1 , · ) = ↑ [ h ] ∗ h ∗ a q (0 , · ) and d q (2 , · ) = ↑ [ g ] ∗ a q (1 , · ) = ↑ [ g ] ∗ h ∗ a q (0 , · ). By iteration, we obtain the sequences of detail coefficients at each scale j = j (cid:48) , i.e.,(4.10) d q ( j (cid:48) , k ) = (cid:0) g j (cid:48) ∗ a q (0 , · ) (cid:1) ( k ) , d q ( j (cid:48) , k ) = d q ( j (cid:48) , j (cid:48) k ) , where g j (cid:48) = ↑ j (cid:48)− [ g ] ∗ (cid:0) j (cid:48) − ∗ j =0 ↑ j (cid:48) [ h ] (cid:1) . Now let γ h q q ( s, t ), 1 ≤ q ≤ q ≤ m , be the covariances of univariate fBm of indices h q q = α q q / γ q q ( s, t ) = ρ q q σ q σ q {| s | α q q + | t | α q q − | s − t | α q q } . Then (cf. [22]; note that a change in time scale would only result in a multiplicative constant,which we assume to be absorbed in ρ q q and which cancels out in the final expression for r q q ) E [ d q ( j, k ) d q ( j (cid:48) , k + τ )] / ( σ q σ q ) == (cid:88) p (cid:88) q g j ( p ) g j (cid:48) ( q ) E [ a q (0 , k − p ) a q (0 , k + τ − q )]= − ρ q q (cid:88) p (cid:88) q g j ( p ) g j (cid:48) ( q ) | − τ + q − p | α q q + ρ q q (cid:88) p g j ( p ) (cid:88) q g j (cid:48) ( q ) | k − p | α q q + ρ q q (cid:88) q g j (cid:48) ( q ) (cid:88) p g j ( p ) | k + τ − q | α q q = − ρ q q (cid:88) p (cid:88) q g j ( p ) g j (cid:48) ( q ) | τ − q + p | α q q ( p (cid:48) = p − q )= − ρ q q (cid:88) p (cid:48) (cid:88) q g j ( p (cid:48) + q ) g j (cid:48) ( q ) | τ − p (cid:48) | α q q = − ρ q q (cid:88) p (cid:48) (cid:88) q g j ( p (cid:48) − q )ˇ g j (cid:48) ( q ) | τ − p (cid:48) | α q q = − ρ q q (cid:88) p ( g j ∗ ˇ g j (cid:48) )( p ) | τ − p | α q q = − ρ q q (( g j ∗ ˇ g j (cid:48) ) ∗ η q q )( τ ) , where ˇ g j ( k ) = g j ( L − k ) , k = 1 , . . . , L,η q q ( τ ) = | τ | α q q . Consequently,(4.12) r q q ( j, k ; j (cid:48) , k (cid:48) ) = ρ q q (( g j ∗ ˇ g j (cid:48) ) ∗ η q q )( k (cid:48) − k ) (cid:112) (( g j ∗ ˇ g j ) ∗ η ii )(0) (( g j (cid:48) ∗ ˇ g j (cid:48) ) ∗ η ll )(0) ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 19 and, for the dyadic wavelet transform,(4.13) ˜ r q q ( j, k ; j (cid:48) , k (cid:48) ) = r q q ( j, j k ; j (cid:48) , j (cid:48) k (cid:48) ) . For given values of α q q and ρ q q , these expressions can be easily evaluated numerically.4.3.2. Practical estimation of (co)variances of ˆ α q q and of ˆ δ q q . Evaluating (4.7) and (4.8)for HfBm in practice requires the unknown parameter values α q q and ρ q q in (4.11),and hence we replace them by their estimates (cid:98) α q q and (cid:98) ρ q q . The former are definedin (3.21), and estimates ˆ ρ q q for ρ q q for q (cid:54) = q can be readily obtained as the cross-correlation coefficients of the difference processes (fGn) of the process components X q and X q However, note that the expressions for the (co)variances of ˆ α q q in the previous sectionsare derived assuming knowledge of the true parameter values and can only be expected tobe approximations when these are replaced by estimates. This will be studied numericallyin the next section.4.3.3. Assessment of the estimated (co)variances of ˆ α q q and of ˆ δ q q by means of MonteCarlo experiments. Monte Carlo studies were conducted following the protocol and settingsdescribed in Section 3.3.3, aiming to evaluate the quality of the estimated approximations(4.7) and (4.8) for the (co)variances of ˆ α q q and of ˆ γ q q . The simulations involved 1000independent realizations of each of two general instances of HfBm with m = 2 components,one using the true values of the parameters α q q and ρ q q , and the other, their estimatesˆ α q q and ˆ ρ q q . Four different sample sizes, n = { , , , } and three different values ρ = { . , . , . } are investigated for the set of exponents [ α , α , α ] = [0 . , . , . n = 2 and weak correlation ρ = 0 .
3, the quality of the approximations (4.7) and (4.8) is very good for the variancesof exponents α q , q = 1 ,
2, and satisfactory for the “cross” exponent α , the covarianceparameters and the connectivity parameter δ . Second, when the sample size n and corre-lation level ρ increase, the approximation of variances and covariances becomes excellent,with maximum errors of the order of 5% for n = 2 and strong correlation ρ = 0 . α q q and ˆ ρ . They indicate that replacing the true parameter values α q q and ρ with estimates has very little impact on the quality of approximations (4.7) and(4.8). Indeed, the average values of the (co)variance estimates are essentially equal to thoseobtained when using true parameter values.5. Statistical test for fractal connectivity
Procedure.
The mathematical and computational results in Sections 3 and 4 enableus to construct component-wise fractal connectivity tests, i.e., for the hypotheses H : δ q q = α q q + α q q − α q q = 0 , q (cid:54) = q . We assume ρ q q (cid:54) = 0 for q (cid:54) = q . As a consequence of Theorem 3.1, the distribution of ˆ δ q q under H can be approximated over finite samples byˆ δ q q ∼ N (cid:16) , Var (cid:104) ˆ δ q q (cid:105)(cid:17) under H , n ρ Ratio of (cid:112) − / − theo/MC est/MC0 . α ] 0 .
95 0 .
97 0 .
95 0 .
99 0 .
94 0 .
97 0 .
94 0 . α ] 0 .
82 0 .
85 0 .
90 0 .
93 0 .
83 0 .
86 0 .
90 0 . α , (cid:98) α ] 1 .
12 1 .
07 0 .
97 1 .
32 1 .
09 1 .
06 0 .
97 1 . α ] 0 .
78 0 .
82 0 .
88 0 .
91 0 .
80 0 .
82 0 .
88 0 . . α ] 0 .
98 0 .
98 0 .
96 0 .
98 0 .
98 0 .
97 0 .
96 0 . α ] 0 .
94 0 .
92 1 .
00 1 .
00 0 .
93 0 .
92 1 .
00 1 . α , (cid:98) α ] 1 .
00 0 .
93 0 .
97 1 .
05 0 .
99 0 .
92 0 .
97 1 . δ ] 0 .
88 0 .
87 0 .
98 0 .
97 0 .
88 0 .
87 0 .
98 0 . . α ] 0 .
97 0 .
95 0 .
98 0 .
98 0 .
96 0 .
95 0 .
98 0 . α ] 1 .
00 0 .
98 1 .
00 1 .
01 0 .
99 0 .
98 1 .
00 1 . α , (cid:98) α ] 1 .
01 0 .
99 1 .
01 1 .
02 1 .
00 0 .
98 1 .
01 1 . δ ] 0 .
93 0 .
92 0 .
95 0 .
97 0 .
95 0 .
93 0 .
96 0 . Table 4.
Estimation of Var (cid:2) α ( · ) (cid:3) . Square roots of ratios of meanof (co)variances computed using (4.7) and of Monte Carlo (co)variances:(4.7) evaluated using theoretical values α , α , α , ρ (left columns,labeled “theo/MC”) and estimates ˆ α , ˆ α , ˆ α , ˆ ρ (right columns, la-beled “est/MC”). ( ( α , α , δ , ρ ) = (0 . , . , , ρ ), j = 2 and j = { , , , } , n = { , , , } ).where, in turn, Var (cid:104) ˆ δ q q (cid:105) can be approximated by (4.8). Therefore, a simple two-sidedtest with significance level s = P (reject H | H true) can be defined as(5.1) d s = , if | ˆ δ q q | > (cid:114) Var (cid:104) ˆ δ q q (cid:105) φ − (1 − s / , otherwise , where φ − ( · ) is the inverse cumulative distribution function of the standard Normal distri-bution. In addition, the p -value of the test statistic (i.e., the probability of observing anabsolute value at least as large as | ˆ δ q q | for the test statistic under H ) is given by(5.2) p ( | ˆ δ q q | ) := 2 φ (cid:16) − | ˆ δ q q | (cid:46)(cid:114) Var (cid:104) ˆ δ q q (cid:105)(cid:17) . This test can be performed by evaluating (5.1) with an estimate for Var (cid:104) ˆ δ q q (cid:105) obtainedfrom the procedure detailed in Section 4.1.5.2. Monte Carlo assessment of the test performance.
We assess the performanceof the test by applying it to 1000 independent realizations of HfBm with exponent val-ues [ α , α ] = [0 . , .
6] and exponent values α detailed below for sample sizes n = { , , , } and correlation levels ρ = { . , . , . } . For simplicity of illustrationand without loss of generality, we consider again HfBm with m = 2 components. For eachrealization, the test decision (5.1) and the p -value (5.2) are evaluated using (4.8) with ap-proximations (4.7) to obtain an estimate of the Var (cid:104) ˆ δ (cid:105) . Estimates of the expected valuesof the test decisions and p -values, denoted by ˆ d s and ˆ p , are then obtained as the averagesover realizations of test decisions and p -values (5.1) and (5.2). ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 21
We now compare the performance of the proposed test, denoted hereinafter HFBM (notto be confused with the stochastic process HfBm), to that of the test put forward in [49].The latter relies on the intuition that the wavelet coherence function of two components ofa multivariate Gaussian scale invariant random process approximately behaves asΓ q q ( j ) = S ( q q ) n j ( j ) / (cid:113) S ( q q ) n j ( j ) S ( q q ) n j ( j ) (cid:39) ρ q q j ( α q q − α q q − α q q ) . The test itself, denoted WCF (for wavelet coherence function), is formulated without therigorous statistical framework developed above. Rather, it is built on the observation thatΓ q q ( j ) is the Pearson product-moment correlation coefficient of the time series d q ( j, · ) and d q ( j, · ), and hence that the Fisher’s z statistics of Γ q q ( j ), j = j , . . . , j , are approximatelyGaussian, with known variances and, in the case of fractal connectivity, with equal meansacross scales. The test for fractal connectivity is then formulated as the UMPI test for theequality of means of Gaussian random variables, cf. [49] for details.5.2.1. Performance under H . We first consider the case that H is true, i.e., ( α + α ) / α = 0 .
4. The significance level is set to s = 0 .
1, results are reported in Table 5 for theproposed test (top) and for the test in [49] (bottom). Note that, under H , averages of testdecisions ˆ d s should equal the preset significance level s , and averages of p -values ˆ p shouldequal . HFBM rejects H with slightly larger probability than the prescribed value s = 0 . d s and s never exceed 5%; similarly,average p -values are slightly below 0 .
5. For large sample size and large ρ , average testdecisions and p -values are very close to the theoretical values s = 0 . p = . Theseremarks are consistent with the results reported in Table 4, where a small but systematicunderestimation of Var (cid:104) ˆ δ (cid:105) for small sample sizes and ρ is observed.In contrast, the empirical significances ˆ d s of WCF strongly differ from the preset value s by values of up to 16%, and this difference is especially pronounced for large samplesizes for which one would expect the test to perform better. One reason for this poorperformance may lie in the fact that the test in [49] was designed for fGn, rather than fBm.Note that the asymptotic calculations developed above can be adapted to the easier caseof fGn (and, in principle, any other Gaussian process with stationary increments) withoutdifficulty by simply changing the covariance function γ q q ( s, t ) in the calculations leadingto the expressions (4.12) and (4.13).5.2.2. Test power.
We assess the power of the test under the alternative hypotheses H : δ = ( α + α ) / − α (cid:54) = 0 with δ = { . , . , . , . } . Yet, a direct comparison ofthe power of HFBM and WCF is only meaningful for identical rejection probabilities under H . Indeed, a test for which ˆ d s > s under H is expected to display an artificially largepower. Since the performance of HFBM and WCF differ, as discussed above (cf. Table 5),we therefore adjust, for each sample size and correlation level, the prescribed significanceto the value ˜ s for which the average rejection rate under H equals ˆ d ˜ s = s = 0 .
1. Using thisadjusted level of significance ˜ s , the power of the test is then estimated as the average of thetest decisions d ˜ s when H is true. Results are reported in Table 6 and yield the followingconclusions. First, the power of each test systematically increases with the magnitudes ofthe deviation from δ = 0, of the correlation level ρ and of the sample size n , as expected.Second, HFBM is systematically and significantly more powerful. Indeed, it enables us todetect a non-zero value for δ up to two times as often as WCF. For instance, for the small sample size of n = 2 and the low correlation level of ρ = 0 .
5, it permits the detectionof a deviation of 0 . δ = 0 with probability 0 .
5, as compared to aprobability of 0 .
22 for the test in [49].Overall, these results confirm that the proposed developments can be relevantly appliedin the assessment of scaling and fractal connectivity in multivariate time series.
HFBM – H : δ ≡ s = 0 . ρ = 0 . ρ = 0 . ρ = 0 . d s ˆ p
100 ˆ d s ˆ p
100 ˆ d s ˆ pn = 2 . .
45 14 . .
48 14 . . n = 2 . .
45 10 . .
47 11 . . n = 2 . .
45 11 . .
46 11 . . n = 2 . .
46 10 . .
48 11 . . H : δ ≡ s = 0 . ρ = 0 . ρ = 0 . ρ = 0 . d s ˆ p
100 ˆ d s ˆ p
100 ˆ d s ˆ pn = 2 . .
44 15 . .
46 16 . . n = 2 . .
42 16 . .
43 17 . . n = 2 . .
40 18 . .
41 22 . . n = 2 . .
36 22 . .
36 22 . . Table 5.
Test significance.
Mean test decisions and p -values for dif-ferent values of n and ρ under H : δ = ( α + α ) / − α = 0([ α , α ] = [0 . , . j = 2 and j = { , , , } , n = { , , , } )for the proposed test (top) and for the test in [49] (bottom). HFBM – H : δ (cid:54) = 0, ˆ d s = 0 . ρ = 0 . ρ = 0 . ρ = 0 . δ .
05 0 . .
15 0 . .
05 0 . .
15 0 . .
05 0 . .
15 0 . n = 2 . . . . . . . . . . . . n = 2 . . . . . . . . . . . . n = 2 . . . . . . . . . . . . n = 2 . . . . . . . . . . . . H : δ (cid:54) = 0, ˆ d s = 0 . ρ = 0 . ρ = 0 . ρ = 0 . δ .
05 0 . .
15 0 . .
05 0 . .
15 0 . .
05 0 . .
15 0 . n = 2 . . . . . . . . . . . . n = 2 . . . . . . . . . . . . n = 2 . . . . . . . . . . . . n = 2 . . . . . . . . . . . . Table 6.
Test power for adjusted significance ˆ d s = 0 . . Mean testdecisions and p -values for different values of n , ρ and α / alternativehypotheses H : δ = ( α + α ) / − α (cid:54) = 0 ([ α , α ] = [0 . , . j = 2and j = { , , , } , n = { , , , } ) for the proposed test (top) andfor the test in [49] (bottom). ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 23 Conclusion
The present contribution introduces a versatile class of multivariate stochastic processes,named Hadamard fractional Brownian motion (HfBm), for scale invariance modeling. HfBmprovides a stochastic framework within which cross-component scaling laws are not directlycontrolled by the scaling along the main diagonal, namely, HfBm is not necessarily fractallyconnected.Interestingly, the theoretical study of HfBm reveals that exact entry-wise scaling on bothauto- and cross-components and departures from fractal connectivity are mathematicallyincompatible. In other words, there is a dichotomy in multivariate scaling modeling: eitherthere is exact entry-wise scaling in every component combined with fractal connectivity, ordepartures from fractal connectivity are allowed at the price of approximate (i.e., asymp-totic) scaling on the cross-components.Our main mathematical results consist of an asymptotically normal, wavelet-based linearregression estimator for the scaling exponents, as well as asymptotically valid confidenceintervals with convenient mathematical expressions. Furthermore, the Taylor expansionsused in the development of the asymptotic confidence intervals lead to the construction ofpractical procedures for the numerical calculation of the variance of the estimates. Theseapproximate calculations enable the study of the ubiquitous issue of the impact of neglect-ing intra- or inter-scale correlations amongst wavelet coefficients in the computations ofvariances and covariances for the estimates. We also devised an asymptotically normal hy-pothesis test for fractal connectivity. Again, a major feature of the designed test procedureis the fact that it can be applied to a single observed HfBm data path.For both fractally and non-fractally connected instances, simulations demonstrate thesatisfactory performance of the estimators of the scaling and fractal connectivity parameters,even for small sample size data. The estimation bias is shown to be negligible, and thevariance decreases according to the inverse of the sample size. In addition, the practicalcomputations of approximated variances and covariances of the estimates are shown to beof excellent quality, irrespective of sample size, and the Monte Carlo significance levels andpowers are very close to their theoretical counterparts.The tools developed in the present contribution pave the way for novel analysis and mod-eling perspectives on multivariate scaling in real-world data, in the spirit of the accomplish-ments in [16]. Routines for the synthesis of HfBm, as well as for estimation, computation ofconfidence intervals and fractal connectivity testing will be made publicly available at timeof publication.
References [1] P. Abry, R. Baraniuk, P. Flandrin, R. Riedi, and D. Veitch. Multiscale nature of network traffic.
IEEESignal Processing Magazine , 19(3):28–46, 2002.[2] P. Abry and G. Didier. Wavelet estimation for operator fractional Brownian motion: supplementarymaterial. To appear in
Bernoulli , pages 1–15, 2015.[3] P. Abry and G. Didier. Wavelet estimation for operator fractional Brownian motion. To appear in
Bernoulli , pages 1–30, 2016.[4] P. Abry and P. Flandrin. On the initialization of the discrete wavelet transform.
IEEE Signal ProcessingLetters , 1(2):32–34, 1994.[5] P. Abry, P. Flandrin, M. Taqqu, and D. Veitch. Wavelets for the analysis, estimation and synthesis ofscaling data. In
Self-similar Network Traffic and Performance Evaluation , pages 39–88. Wiley, spring2000. [6] P. Abry, S. G. Roux, H. Wendt, P. Messier, A. G. Klein, N. Tremblay, P. Borgnat, S. Jaffard, B. Vedel,J. Coddington, and L. Daffner. Multiscale anisotropic texture analysis and classification of photographicprints: art scholarship meets image processing algorithms.
IEEE Signal Processing Magazine , 32(4):18–27, July 2015.[7] P. Abry and D. Veitch. Wavelet analysis of long-range dependent traffic.
IEEE Transactions on Infor-mation Theory , 44(1):2–15, 1998.[8] P. Abry, H. Wendt, and S. Jaffard. When Van Gogh meets Mandelbrot: multifractal classification ofpainting’s texture.
Signal Processing , 93(3):554–572, 2013.[9] S. Achard, D. S. Bassett, A. Meyer-Lindenberg, and E. Bullmore. Fractal connectivity of long-memorynetworks.
Physical Reviews E , 77(3):036104, 2008.[10] S. Akselrod, D. Gordon, F. A. Ubel, D. C. Shannon, A. C. Berger, and R. J. Cohen. Power spectrumanalysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control.
Science ,213(4504):220–222, 1981.[11] P.-O. Amblard and J.-F. Coeurjolly. Identification of the multivariate fractional Brownian motion.
IEEETransactions on Signal Processing , 59(11):5152–5168, 2011.[12] J.-M. Bardet. Testing for the presence of self-similarity of Gaussian time series having stationary incre-ments.
Journal of Time Series Analysis , 21(5):497–515, 2000.[13] J.-M. Bardet. Statistical study of the wavelet analysis of fractional Brownian motion.
IEEE Transactionson Information Theory , 48(4):991–999, 2002.[14] J.-M. Bardet, G. Lang, G. Oppenheim, A. Philippe, S. Stoev, and M.S. Taqqu. Semi-parametric es-timation of the long-range dependence parameter: A survey. In P. Doukhan, G. Oppenheim, andM. S. Taqqu, editors,
Theory and applications of Long-range dependence , pages 557–577, Boston, 2003.Birkh¨auser.[15] P.J. Brockwell and R.A. Davis.
Time Series: Theory and Methods . Springer, 2 edition, 1991.[16] P. Ciuciu, P. Abry, and B. J. He. Interplay between functional connectivity and scale-free dynamics inintrinsic fMRI networks.
NeuroImage , 95(186):248–263, 2014.[17] P. Ciuciu, G. Varoquaux, P. Abry, S. Sadaghiani, and A. Kleinschmidt. Scale-free and multifractaldynamic properties of fmri signals during rest and task.
Frontiers in Physiology , 3(186):1–18, 2012.[18] J.-F. Coeurjolly, P.-O. Amblard, and S. Achard. Wavelet analysis of the multivariate fractional Brownianmotion.
ESAIM: Probability and Statistics , 17:592–604, 2013.[19] H. Cram´er and M. R. Leadbetter.
Stationary and Related Stochastic Processes: Sample Function Prop-erties and their Applications . Courier Dover Publications, 1967.[20] G. Didier, H. Helgason, and P. Abry. Demixing multivariate-operator self-similar processes. In
Proc.IEEE Int. Conf. Acoust., Speech, and Signal Proc. (ICASSP) , pages 3671–3675, Brisbane, Australia,2015.[21] G. Didier and V. Pipiras. Integral representations and properties of operator fractional Brownian mo-tions.
Bernoulli , 17(1):1–33, 2011.[22] P. Flandrin. Wavelet analysis and synthesis of fractional Brownian motion.
IEEE Transactions onInformation Theory , 38:910 – 917, March 1992.[23] B. J. He. Scale-free brain activity: past, present, and future.
Trends in Cognitive Science , to appear,2014.[24] B. J. He, J. Zempel, A. Snyder, and M. Raichle. The temporal structures and functional significance ofscale-free brain activity.
Neuron , 66:353–69, 2010.[25] H. Helgason, V. Pipiras, and P. Abry. Fast and exact synthesis of stationary multivariate Gaussian timeseries using circulant embedding.
Signal Processing , 91(5):1123 – 1133, 2011.[26] H. Helgason, V. Pipiras, and P. Abry. Synthesis of multivariate stationary series with prescribed mar-ginal distributions and covariance using circulant matrix embedding.
Signal Processing , 91(8):1741 –1758, 2011.[27] J. Istas and G. Lang. Quadratic variations and estimation of the local H¨older index of a Gaussianprocess.
Annales de l’Institut Henri Poincar´e , 33(4):407–436, 1997.[28] C. R. Johnson Jr., E. Hendriks, I. J. Berezhnoy, E. Brevdo, S. M. Hughes, I. Daubechies, J. Li,E. Postma, and J. Z. Wang. Processing for artist identification: Computerized analysis of Vincent vanGogh’s painting brushstrokes.
IEEE Signal Processing Magazine (Special Section - Signal Processing inVisual Cultural Heritage) , 25:37–48, 2008.[29] G. A. Jones and J. M. Jones.
Elementary Number Theory . Berlin: Springer-Verlag, 1998.
ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 25 [30] S. Kechagias and V. Pipiras. Definitions and representations of multivariate long-range dependent timeseries.
Journal of Time Series Analysis , 36(1):1–25, 2015.[31] K. Kiyono, Z. R. Struzik, N. Aoyagi, and Y. Yamamoto. Multiscale probability density function analy-sis: non-Gaussian and scale-invariant fluctuations of healthy human heart rate.
IEEE Transactions onBiomedical Engineering , 53(1):95–102, Jan. 2006.[32] L. Kristoufek. Mixed-correlated ARFIMA processes for power-law cross-correlations.
Physica A- ,392(24):6484–6493, 2013.[33] L. Kristoufek. Can the bivariate Hurst exponent be higher than an average of the separate Hurstexponents?
Physica A , 431:124–127, 2015.[34] B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model selection.
Annalsof Statistics , pages 1302–1338, 2000.[35] E. Lehmann and G. Casella.
Theory of Point Estimation . Springer, 1998.[36] I. N. Lobato. A semiparametric two-step estimator in a multivariate long memory model.
Journal ofEconometrics , 90(1):129–153, 1999.[37] R. Lopes and N. Betrouni. Fractal and multifractal analysis: a review.
Medical Image Analyis , 13:634–649, 2009.[38] S. Mallat.
A Wavelet Tour of Signal Processing . Academic Press, San Diego, CA, 1998.[39] B. Mandelbrot. A multifractal walk down wall street.
Scientific American , 280(2):70–73, Feb. 1999.[40] B.B. Mandelbrot. Intermittent turbulence in self-similar cascades: divergence of high moments anddimension of the carrier.
Journal of Fluid Mechanics , 62:331–358, 1974.[41] E. Moulines, F. Roueff, and M. Taqqu. On the spectral density of the wavelet coefficients of long-memorytime series with application to the log-regression estimation of the memory parameter.
Journal of TimeSeries Analysis , 28(2):155–187, 2007.[42] F. S. Nielsen. Local Whittle estimation of multi-variate fractionally integrated processes.
Journal ofTime Series Analysis , 32(3):317–335, 2011.[43] P. Robinson. Multivariate Local Whittle estimation in stationary systems.
Annals of Statistics ,36(5):2508–2530, 2008.[44] G. Samorodnitsky and M. S. Taqqu.
Stable non-Gaussian random processes . Chapman and Hall, NewYork, 1994.[45] S. Stoev, V. Pipiras, and M. Taqqu. Estimation of the self-similarity parameter in linear fractionalstable motion.
Signal Processing , 82(12):1873–1901, 2002.[46] M. S. Taqqu. Fractional Brownian motion and long range dependence. In
Theory and Applications ofLong-Range Dependence (P. Doukhan, G. Oppenheim and M. S. Taqqu, eds.) , pages 5–38. Birkh¨auser,Boston, 2003.[47] D. Veitch and P. Abry. A wavelet-based joint estimator of the parameters of long–range dependence.
IEEE Transactions on Information Theory , 45(3):878–897, April 1999.[48] C. Vignat. A generalized Isserlis theorem for location mixtures of Gaussian random vectors.
Statisticsand Probability Letters , 82(1):67–71, 2012.[49] H. Wendt, A. Scherrer, P. Abry, and S. Achard. Testing fractal connectivity in multivariate long memoryprocesses. In
Proc. IEEE Int. Conf. Acoust., Speech, and Signal Proc. (ICASSP) , pages 2913–2916,Taipei, Taiwan, 2009.
Appendix: proofs
This appendix comprises three parts, i.e., A, B and C, which contain the proofs for Sec-tions 3.2, 3.3 and 4.1, respectively. In B and C, we assume throughout that the assumptionsof Theorems 3.1 and 4.1, respectively, hold.In the proofs, whenever convenient we will use the shorthand(.1) a = a ( n ) . For notational simplicity, we will assume throughout that σ q = 1, q = 1 , . . . , m . Since themain diagonal entries of an HfBm behave like a perturbed (univariate) fBm, throughoutthe appendix we only provide proofs for cross-entry components, i.e., when the indices q l are pairwise distinct, l = 1 , , ,
4. Whenever convenient we will use l = 1 , , , q l , respectively, and also write h q l q l = h l , h q l q p = h lp , In addition, without loss of generality it will be assumed that(.2) (cid:37) (12) ( a j ) > , n ∈ N , j ∈ N (see 3.8), since otherwise we can consider − (cid:37) (12) ( a j ) instead. We will also write log insteadof log , for visual clarity. In the proofs, C represents a generic constant whose value maychange from one line to the next. Appendix A. Section 3.2Proof of Proposition 3.1 : For simplicity, we will write Ξ jj (cid:48) ( · ) instead of Ξ jj (cid:48) ,a ( · )throughout the proof.To show ( i ), the change of variable s = 2 − j t − k in (3.5) and the harmonizable represen-tation of HfBm yield(A.1) Ξ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) = (cid:90) R e ia (2 j k − j (cid:48) k (cid:48) ) x f ( x ) (cid:98) ψ ( a j x ) (cid:98) ψ ( a j (cid:48) x ) dx, where f ( x ) := (cid:16) | x | − h q q +1 / g q q ( x ) (cid:17) q ,q =1 ,...,m . Statement ( ii ) is a consequence of theformula (A.1) with(A.2) k = k (cid:48) = 0 , j = j (cid:48) . Statement ( iii ) follows from Lemma C.1, ( ii ), below under condition (3.6).Turning to ( iv ), establishing (3.16) is equivalent to showing that2 − ( j + j (cid:48) ) / n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 (cid:16) Ξ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) (A.3)+ 2 − ( j + j (cid:48) ) / n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) → − ( j + j (cid:48) ) / gcd(2 j , j (cid:48) ) ∞ (cid:88) z = −∞ Φ jj (cid:48) ( z gcd(2 j , j (cid:48) )) ,n → ∞ .Consider the entry q = 1, q = 2 of the matrix-valued function Ξ jj (cid:48) as in (A.1). Basedon a change of variable y = ax , recast the first term of the sum on the left-hand side of(A.3) as(A.4) (cid:90) R e iyz | y | − (2 h +1) ρ (cid:110) g (cid:16) ya (cid:17) − (cid:111) (cid:98) ψ (2 j y ) (cid:98) ψ (2 j (cid:48) y ) dy =: (cid:90) R e iyz h a ( y ) dy. Fix
A > (cid:98) ψ ( l ) ( x ) = d l dx l ( (cid:60) (cid:98) ψ ( x ) + i (cid:61) (cid:98) ψ ( x )), l ≥
0. By (3.1), (3.4) and a Taylorexpansion with Lagrange residual of the real and imaginary parts of (cid:98) ψ , there exist functions λ , λ on [ − A, A ] such that (cid:98) ψ ( x ) = (cid:16) d N ψ dx N ψ (cid:60) (cid:98) ψ ( x ) (cid:12)(cid:12)(cid:12) λ ( x ) + i d N ψ dx N ψ (cid:61) (cid:98) ψ ( x ) (cid:12)(cid:12)(cid:12) λ ( x ) (cid:17) x N ψ N ψ ! . ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 27
Therefore, and extending this reasoning to (cid:98) ψ (cid:48) ( x ), (cid:98) ψ (cid:48)(cid:48) ( x ),(A.5) | (cid:98) ψ ( l ) ( x ) | = O ( | x | N ψ − l ) , x → , l = 0 , , . For h a as in (A.4), we now show that(A.6) h a , h (cid:48) a , h (cid:48)(cid:48) a are differentiable and h a (0) = h (cid:48) a (0) = h (cid:48)(cid:48) a (0) = 0 . We will only develop expressions for y >
0, since analogous developments hold for y < h a as in (A.4) as(A.7) h a ( y ) = y − (2 h +1) ρ ϑ a ( y ) . Hence,(A.8) h (cid:48) a ( y ) = ρ {− (2 h + 1) y − (2 h +2) ϑ a ( y ) + y − (2 h +1) ϑ (cid:48) a ( y ) } ,h (cid:48)(cid:48) a ( y ) = ρ { (2 h + 1)(2 h + 2) y − (2 h +3) ϑ a ( y ) − h + 1) y − (2 h +2) ϑ (cid:48) a ( y )(A.9) + y − (2 h +1) ϑ (cid:48)(cid:48) a ( y ) } . Note that, by conditions (3.1) and (2.7),(A.10) 2 N ψ + (cid:36) − (2 h + 1 + l ) − ≥ (cid:36) − h max − l > , l = 0 , , . Then, h a is also smooth around zero and h a (0) = 0. Moreover, by (A.5) and (A.10) with l = 0, for fixed n and | y | ≤ aε ,(A.11) | h a ( y ) | y ≤ Cy y − (2 h +1) (cid:16) ya (cid:17) (cid:36) y N ψ = Ca (cid:36) y N ψ + (cid:36) − (2 h +2) → , as y → + , where the constant C > n . Similarly, by (2.6), (A.5) and(A.10) with l = 1, | h (cid:48) a ( y ) | y ≤ C (cid:48) y (cid:110) y − (2 h +2) (cid:16) ya (cid:17) (cid:36) y N ψ + y − (2 h +1) (cid:104)(cid:16) ya (cid:17) (cid:36) − a y N ψ + (cid:16) ya (cid:17) (cid:36) y N ψ − (cid:105)(cid:111) (A.12) ≤ C (cid:48) y N ψ − (2 h +2)+ (cid:36) − a (cid:36) → . This proves (A.6), as desired. Next, note that ∂∂x e itx ψ ( t ) = ite itx ψ ( t ) , ∂ ∂x e itx ψ ( t ) = − t e itx ψ ( t )and itψ ( t ), t ψ ( t ) ∈ L ( R ) by the continuity of ψ and condition (3.2). Therefore, by thedominated convergence theorem, (cid:98) ψ (cid:48) ( x ) = C (cid:82) R e itx itψ ( t ) dt , (cid:98) ψ (cid:48)(cid:48) ( x ) = C (cid:48) (cid:82) R e itx ( − t ) ψ ( t ) dt for appropriate constants C, C (cid:48) ∈ R . Consequently, by (3.2),(A.13) max l =0 , , sup x ∈ R | (cid:98) ψ ( l ) ( x ) | ≤ C (cid:90) R | t | l | ψ ( t ) | dt < ∞ . So, fix z (cid:54) = 0. By (2.5) and (A.13),(A.14) lim | y |→∞ (cid:12)(cid:12)(cid:12) h a ( y ) (cid:12)(cid:12)(cid:12) = 0 = lim | y |→∞ (cid:12)(cid:12)(cid:12) h (cid:48) a ( y ) (cid:12)(cid:12)(cid:12) . Thus, in view of (A.4), (A.6), (A.8), (A.9), (A.10) (with l = 2) and (A.14), by integratingby parts twice we obtain (cid:12)(cid:12)(cid:12) Ξ jj (cid:48) ( az ) a h − Φ jj (cid:48) ( z ) (cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12) z (cid:90) R e izy h (cid:48)(cid:48) a ( y ) dy (cid:12)(cid:12)(cid:12) (A.15) ≤ Cz (cid:90) R {| y | − (2 h +3) | ϑ a ( y ) | + | y | − (2 h +2) | ϑ (cid:48) a ( y ) | + | y | − (2 h +1) | ϑ (cid:48)(cid:48) a ( y ) |} dy ≤ C (cid:48) z , where the last inequality is a consequence of (2.5), (A.5) and (A.13).Now consider the first summation term in (A.3). We proceed as in the proof of Proposition3.3, ( i ), in [2] to establish that(A.16) 1 n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 (cid:16) Ξ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) → , n → ∞ . We outline the main steps for the reader’s convenience. By Theorem 1.8 in Jones andJones [29], p. 10, the set of values r ∈ Z to the equation a j k − a j (cid:48) k (cid:48) = r , k, k (cid:48) ∈ Z , isgiven by gcd( a j , a j (cid:48) ) Z =: R . Therefore, a pair ( k, k (cid:48) ) ∈ Z is a solution to(A.17) a j k − a j (cid:48) k (cid:48) = gcd( a j , a j (cid:48) ) w for some w ∈ Z if and only if it is a solution to(A.18) 2 j k − j (cid:48) k (cid:48) = gcd(2 j , j (cid:48) ) w for the same w . Therefore, we can replace n with n ∗ in Lemmas B.2 and B.3, [2], andreexpress the first summation term on the left-hand side of (A.16) as(A.19) (cid:88) r ∈R∩ B jj (cid:48) ( n ∗ ) ξ r ( n ∗ ) n ∗ (cid:16) Ξ jj (cid:48) ( ar ) − Φ jj (cid:48) ( ar ) a h (cid:17) . In (A.19), B jj (cid:48) ( n ∗ ) is the range for r such that the pairs ( k, k (cid:48) ) satisfying (A.18) lie in theregion(A.20) 1 ≤ k ≤ j (cid:48) n ∗ , ≤ k (cid:48) ≤ j n ∗ , and ξ r ( n ∗ ) is the number of such solution pairs ( k, k (cid:48) ) given some r . Moreover, for anysufficiently large n , let k ∈ { , . . . , j (cid:48) n ∗ } be the smallest number such that ( k , k (cid:48) ( k )) ∈ N solves (A.17) (for some w ∈ Z ), where(A.21) k (cid:48) ( k ) := 2 j j (cid:48) k − gcd(2 j , j (cid:48) ) w j (cid:48) . From the proof of Lemma B.2, [2], the set A of such solutions to (A.18) has the form(A.22) A = (cid:110) ( k, k (cid:48) ) ∈ Z : k = k + 2 j (cid:48) gcd(2 j , j (cid:48) ) Z , k (cid:48) is given by (A.21) (cid:111) . In light of (A.22), define the function k ( z ) = k + j (cid:48) gcd(2 j , j (cid:48) ) z , z ∈ Z . In particular,( k (0) , k (cid:48) ( k (0))) is a solution pair for (A.18). Let R (cid:51) x = gcd(2 j , j (cid:48) )( n ∗ − k / j (cid:48) ). Then, by(A.22), ( k ( (cid:98) x (cid:99) ) , k (cid:48) ( k ( (cid:98) x (cid:99) ))) is the rightmost solution for (A.17) within the first-entry range ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 29 k = 1 , . . . , j (cid:48) n ∗ . Moreover, given r , the number of solution pairs in the region (A.20) is ξ r ( n ∗ ) = (cid:98) x (cid:99) + 1, where(A.23) (cid:98) x (cid:99) n − ∗ → gcd(2 j , j (cid:48) ) , n → ∞ . In addition, by (A.15), and (3.14),(A.24) (cid:12)(cid:12)(cid:12) Ξ jj (cid:48) ( ar ) − Φ jj (cid:48) ( ar ) a h (cid:12)(cid:12)(cid:12) ≤ Cr , r (cid:54) = 0 , lim n →∞ (cid:12)(cid:12)(cid:12) Ξ jj (cid:48) ( ar ) − Φ jj (cid:48) ( ar ) a h (cid:12)(cid:12)(cid:12) = 0 , r ∈ Z . By expression (A.19), the dominated convergence theorem and (A.23), the limit (A.16)holds.Next recall that, up to a change of sign, Φ jj (cid:48) ( z ) corresponds to the cross moment of thewavelet transform of a fBm. Thus, by an analogous procedure, we also obtain(A.25) 1 n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) → gcd(2 j , j (cid:48) ) ∞ (cid:88) z = −∞ Φ jj (cid:48) ( z gcd(2 j , j (cid:48) )) , n → ∞ . This establishes (A.3).To show statement ( v ), first note that o (cid:16) a { h + h ,h + h }− h + h ) n ∗ (cid:17)(cid:46) a ( h + h + h + h ) − h + h ) n ∗ = o (1) , which follows from (3.6) and the fact that(A.26) 2 max { h + h , h + h } ≤ h + h + h + h . Thus, by expression (C.42) for Cov[ W (12) n ( a j ) , W (34) n ( a j (cid:48) )] (established in the proof ofLemma C.5 below), √ n a,j a δ √ n a,j (cid:48) a δ Cov (cid:104) W (12) n ( a j ) , W (34) n ( a j (cid:48) ) (cid:105) = 2 j + j (cid:48) n ∗ a δ + δ Cov (cid:104) W (12) n ( a j ) , W (34) n ( a j (cid:48) ) (cid:105) = 2 − ( j + j (cid:48) )2 [Φ jj (0)Φ j (cid:48) j (cid:48) (0)(1 + O ( a − (cid:36) )) ] · (cid:110) a h + h ) a h + h + h + h n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )(A.27) + a h + h ) a h + h + h + h n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:111) + o (1) . By (A.26) and (A.27), statement ( v ) holds.The argument for showing ( vi ) is an adaptation of the proof of Theorem 3.1 in [2] (seealso [12], pp.510-513; [13], p.997; and [27], Lemma 2). For notational simplicity, we onlywrite the proof for m = 2, where the entries are indexed 1 and 2, and q , q = 1 , Y n := (cid:16) d ( a j , , d ( a j , , . . . , d ( a j , n a,j ) , d ( a j , n a,j ); . . . ; (A.28) d ( a j , , d ( a j , , . . . , d ( a j , n a,j ) , d ( a j , n a,j ) (cid:17) ∈ R Υ( n ) , where(A.29) Υ( n ) := 2 j (cid:88) j = j n a,j . Let(A.30) θ = ( θ j , . . . , θ j ) ∈ R J , where J = j − j + 1 and θ j = ( θ j, , θ j, , θ j, ) ∗ ∈ R , j = j , . . . , j . Now form theblock-diagonal matrix D n defined by(A.31)diag (cid:16) n a,j (cid:114) j Ω n,j , . . . , n a,j (cid:114) j Ω n,j (cid:124) (cid:123)(cid:122) (cid:125) n a,j ; . . . ; 1 n a,j (cid:114) j Ω n,j , . . . , n a,j (cid:114) j Ω n,j (cid:124) (cid:123)(cid:122) (cid:125) n a,j (cid:17) , where(A.32) Ω n,j = θ j, E [ d ( a j , ] a δ θ j, E [ d ( a j , d ( a j , a δ θ j, E [ d ( a j , d ( a j , θ j, E [ d ( a j , ] , j = j , . . . , j In (A.32), it can be understood that the slow growth factors for the main diagonal termsare a δ = a δ ≡
1. We would like to establish the limiting distribution of the statistic T n = j (cid:88) j = j θ ∗ j √ j vec S (cid:104)(cid:16) a δ q q (cid:17) q ,q =1 , ◦ W n ( a j ) (cid:105) = j (cid:88) j = j θ ∗ j √ j n a,j (cid:16) n a,j (cid:88) k =1 d ( a j , k ) E (cid:2) d ( a j , (cid:3) , a δ n a,j (cid:88) k =1 d ( a j , k ) d ( a j , k ) E [ d ( a j , d ( a j , , n a,j (cid:88) k =1 d ( a j , k ) E (cid:2) d ( a j , (cid:3) (cid:17) ∗ = Y ∗ n D n Y n , where it suffices to consider θ in (A.30) such that(A.33) θ ∗ Σ( H ) θ > H ) in (A.33) is obtained from (3.17) and canbe written in block form as Σ( H ) =: ( G jj (cid:48) ) j,j (cid:48) = j ,...,j , corresponding to block entries of thevector θ = ( θ j , . . . , θ j ) ∗ . Let Γ Y n = Cov [ Y n , Y n ], and consider the spectral decompositionΓ / Y n D n Γ / Y n = O Λ O ∗ , where Λ is diagonal with real, and not necessarily positive, eigenvalues(A.34) ξ i ( a j ) , i = 1 , . . . , Υ( n ) , and O is an orthogonal matrix. Now let Z ∼ N (0 , I Υ( n ) ). Then, T n d = Z ∗ Γ / Y n D n Γ / Y n Z = Z ∗ O Λ O ∗ Z d = Z ∗ Λ Z =: Υ( n ) (cid:88) i =1 ξ i ( a j ) Z i . Assume for the moment that(A.35) max i =1 ,..., Υ( n ) | ξ i ( a j ) | = o (cid:16)(cid:16) an (cid:17) / (cid:17) . ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 31
By (A.33) and (3.17), na Var [ T n ] = j (cid:88) j = j j (cid:88) j (cid:48) = j θ ∗ j (cid:110)(cid:114) na j (cid:114) na j (cid:48) Cov (cid:104) vec S (cid:104)(cid:16) a δ q q (cid:17) q ,q =1 , ◦ W n ( a j ) (cid:105) , vec S (cid:104)(cid:16) a δ q q (cid:17) q ,q =1 , ◦ W n ( a j (cid:48) ) (cid:105) (cid:105)(cid:111) θ j (cid:48) → j (cid:88) j = j j (cid:88) j (cid:48) = j θ ∗ j G jj (cid:48) θ j (cid:48) > . Therefore, there exists a constant
C > n , na Var [ T n ] ≥ C > i =1 ,..., Υ( n ) | ξ i ( a j ) | (cid:112) Var [ T n ] ≤ C (cid:48) (cid:16) na (cid:17) / max i =1 ,..., Υ( n ) | ξ i ( a j ) | → , n → ∞ . The claim (3.19) is now a consequence of Lemma B.4 in [2].So, we need to show (A.35). The first step is to establish the bound(A.36) sup u ∈ S Υ( n ) − | u ∗ Γ / Y n D n Γ / Y n u | ≤ C max j = j ,...,j n a,j (cid:107) Ω n,j (cid:107) sup u ∈ S Υ( n ) − u ∗ Γ Y n u . Let u ∈ S Υ( n ) − and let v = Γ / u . We can break up the vector v into two-dimensionalsubvectors v · , · to reexpress v = ( v j , , . . . , v j ,n a,j ; . . . ; v j , , . . . , v j ,n a,j ) ∗ . Based on theblock-diagonal structure of D n expressed in (A.31), | u ∗ Γ / Y n D n Γ / Y n u | = | v ∗ D n v | = (cid:12)(cid:12)(cid:12) j (cid:88) j = j n j (cid:88) l =1 v ∗ j,l Ω n,j n a,j √ j v j,l (cid:12)(cid:12)(cid:12) ≤ C j (cid:88) j = j n a,j (cid:88) l =1 n a,j √ j (cid:107) Ω n,j (cid:107) (cid:107) v j,l (cid:107) (A.37) ≤ C (cid:16) max j = j ,...,j n a,j √ j (cid:107) Ω n,j (cid:107) (cid:17) j (cid:88) j = j n a,j (cid:88) l =1 (cid:107) v j,l (cid:107) = C (cid:16) max j = j ,...,j n a,j √ j (cid:107) Ω n,j (cid:107) (cid:17) u ∗ Γ Y n u , where the constant C comes from a change of matrix norms and only depends on the fixeddimension m = 2. By taking sup u ∈ S Υ( n ) − on both sides of (A.37), we arrive at (A.36).The second step towards showing (A.35) consists of analyzing the asymptotic behaviorof the right-hand side of (A.36), as n → ∞ . For this, we will assume the result of LemmaC.1 below. So, note that(A.38) max j = j ,...,j n a,j (cid:107) Ω n,j (cid:107) ∼ C an a { h ,h ,h } a δ , n → ∞ . Moreover, by (C.14) below, the maximum eigenvalue of Γ Y n is bounded by (cid:107) Γ Y n (cid:107) ≤ Ca { h ,h } , where (cid:107) · (cid:107) is the matrix Euclidean norm. Therefore, in view of(A.38) and (C.14), the right-hand side of (A.36) is bounded by C a h max − h min)+1 n a δ . Inturn, the latter expression divided by (cid:112) an is bounded by C (cid:16) a h max − h min)+1 n (cid:17) / a δ . Bycondition (3.6), this implies (A.35), and as a result, also (3.19). (cid:3) Remark A.1.
As a consequence of expression (A.27) and the fact that a δ √ n a,j → W (12) n ( a j ) L ( P ) −→ , n → ∞ . Appendix B. Section 3.3Proof of Theorem 3.1 : Fix q , q = 1 , . . . , m . Based on (3.22), rewrite1 a δ q q (cid:114) na ( (cid:98) α q q − α q q )= j (cid:88) j = j w j a δ q q (cid:114) na log | W ( q q ) n ( a j ) | + j (cid:88) j = j w j a δ q q (cid:114) na (cid:110) log (cid:12)(cid:12)(cid:12) E (cid:104) S ( q q ) n ( a j ) (cid:105) (cid:12)(cid:12)(cid:12) a α q q − log | Φ jjq q (0) | (cid:111) (B.1) + 1 a δ q q (cid:114) na (cid:16) j (cid:88) j = j w j log | Φ jjq q (0) | − α q q (cid:17) . By Lemma C.1, ( ii ) (below), and an application of the mean value theorem, for some θ ( n ) > (cid:12)(cid:12)(cid:12) E (cid:104) S ( q q ) n ( a j ) (cid:105) (cid:12)(cid:12)(cid:12) /a α q q and | Φ jjq q (0) | the second term in the sum (B.1)can be bounded in absolute value by j (cid:88) j = j w j a δ q q (cid:114) na (cid:110) θ ( n ) (cid:12)(cid:12)(cid:12) (cid:12)(cid:12)(cid:12) E (cid:104) S ( q q ) n ( a j ) (cid:105) (cid:12)(cid:12)(cid:12) a α q q − | Φ jjq q (0) | (cid:12)(cid:12)(cid:12)(cid:111) = j (cid:88) j = j w j a δ q q (cid:114) na | Φ jjq q (0) | + o (1) C j a (cid:36) ≤ C (cid:48) (cid:114) na δ q q +1+2 (cid:36) → , as n → ∞ , where the limit is a consequence of condition (3.6). Also note that, after achange of variable 2 j x = y in the expression for Φ jjq q (0) (see (3.14)),Φ jjq q (0) = 2 jα q q ρ q q (cid:90) R | y | − ( α q q +1) | (cid:98) ψ ( y ) | dy =: 2 jα q q c q q ∈ R . Therefore, by (3.22), the third term in the sum (B.1) can be written as1 a δ q q (cid:114) na (cid:16) j (cid:88) j = j w j ( jα q q + log | c q q | ) − α q q (cid:17) = 0 . So, in regard to the first term in the sum (B.1), consider the weight matrix M ∈ M ( m ( m +1)2 , m ( m +1)2 J, R ) defined by (cid:16) j / w j I m ( m +1)2 ; 2 ( j +1) / w j +1 I m ( m +1)2 ; . . . ; 2 j / w j I m ( m +1)2 (cid:17) , where I m ( m +1)2 is a m ( m +1)2 × m ( m +1)2 identity matrix and J is given by (3.20). We wouldlike to show that(B.2) M (cid:16) vec S (cid:104)(cid:16) √ n a,j a δ q q (cid:17) q ,q =1 ,...,m ◦ log ◦| W n ( a j ) | (cid:105)(cid:17) j = j ,...,j d → N m ( m +1)2 (0 , M GM ∗ ) , ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 33 where log ◦| A | := (cid:16) log | A q q | (cid:17) q ,q =1 ,...,m for any m × m real matrix A , and the term post-multiplying the matrix M in (B.2) is a m ( m +1)2 J -dimensional random vector.For any 0 < r <
1, fix a pair q , q and an octave j , which specifies one of the en-tries of the random vector on the left-hand side of (B.2). Define the set A n = { ω :min q ,q W ( q q ) n ( a j ) > r } . Under (.2), Proposition 3.1, ( vi ), implies that P ( A n ) → n → ∞ . Thus, for large enough n , in the set A n the mean value theorem gives thealmost sure expression(B.3) R (cid:51) j/ w j √ n a,j a δ q q log | W ( q q ) n ( a j ) | = 2 j/ w j √ n a,j a δ q q ( W ( q q ) n ( a j ) − θ + ( W ( q q ) n ( a j ))for a random variable θ + ( W ( q q ) n ( a j )) between W ( q q ) n ( a j ) and 1. Since W ( q q ) n ( a j ) P → θ + ( W ( q q ) n ( a j )) P →
1. By considering (B.3) for all 1 ≤ q ≤ q ≤ m and j = j , . . . , j ,relation (B.2) is now a consequence of Proposition 3.1, ( vi ), and Slutsky’s theorem. (cid:3) Appendix C. Section 4.1
C.1.
General results.
For j = j , . . . , j and n ∈ N , consider the jointly Gaussian vector Y n = (cid:16) d ( a j , (cid:112) (cid:37) (12) ( a j ) , . . . , d ( a j , n a,j ) (cid:112) (cid:37) (12) ( a j ) , . . . , d ( a j , (cid:112) (cid:37) (12) ( a j ) , . . . , d ( a j , n a,j ) (cid:112) (cid:37) (12) ( a j ) , (C.1) d ( a j , (cid:112) (cid:37) (12) ( a j ) , . . . , d ( a j , n a,j ) (cid:112) (cid:37) (12) ( a j ) , . . . , d ( a j , (cid:112) (cid:37) (12) ( a j ) , . . . , d ( a j , n a,j ) (cid:112) (cid:37) (12) ( a j ) (cid:17) ∗ ∈ R Υ( n ) (see (A.29) for the definition of Υ( n )). Let(C.2) Γ Y n = E [ Y n Y ∗ n ] = O Λ O ∗ = O diag( λ , Y , . . . , λ Υ( n ) , Y ) O ∗ , O ∈ O (Υ( n )) , be the associated covariance matrix and its matrix spectral decomposition. The followinglemma provides the finite-sample and asymptotic properties of the covariance structure ofwavelet coefficients, both from Ξ jj (cid:48) and Y n . Note that such covariance structure does notin general correspond to a multivariate stationary stochastic process when multiple octaves j are considered. Lemma C.1.
For j = j , . . . , j , q , q = 1 , . . . , m , and n ∈ N , the following statementshold. ( i ) Consider Ξ jj (cid:48) ,a ( a ( n )(2 j k − j (cid:48) k (cid:48) )) = E (cid:104) D ( a ( n )2 j , k ) D ( a ( n )2 j (cid:48) , k (cid:48) ) ∗ (cid:105) (see (3.12) ).For j = j (cid:48) and n ∈ N , there is a continuous spectral density f ψ,n ( x ) q q such thatwe can write (C.3) Ξ jj,aq q ( a ( n )2 j ( k − k (cid:48) ))( a ( n )2 j ) h q q = (cid:90) π − π e i ( k − k (cid:48) ) x f ψ,n ( x ) q q dx, k, k (cid:48) ∈ Z . Moreover, (C.4) | f ψ,n ( x ) q q | ≤ C, x ∈ ( − π, π ] , for a constant C > that does not depend on n . ( ii ) Let Φ jj (cid:48) ·· ( z ) be as in (3.14) , z ∈ Z , and fix any < ξ < . Then, (C.5) (cid:12)(cid:12)(cid:12) Ξ jj (cid:48) ,aq q ( za ( n )) a ( n ) h q q − Φ jj (cid:48) q q ( z ) (cid:12)(cid:12)(cid:12) ≤ Ca ( n ) (cid:36) , n → ∞ , where (cid:36) and β are defined in expressions (2.6) and (3.3) , respectively. ( iii ) Let Y n be as in (C.1) , and let λ i, Y , i = 1 , . . . , Υ( n ) , be the eigenvalues of thecovariance matrix Γ Y n (see (A.29) and (C.2) ). Then, for some C > , (C.6) max i =1 ,..., Υ( n ) λ i, Y ≤ Ca ( n ) { h ,h }− h ) . Proof.
We will use the shorthand q = 1 and q = 2 throughout the proof.We first show ( i ). By making the change of variable y = a j x in (A.1), we obtain (C.3)with(C.7) f ψ,n ( x ) := ∞ (cid:88) l = −∞ | (cid:98) ψ ( x + 2 πl ) | | x + 2 πl | h +1 ρ g (cid:16) x + 2 πla j (cid:17) . Moreover, by (3.1)–(3.3), (cid:98) ψ is continuous, and hence, by (2.5), (3.3), (3.4) and the dominatedconvergence theorem, the periodic function f ψ,n ( x ) in (C.7) is also continuous on [ − π, π ],and hence, bounded. This shows (C.3) and (C.4).We now turn to ( ii ). It suffices to show that, for any 0 < ξ < (cid:12)(cid:12)(cid:12) Ξ jj (cid:48) ,a ( za ) a h − Φ jj (cid:48) ( z ) (cid:12)(cid:12)(cid:12) ≤ Ca min { (cid:36) , (1 − ξ )(2 h +2 β ) } , n → ∞ , since for small enough ξ , conditions (2.7) and (3.3) imply that(1 − ξ )(2 h + 2 β ) > h min + 2 ≥ (cid:36) . In fact, by (A.1), (3.14) and a change of variable y = ax ,Ξ jj (cid:48) ,a ( za ) a h − Φ jj (cid:48) ( z )= 1 a h (cid:90) R e izax | x | − (2 h +1) ρ { g ( x ) − } (cid:98) ψ ( a j x ) (cid:98) ψ ( a j (cid:48) x ) dx (C.9) = (cid:110) (cid:90) | y |≤ a − ξ + (cid:90) | y | >a − ξ (cid:111) e izy | y | − (2 h +1) ρ (cid:110) g (cid:16) ya (cid:17) − (cid:111) (cid:98) ψ (2 j y ) (cid:98) ψ (2 j (cid:48) y ) dy for any 0 < ξ <
1. For large enough n , by (2.6), (2.7) and (3.1)–(3.3), the absolute value ofthe first integral in the sum (C.9) can be bounded by(C.10) Ca (cid:36) (cid:90) | y |≤ a − ξ | y | − (2 h +1)+ (cid:36) | (cid:98) ψ (2 j y ) (cid:98) ψ (2 j (cid:48) y ) | dy ≤ C (cid:48) a (cid:36) . On the other hand, by (3.2) the absolute value of the second term in the sum (C.9) can bebounded by(C.11) C (cid:90) | y | >a − ξ | y | − (2 h +1) − β dy ≤ C (cid:48) a (1 − ξ )(2 h +2 β ) . Expressions (C.10) and (C.11) yield (C.5).
ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 35
We turn to ( iii ). For notational simplicity, consider only two octaves j, j (cid:48) , whence J = j − j + 1 = 2. Then, from (A.29),Υ( n ) = 2( n a,j + n a,j ) . Fix v ∈ C Υ( n ) . For notational simplicity, divide the summation range k = 1 , . . . , Υ( n ) intothe subranges K = { , . . . , n a,j } , K = { n a,j + 1 , . . . , n a,j + n a,j (cid:48) } ,K = { n a,j + n a,j (cid:48) + 1 , . . . , n a,j + n a,j (cid:48) } , K = { n a,j + n a,j (cid:48) + 1 , . . . , n a,j + n a,j (cid:48) ) } . Define the octave and index functions j ( k ) = j { K ∪ K } ( k ) + j (cid:48) { K ∪ K } ( k ) , q ( k ) = 1 1 { K ∪ K } ( k ) + 2 1 { K ∪ K } ( k ) , respectively, which reflects the order of appearance of different octave and index values inthe vector (C.1). By (A.1), v ∗ Γ Y n v = Υ( n ) (cid:88) k =1 Υ( n ) (cid:88) k =1 v k v k (cid:16) Γ Y n (cid:17) k ,k (C.12) = (cid:90) R Υ( n ) (cid:88) k =1 Υ( n ) (cid:88) k =1 v k v k e ia (2 j ( k k − j ( k k ) x (cid:98) ψ ( a j ( k ) x ) (cid:98) ψ ( a j ( k ) x ) (cid:112) (cid:37) (12) ( a j ( k ) ) (cid:37) (12) ( a j ( k ) ) f q ( k ) q ( k ) ( x ) dx By expanding the double summation in (C.12) into each pair in the Cartesian product { K , K , K , K } , we end up with 16 double summation terms under the sign of the inte-gral, i.e., 8 pairs of conjugates. To the cross terms, namely, those involving distinct sum-mation ranges in the index k , we can apply the elementary inequality xy + xy ≤ | x | + | y | .One such pair, associated with the summation regions K × K and K × K , is (cid:90) R (cid:16) (cid:88) k ∈ K v k e ia j k x (cid:98) ψ ( a j x ) (cid:112) (cid:37) (12) ( a j ) (cid:88) k ∈ K v k e − ia j (cid:48) k x (cid:98) ψ ( a j (cid:48) x ) (cid:112) (cid:37) (12) ( a j (cid:48) )+ (cid:88) k ∈ K v k e ia j (cid:48) k x (cid:98) ψ ( a j (cid:48) x ) (cid:112) (cid:37) (12) ( a j (cid:48) ) (cid:88) k ∈ K v k e − ia j k x (cid:98) ψ ( a j x ) (cid:112) (cid:37) (12) ( a j ) (cid:17) f ( x ) dx ≤ (cid:90) R (cid:110)(cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e ia j kx (cid:98) ψ ( a j x ) (cid:112) (cid:37) (12) ( a j ) (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e − ia j (cid:48) kx (cid:98) ψ ( a j (cid:48) x ) (cid:112) (cid:37) (12) ( a j (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:111) f ( x ) dx, since f ( x ) = f ( x ), and analogous bounds hold for the remaining terms. Therefore,(C.12) is bounded by (cid:90) R (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e ia j kx (cid:12)(cid:12)(cid:12) | (cid:98) ψ ( a j x ) | (cid:37) (12) ( a j ) ( f ( x ) + 3 f ( x )) dx + (cid:90) R (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e ia j (cid:48) kx (cid:12)(cid:12)(cid:12) | (cid:98) ψ ( a j (cid:48) x ) | (cid:37) (12) ( a j (cid:48) ) ( f ( x ) + 3 f ( x )) dx + (cid:90) R (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e ia j kx (cid:12)(cid:12)(cid:12) | (cid:98) ψ ( a j x ) | (cid:37) (12) ( a j ) ( f ( x ) + 3 f ( x )) dx + (cid:90) R (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e ia j (cid:48) kx (cid:12)(cid:12)(cid:12) | (cid:98) ψ ( a j (cid:48) x ) | (cid:37) (12) ( a j (cid:48) ) ( f ( x ) + 3 f ( x )) dx. By the change of variable y = a j ( k ) x in each integral, breaking them up (in R ) intosubregions of length 2 π and using the periodicity of Fourier sums, we obtain (cid:90) π − π (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) ∞ (cid:88) l = −∞ (cid:110) ( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:17) +3( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:17)(cid:111) | (cid:98) ψ ( y + 2 πl ) | (cid:37) (12) ( a j ) dy + (cid:90) π − π (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) ∞ (cid:88) l = −∞ (cid:110) ( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:48) (cid:17) +3( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:48) (cid:17)(cid:111) | (cid:98) ψ ( y + 2 πl ) | (cid:37) (12) ( a j (cid:48) ) dy + (cid:90) π − π (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) ∞ (cid:88) l = −∞ (cid:110) ( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:17) +3( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:17)(cid:111) | (cid:98) ψ ( y + 2 πl ) | (cid:37) (12) ( a j ) dy + (cid:90) π − π (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) ∞ (cid:88) l = −∞ (cid:110) ( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:48) (cid:17) +3( a j ) h | y + 2 πl | − (2 h +1) ρ g (cid:16) y + 2 πla j (cid:48) (cid:17)(cid:111) | (cid:98) ψ ( y + 2 πl ) | (cid:37) (12) ( a j (cid:48) ) dy (C.13) ≤ Ca { h ,h }− h ) (cid:90) π − π (cid:110)(cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) (cid:88) k ∈ K v k e iky (cid:12)(cid:12)(cid:12) (cid:111) dy = Ca { h ,h }− h ) v ∗ v for some C >
0. The inequality (C.13) is a consequence of (C.3) (from Lemma C.1, ( i ))and Proposition 3.1, ( iii ), as applied to (cid:37) (12) ( a j ( k ) ). Thus, the claim (C.6) holds. (cid:3) (cid:3) Remark C.1.
For q = 1 and q = 2, let Y n and Y n be as in (A.28) and (C.1), respectively.Then, Y n = P n diag (cid:16) (cid:113) (cid:37) (12) ( a j ) , . . . , (cid:113) (cid:37) (12) ( a j ) (cid:124) (cid:123)(cid:122) (cid:125) n a,j , . . . , (cid:113) (cid:37) (12) ( a j ) , . . . , (cid:113) (cid:37) (12) ( a j ) (cid:124) (cid:123)(cid:122) (cid:125) n a,j , (cid:113) (cid:37) (12) ( a j ) , . . . , (cid:113) (cid:37) (12) ( a j ) (cid:124) (cid:123)(cid:122) (cid:125) n a,j , . . . , (cid:113) (cid:37) (12) ( a j ) , . . . , (cid:113) (cid:37) (12) ( a j ) (cid:124) (cid:123)(cid:122) (cid:125) n a,j (cid:17) Y n =: P n N n Y n ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 37 for some permutation matrix P n . Moreover, since Γ Y n = P n N n Γ Y n N n P ∗ n is a real symmetricmatrix, by Lemma C.1, ( ii ) and ( v ), we obtain the bound(C.14) (cid:107) Γ Y n (cid:107) = (cid:107)P n N n Γ Y n N n P ∗ n (cid:107) ≤ Ca ( n ) { h ,h } for the maximum eigenvalue of Γ Y n , where (cid:107) · (cid:107) is the matrix Euclidean norm.For any n ∈ N , consider the Gaussian vector Y n as in (C.1) but with only one octave j .It will be convenient to reexpress the sum W (12) n ( a j ) as in (3.10) (with q = 1 and q = 2)based on a quadratic form. Define the permutation matrix(C.15) Π n = (cid:18) I n a,j I n a,j (cid:19) ∈ O (2 n a,j ) . Consider the spectral decomposition (C.2). Then, W (12) n ( a j ) = 12 n a,j Y ∗ n Π n Y n (for one octave j )(C.16) d = 1 n a,j Z ∗ ( O Λ / O ∗ )Π n ( O Λ / O ∗ )2 Z d = 1 n a,j Z ∗ Λ / O ∗ Π n O Λ / Z , where Z ∼ N (0 , I n a,j ). Let(C.17) S n = Λ / O ∗ Π n O Λ / , which is a real symmetric matrix. Write out its spectral decomposition(C.18) S n = O S Λ S O ∗ S , O ∈ O (2 n a,j ) , Λ S = diag( λ ,S , . . . , λ n a,j ,S ) . Expression (C.16) becomes W (12) n ( a j ) = 1 n a,j Z ∗ O S Λ S O ∗ S Z d = 1 n a,j Z ∗ Λ S Z (C.19) = (cid:88) i ∈ i + ( n ) λ i,S n a,j Z i + (cid:88) i ∈ i − ( n ) λ i,S n a,j Z i =: (cid:88) i ∈ i + ( n ) η i,n Z i − (cid:88) i ∈ i − ( n ) η i,n Z i , where i + ( n ) and i − ( n ) are the indices for which λ i,S is nonnegative or negative, respectively.In particular, note that Var (cid:104) W (12) n ( a ( n )2 j ) (cid:105) = (cid:107) η n (cid:107) , where(C.20) η n := ( η ,n , . . . , η n a,j ,n ) ∗ and η i,n , i = 1 , . . . , n a,j , are as in expression (C.19).The following lemma is an exponential inequality for χ -like distributions, and corre-sponds to Lemma 1 in [34]. It will be used in the ensuing Lemma C.3 to establish abound on the rate of convergence to zero of the probability P ( W (12) n ( a j ) ≤ r ), where −∞ < r < / Lemma C.2. [34] Let Z , . . . , Z n i.i.d. ∼ N (0 , and η , . . . , η n ≥ , not all zero. Let (cid:107) η (cid:107) and (cid:107) η (cid:107) ∞ be the Euclidean square and sup norms of the vector η = ( η , . . . , η n ) ∗ . Also,define the random variable X = (cid:80) ni =1 η i,n ( Z i − . Then, for every x > , (C.21) P ( X ≥ (cid:107) η (cid:107) √ x + 2 (cid:107) η (cid:107) ∞ x ) ≤ exp( − x ) , (C.22) P ( X ≤ − (cid:107) η (cid:107) √ x ) ≤ exp( − x ) . Lemma C.3.
Let W (12) n ( a j ) be as in (C.19) , and fix −∞ < r < / . Then, for any < ξ < , (C.23) P ( W (12) n ( a ( n )2 j ) ≤ r ) ≤ exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) for large enough n .Proof. Expression (C.19) implies that P ( W (12) n ( a j ) ≤ r ) = P (cid:16) (cid:88) i ∈ i + ( n ) η i,n ( Z i − ≤ − r + (cid:88) i ∈ i − ( n ) η i,n ( Z i − (cid:17) . For notational simplicity, let X n a,j = (cid:80) i ∈ i + ( n ) η i,n ( Z i −
1) and Y n a,j = (cid:80) i ∈ i − ( n ) η i,n ( Z i − r < δ <
1. By the independence of X n a,j and Y n a,j ,(C.24) P ( X n a,j ≤ − r + Y n a,j ) = (cid:110) (cid:90) − δ −∞ + (cid:90) ∞ − δ (cid:111) P ( X n a,j ≤ − r + y ) f Y na,j ( y ) dy, where f Y na,j is the density function of Y n a,j . The first integral in (C.24) is bounded fromabove by P ( X n a,j ≤ − δ + r ) P ( Y n a,j ≤ − δ ). Let η ,n = ( η i,n ) i ∈ i + ( n ) , η ,n = ( η i,n ) i ∈ i − ( n ) . Then, by (C.42) below (under (C.27)),(C.25)max {(cid:107) η ,n (cid:107) , (cid:107) η ,n (cid:107) } ≤ (cid:107) η n (cid:107) ∼ C a h + h ) − h +1 n , (cid:107) η ,n (cid:107) ∞ ≤ (cid:107) η n (cid:107) ∞ ≤ (cid:107) η n (cid:107) , for a constant C ∈ R . Moreover, since − δ + r <
0, then P ( X n a,j ≤ − δ + r ) = P (cid:16) (cid:88) i ∈ i + ( n ) η i,n ( Z i − ≤ − δ + r (cid:17) = P (cid:16) (cid:88) i ∈ i + ( n ) η i,n ( Z i − ≤ − (cid:107) η ,n (cid:107) (cid:115) ( δ − r ) (cid:107) η ,n (cid:107) (cid:17) ≤ exp (cid:110) − ( δ − / (cid:107) η ,n (cid:107) (cid:111) (C.26) ≤ exp (cid:110) − C na h + h ) − h +1 (cid:111) . In (C.26), C does not depend on r and the last two inequalities are a consequence of (C.22)and (C.25) below, respectively, where(C.27) 2 max { h q + h q , h q q } = 2( h q + h q ) , q , q = 1 , . . . , m, ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 39 stems from (2.9). On the other hand, for any 0 < ξ < n the secondintegral in (C.24) is bounded from above by(C.28) (cid:90) ∞ − δ f Y na,j ( y ) dy = P ( Y n a,j > − δ ) = P (cid:16) (cid:88) i ∈ i − ( n ) η i,n ( Z i − > − δ (cid:17) . Suppose, without loss of generality, that (cid:107) η ,n (cid:107) > , n ∈ N (otherwise, P ( (cid:80) i ∈ i − ( n ) η i,n ( Z i − > − δ ) = 0 for n such that (cid:107) η ,n (cid:107) = 0). Therefore,by (C.25), (cid:107) η ,n (cid:107) (cid:16) na h + h ) − h +1 (cid:17) − ξ + (cid:107) η ,n (cid:107) ∞ (cid:16) na h + h ) − h +1 (cid:17) − ξ → , for any 0 < ξ <
1, as n → ∞ . Consequently, for large enough n , (C.28) is bounded by P (cid:16) (cid:88) i ∈ i − ( n ) η i,n ( Z i − > (cid:107) η ,n (cid:107) (cid:16) na h + h ) − h +1 (cid:17) − ξ +2 (cid:107) η ,n (cid:107) ∞ (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:17) (C.29) ≤ exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) , where the last inequality follows from (C.21). From (C.26) and (C.29), since C does notdepend on r , we obtain (C.23). (cid:3) (cid:3) C.2.
Asymptotic covariances.
We will establish Theorem 4.1 at the end of this section,after proving Lemmas C.4–C.7. The latter establish the asymptotic behavior of the first fourmoments of the R -valued random variables W ( q q ) n ( a ( n )2 j ), W ( q q ) n ( a ( n )2 j ). So, considerthe function(C.30) log | S ( q q ) n ( a ( n )2 j ) | {| W ( q q n ( a ( n )2 j ) | >r } , ≤ q ≤ q ≤ m, < r < , where the truncation works as a regularization of the log function around the origin. In theevent | W ( q q ) n ( a ( n )2 j ) | ≤ r , we interpret thatlog | S ( q q ) n ( a ( n )2 j ) | {| W ( q q n ( a ( n )2 j ) | >r } = 0 a.s.Throughout this section, we will make use of the Isserlis theorem (e.g., [48]). For a zeromean, Gaussian random vector Z = ( Z , . . . , Z m ) ∗ ∈ R m , the theorem states that(C.31) E [ Z . . . Z k ] = (cid:88) (cid:89) E [ Z i Z j ] , E [ Z . . . Z k +1 ] = 0 , k = 1 , . . . , (cid:98) m/ (cid:99) . The notation (cid:80) (cid:81) stands for adding over all possible k -fold products of pairs E [ Z i Z j ],where the indices partition the set 1 , . . . , k .The following lemma shows that the high order centered cross moments of W (12) n ( a ( n )2 j ), W (34) n ( a ( n )2 j ) are negligible with respect to their low order counterparts. Lemma C.4.
Let κ , κ ∈ N ∪ { } , κ + κ ≥ , and fix < r < / . Then, as n → ∞ , E (cid:34)(cid:16) W (12) n ( a ( n )2 j ) − a ( n ) δ (cid:17) κ (cid:16) W (34) n ( a ( n )2 j (cid:48) ) − a ( n ) δ (cid:17) κ { min {| W (12) n ( a ( n )2 j (cid:48) ) | , | W (34) n ( a ( n )2 j (cid:48) ) |} >r } (cid:35) (C.32) = O (cid:104)(cid:16) a ( n ) n (cid:17) (cid:105) . Proof.
We first show that(C.33) E (cid:34)(cid:16) W (12) n ( a j ) − a δ (cid:17) κ (cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) κ (cid:35) = O (cid:104)(cid:16) an (cid:17) (cid:105) ,κ , κ ∈ N , κ + κ ≥
3. We will only establish (C.33) for κ = 1 and κ = 2, since theremaining cases can be tackled by a similar argument.The left-hand side of (C.33) can be rewritten as1 a δ +2 δ n a,j n a,j (cid:48) E (cid:104) n a,j (cid:88) k =1 n a,j (cid:48) (cid:88) k =1 n a,j (cid:48) (cid:88) k =1 (cid:16) d ( a j , k ) d ( a j , k ) (cid:37) (12) ( a j ) − (cid:17) (C.34) (cid:16) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:37) (34) ( a j (cid:48) ) − (cid:17)(cid:16) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:37) (34) ( a j (cid:48) ) − (cid:17)(cid:105) , Starting from the assumption (.2), for notational simplicity denote the genericterms under the summation sign in (C.34) as X = d ( a j , k ) / (cid:112) (cid:37) (12) ( a j ), X = d ( a j , k ) / (cid:112) (cid:37) (12) ( a j ), X = d ( a j (cid:48) , k (cid:48) ) / (cid:112) (cid:37) (34) ( a j (cid:48) ), X = d ( a j (cid:48) , k (cid:48) ) / (cid:112) (cid:37) (34) ( a j (cid:48) ), X = d ( a j (cid:48) , k (cid:48) ) / (cid:112) (cid:37) (34) ( a j (cid:48) ), X = d ( a j (cid:48) , k (cid:48) ) / (cid:112) (cid:37) (34) ( a j (cid:48) ). Then, E [( X X − X X − X X − E [ X X X X X X ] − { E [ X X X X ] + E [ X X X X ] + E [ X X X X ] } + 2 . By applying the Isserlis relation (C.31) to the six-fold and four-fold products in (C.35), thelatter expression becomes E [ X X ] E [ X X ] E [ X X ] + E [ X X ] E [ X X ] E [ X X ]+ E [ X X ] E [ X X ] E [ X X ] + E [ X X ] E [ X X ] E [ X X ]+ E [ X X ] E [ X X ] E [ X X ] + E [ X X ] E [ X X ] E [ X X ](C.36) + E [ X X ] E [ X X ] E [ X X ] + E [ X X ] E [ X X ] E [ X X ] . The asymptotic behavior of the summation of each term in the seven-fold sum (C.36) canbe established in the same way, so we only study the first one. Up to a constant, the latteris asymptotically equivalent to a (2 h +2 h +2 h ) − (2 h +4 h ) a δ +2 δ n a,j n a,j (cid:48) n a,j (cid:88) k =1 n a,j (cid:48) (cid:88) k (cid:48) =1 n a,j (cid:48) (cid:88) k (cid:48) =1 E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h E (cid:104) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 41 = a h + h + h ) a h + h +2( h + h ) n a,j n a,j (cid:48) n a,j (cid:48) (cid:88) k (cid:48) =1 n a,j (cid:48) (cid:88) k (cid:48) =1 E (cid:104) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h (C.37) · (cid:110) n a,j (cid:88) k =1 E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h (cid:111) . However, the summation in k in expression (C.37) is bounded by (cid:12)(cid:12)(cid:12) n a,j (cid:88) k =1 E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h (cid:12)(cid:12)(cid:12) ≤ n a,j (cid:88) k =1 (cid:16)(cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:17)(cid:16)(cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:17) ≤ n a,j (cid:88) k =1 (cid:110)(cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) · (cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) (C.38) + (cid:12)(cid:12)(cid:12) Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:111) ≤ C, where C does not depend on k . To justify the last inequality, we only look at the secondterm in (C.38), since the remaining terms can be analyzed in a similar way. It suffices toproceed as in [2], in particular around expression (B.31) in the latter reference. Indeed,starting from (3.2), suppose without loss of generality thatsupp( ψ ) = [0 , . For 0 < ε < /
2, decompose n a,j (cid:88) k =1 (cid:12)(cid:12)(cid:12) Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) = n a,j (cid:88) k =1 (cid:16) (cid:110) max { j, j (cid:48)}| jk − j (cid:48) k (cid:48) | >ε (cid:111) + 1 (cid:110) max { j, j (cid:48)}| jk − j (cid:48) k (cid:48) | ≤ ε (cid:111)(cid:17) (C.39) (cid:12)(cid:12)(cid:12) Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12) . The first sum term in (C.39) only contains a finite number of terms, where such numberdoes not depend on k or k (cid:48) . Moreover,max (cid:110) | Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) | , (cid:12)(cid:12)(cid:12) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) − Φ jj (cid:48) ( a (2 j k − j (cid:48) k (cid:48) )) a h (cid:12)(cid:12)(cid:12)(cid:111) ≤ C, where C does not depend k , k (cid:48) or k (cid:48)(cid:48) . The latter statement follows from the fact that Φ jj (cid:48) is, up to a change of sign, the wavelet variance of a fBm and from (A.24). Moreover, byProposition II.2 in [13] and again by (A.24), the second sum term in (C.39) is bounded by C (cid:80) z ∈ Z \{ } z − . This shows that (C.39) is bounded by a constant not depending on k , k (cid:48) or k (cid:48) . Hence, (C.38) holds.Consequently, up to a constant the absolute value of (C.37) is bounded from above by C (cid:48) a h + h + h )+2 a h + h +2( h + h ) n n a,j (cid:48) n a,j (cid:48) (cid:88) k (cid:48) =1 n a,j (cid:48) (cid:88) k (cid:48) =1 (cid:12)(cid:12)(cid:12) E (cid:104) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h (cid:12)(cid:12)(cid:12) ≤ C (cid:48)(cid:48) (cid:16) an (cid:17) , where we used the fact that a δ δ ≤
1. This establishes (C.33).So, rewrite E (cid:34)(cid:16) W (12) n ( a j ) − a δ (cid:17)(cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) (cid:35) − E (cid:34)(cid:16) W (12) n ( a j ) − a δ (cid:17) { W (12) n ( a j ) >r } (cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) { W (34) n ( a j (cid:48) ) >r } (cid:35) = E (cid:104)(cid:16) W (12) n ( a j ) − a δ (cid:17)(cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) (cid:110) { W (12) n ( a j ) ≤ r } { W (34) n ( a j (cid:48) ) >r } (C.40) +1 { W (12) n ( a j ) >r } { W (34) n ( a j (cid:48) ) ≤ r } + 1 { W (12) n ( a j ) ≤ r } { W (34) n ( a j (cid:48) ) ≤ r } (cid:111)(cid:105) . The asymptotic behavior of every term on the right-hand side of (C.40) can be establishedin the same way, so we only look at the first one. For any 0 < ξ <
1, the Cauchy-Schwarzinequality, Lemma C.3 and expression (C.33) yield (cid:12)(cid:12)(cid:12) E (cid:34)(cid:16) W (12) n ( a j ) − a δ (cid:17) { W (12) n ( a j ) ≤ r } (cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) { W (34) n ( a j (cid:48) ) >r } (cid:35) (cid:12)(cid:12)(cid:12) ≤ (cid:118)(cid:117)(cid:117)(cid:116) E (cid:34)(cid:16) W (12) n ( a j ) − a δ (cid:17) (cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) (cid:35)(cid:113) P ( W (12) n ( a j ) ≤ r ) ≤ C (cid:16) an (cid:17) exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) , where C does not depend on r . Moreover, under (3.6),exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) = o (cid:16) an (cid:17) . ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 43
Therefore, (C.32) holds. (cid:3) (cid:3)
The following lemma expresses, up to a residual term, the cross-covariance (first cross-moment) between sample wavelet variances in terms of the functions (3.14).
Lemma C.5.
Let Φ jj (cid:48) ·· ( z ) and n ∗ be as in (3.14) and (4.1) , respectively. For < r < / , E (cid:104) ( W (12) n ( a ( n )2 j ) − {| W (12) n ( a ( n )2 j ) | >r } ( W (34) n ( a ( n )2 j (cid:48) ) − {| W (34) n ( a ( n )2 j (cid:48) ) | >r } (cid:105) = 2 − ( j + j (cid:48) ) Φ jj (0)Φ j (cid:48) j (cid:48) (0) { O ( a ( n ) − (cid:36) ) } (cid:110) a ( n ) h + h ) − h + h ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:105) + a ( n ) h + h ) − h + h ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:105) (C.41) + o (cid:16) a { h + h ,h + h }− h + h ) n ∗ (cid:17)(cid:111) , as n → ∞ .Proof. As in the proof of Lemma C.4, we first drop the indicator functions on the left-handside of (C.41) and investigate the limit. We will show thatCov (cid:104) W (12) n ( a j ) , W (34) n ( a j (cid:48) ) (cid:105) = 2 − ( j + j (cid:48) ) [Φ jj (0)Φ j (cid:48) j (cid:48) (0)(1 + O ( a − (cid:36) )) ] · (cid:110) a h + h ) − h + h ) n ∗ n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )+ a h + h ) − h + h ) n ∗ n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )(C.42) + o (cid:16) a { h + h ,h + h }− h + h ) n ∗ (cid:17)(cid:111) . In fact, the left-hand side of (C.46) can be written as E (cid:104) ( W (12) n ( a j ) − W (34) n ( a j (cid:48) ) − (cid:105) (C.43) = 1 n a,j n a,j (cid:48) n a,j (cid:88) k =1 n a,j (cid:48) (cid:88) k (cid:48) =1 (cid:110) E (cid:34) d ( a j , k ) d ( a j , k ) (cid:37) (12) ( a j ) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:37) (34) ( a j (cid:48) ) (cid:35) − (cid:111) . By the Isserlis relation (C.31), the first term in the argument of the sum (C.43) can bereexpressed as E (cid:104) d ( a j , k ) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) d ( a j (cid:48) , k (cid:48) ) (cid:105) (cid:37) (12) ( a j ) (cid:37) (34) ( a j (cid:48) ) = (cid:110) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) (cid:37) (12) ( a j ) (cid:37) (34) ( a j (cid:48) )(C.44) + E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) (cid:37) (12) ( a j ) (cid:37) (34) ( a j (cid:48) ) (cid:111) . By Lemma C.1, ( ii ), and (C.44), we can rewrite (C.43) as1 n a,j n a,j (cid:48) n a,j (cid:88) k =1 n a,j (cid:48) (cid:88) k (cid:48) =1 (cid:110) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) (cid:37) (12) ( a j ) (cid:37) (34) ( a j (cid:48) )+ E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) (cid:37) (12) ( a j ) (cid:37) (34) ( a j (cid:48) ) (cid:111) = 2 − ( j + j (cid:48) ) Φ jj (0)Φ j (cid:48) j (cid:48) (0)(1 + O ( a − (cid:36) )) (cid:110) a h + h ) − h + h ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h (cid:105) (C.45)+ a h + h ) − h + h ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h (cid:105)(cid:111) . By adding and subtracting the counterparts Φ jj (cid:48) ·· (2 j k − j (cid:48) k (cid:48) ) for each term, up to the factor2 − ( j + j (cid:48) ) / { Φ jj (0)Φ j (cid:48) j (cid:48) (0)(1 + O ( a − (cid:36) )) } the expression (C.45) can be written as a h + h ) − h + h ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17)(cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) +Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) +Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) +Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:111)(cid:105) + a h + h ) − h + h ) n ∗ (cid:104) n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 (cid:110)(cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) · (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ j (cid:48) j (2 j k − j (cid:48) k (cid:48) ) (cid:17) ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 45 +Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) +Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17) +Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:111)(cid:105) = a h + h ) − h + h ) n ∗ n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )+ a h + h ) − h + h ) n ∗ n ∗ j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) )(C.46) + o (cid:16) a { h + h ,h + h }− h + h ) n ∗ (cid:17) . It remains to justify the order of the error term in (C.46). So, by adapting the proof of(A.16), we obtain a h + h ) − h + h ) n ∗ n ∗ (cid:12)(cid:12)(cid:12) j n ∗ (cid:88) k =1 2 j (cid:48) n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17)(cid:12)(cid:12)(cid:12) = a h + h ) − h + h ) n ∗ o (1) = o (cid:16) a h + h ) − h + h ) n ∗ (cid:17) . The same bound holds for a h + h ) − h + h ) n ∗ n ∗ (cid:12)(cid:12)(cid:12) j n ∗ (cid:88) k =1 2 j (cid:48) n ∗ (cid:88) k (cid:48) =1 Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17)(cid:12)(cid:12)(cid:12) . and a h + h ) − h + h ) n ∗ n ∗ (cid:12)(cid:12)(cid:12) j (cid:48) n ∗ (cid:88) k =1 2 j n ∗ (cid:88) k (cid:48) =1 (cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17)(cid:16) E (cid:104) d ( a j , k ) d ( a j (cid:48) , k (cid:48) ) (cid:105) a h − Φ jj (cid:48) (2 j k − j (cid:48) k (cid:48) ) (cid:17)(cid:12)(cid:12)(cid:12) By extending this analysis to the remaining terms of (C.46), we obtain analogous boundsand the error term o (cid:16) a { h h ,h h }− h h n ∗ (cid:17) , as claimed.In view of (C.45) and (C.46), it suffices to show that the indicator functions on theleft-hand side of (C.41) do not affect the approximation order. In fact, E (cid:104) ( W (12) n ( a j ) − W (34) n ( a j (cid:48) ) − (cid:105) − E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) >r } ( W (34) n ( a j (cid:48) ) − { W (34) n ( a j (cid:48) ) >r } (cid:105) = E (cid:104) ( W (12) n ( a j ) − W (34) n ( a j (cid:48) ) − · (cid:16) { W (12) n ( a j ) >r } { W (34) n ( a j (cid:48) ) ≤ r } (C.47) +1 { W (12) n ( a j ) ≤ r } { W (34) n ( a j (cid:48) ) >r } + 1 { W (12) n ( a j ) ≤ r } { W (34) n ( a j (cid:48) ) ≤ r } (cid:17)(cid:105) . Define(C.48) h ∗ = max q ,q =1 , , , h q q , h ∗ = min q ,q =1 , , , h q q . For 0 < ξ (cid:48) <
1, by the Cauchy-Schwarz inequality, expression (C.33) (from Lemma C.4)and Lemma C.3, the first term on the right-hand side of (C.47) is bounded by (cid:12)(cid:12)(cid:12) E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) >r } ( W (34) n ( a j (cid:48) ) − { W (34) n ( a j (cid:48) ) ≤ r } (cid:105) (cid:12)(cid:12)(cid:12) ≤ a δ + δ (cid:118)(cid:117)(cid:117)(cid:116) E (cid:34)(cid:16) W (12) n ( a j ) − a δ (cid:17) (cid:16) W (34) n ( a j (cid:48) ) − a δ (cid:17) (cid:35)(cid:113) P ( W (34) n ( a j ) ≤ r ) ≤ C (cid:16) a h max − h min +1 n (cid:17) exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:48) (cid:111) = o (cid:16) a { h + h ,h + h }− h + h ) n ∗ (cid:17) , where the constant C does not depend on r and the last equality is a consequence of (3.6).Similar bounds hold for the remaining terms on the right-hand side of (C.47). Therefore,the expression (C.41) follows. (cid:3) (cid:3) The following lemma describes the decay rate of the first individual truncated momentof the wavelet variance.
Lemma C.6.
For any < ξ < and < r < / , max (cid:110)(cid:12)(cid:12)(cid:12) E (cid:104) { W (12) n ( a ( n )2 j ) − } { W (12) n ( a ( n )2 j ) ≤ r } (cid:105) (cid:12)(cid:12)(cid:12) ; (cid:12)(cid:12)(cid:12) E (cid:104) { W (12) n ( a ( n )2 j ) − } { W (12) n ( a ( n )2 j ) >r } (cid:105) (cid:12)(cid:12)(cid:12)(cid:111) (C.49) = O (cid:16) a ( n ) h + h − h +1 / √ n exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111)(cid:17) . Proof.
Notice that0 = E (cid:104) W (12) n ( a j ) − (cid:105) = E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) >r } (cid:105) + E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) ≤ r } (cid:105) . Hence, E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) >r } (cid:105) = − E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) ≤ r } (cid:105) . However, by the Cauchy-Schwarz inequality,(C.50) E (cid:104) ( W (12) n ( a j ) − { W (12) n ( a j ) ≤ r } (cid:105) (cid:12)(cid:12)(cid:12) ≤ (cid:114) Var (cid:104) W (12) n ( a j ) (cid:105)(cid:113) P ( W (12) n ( a j ) ≤ r ) . By expressions (C.50), (C.42) and Lemma C.3, the expression (C.49) follows. (cid:3) (cid:3)
The following lemma establishes the decay of the covariances between truncated terms(C.30) and indicators involving wavelet variances, or between indicators only.
Lemma C.7.
For any < ξ < and < r < / , ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 47 ( i ) Cov (cid:104) {| W (12) n ( a j (cid:48) ) | >r } , {| W (34) n ( a j (cid:48) ) | >r } (cid:105) (C.51) ≤ exp (cid:110) − (cid:104)(cid:16) na h + h ) − h +1 (cid:17) − ξ + (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:105)(cid:111) ;( ii ) Var (cid:104) log | W (12) n ( a ( n )2 j ) | { W (12) n ( a ( n )2 j ) < − r } (cid:105) (C.52) ≤ ( C + log ( r )) exp (cid:110) − (cid:16) na ( n ) h + h ) − h +1 (cid:17) − ξ (cid:111) and Var (cid:104) log W (12) n ( a ( n )2 j )1 { W (12) n ( a ( n )2 j ) >r } (cid:105) (C.53) ≤ log ( r ) exp (cid:110) − (cid:16) na ( n ) h + h ) − h +1 (cid:17) − ξ (cid:111) + o (1);( iii ) Cov (cid:104) log | W (12) n ( a j ) | {| W (12) n ( a j ) | >r } , {| W (34) n ( a j (cid:48) ) | >r } (cid:105) ≤ (cid:115) log ( r ) exp (cid:110) − (cid:16) na ( n ) h + h ) − h +1 (cid:17) − ξ + o (1)(C.54) · exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) . In (C.52) , C > does not depend on r .Proof. We first show (C.51). By the Cauchy-Schwarz inequality, the left-hand side of (C.51)is bounded from above by(C.55) (cid:114)
Var (cid:104) {| W (12) n ( a j ) | >r } (cid:105)(cid:114) Var (cid:104) {| W (34) n ( a j (cid:48) ) | >r } (cid:105) . Moreover, for 0 < ξ < j , Lemma C.3 implies thatVar (cid:104) {| W (12) n ( a j ) | >r } (cid:105) = P ( | W (12) n ( a j ) | > r ) P ( | W (12) n ( a j ) | ≤ r )(C.56) ≤ exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) . The same bound holds for the octave j (cid:48) in (C.55). Thus, (C.51) holds.To prove (C.52), note that Var (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) < − r } (cid:105) is bounded by E (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) < − r } (cid:105) = E (cid:104) log | W (12) n ( a j ) | (cid:16) { W (12) n ( a j ) ≤− / } + 1 {− / 2) + log ( r ) P ( − / < W (12) n ( a j ) < − r ) ≤ C (cid:48) exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) + log ( r ) exp (cid:110) − (cid:16) na h + h ) − h +1 (cid:17) − ξ (cid:111) . This establishes (C.52).To prove (C.53), note that Var (cid:104) log W (12) n ( a j )1 { W (12) n ( a j ) >r } (cid:105) is bounded by E (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) >r } (cid:105) = E (cid:104) log | W (12) n ( a j ) | (cid:16) { r 2) + E (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) ≥ / } (cid:105) . However, for some C > | W (12) n ( a j ) | { W (12) n ( a j ) ≥ / } ≤ CW (12) n ( a j ) , where, by (A.39), W (12) n ( a j ) P → E (cid:104) W (12) n ( a j ) (cid:105) = Var W (12) n ( a j ) + 1 → , n → ∞ . Therefore, by the dominated convergence theorem for convergence in probability,lim n →∞ E (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) ≥ / } (cid:105) = 0 . This establishes (C.53).To show (C.54), again by applying the Cauchy-Schwarz inequality, the left-hand side of(C.54) is bounded from above by(C.57) (cid:114) Var (cid:104) log | W (12) n ( a j ) | {| W (12) n ( a j ) | >r } (cid:105)(cid:114) Var (cid:104) {| W (34) n ( a j (cid:48) ) | >r } (cid:105) . Hence, the bound follows from (C.51) and (C.53). (cid:3) (cid:3) We are now in a position to establish Theorem 4.1. Proof of Theorem 4.1 : Fix 0 < ξ < n ∗ = na j + j (cid:48) . Then,Cov (cid:104) log | S (12) n ( a j ) | {| W (12) n ( a j ) | >r n } , log | S (34) n ( a j (cid:48) ) | {| W (34) n ( a j (cid:48) ) | >r n } (cid:105) = Cov (cid:104) log | W (12) n ( a j ) | {| W (12) n ( a j ) | >r n } , log | W (34) n ( a j (cid:48) ) | {| W (34) n ( a j (cid:48) ) | >r n } (cid:105) + log (cid:12)(cid:12) E (cid:104) d ( a j , d ( a j (cid:48) , (cid:105) (cid:12)(cid:12) Cov (cid:104) log | W (12) n ( a j ) | {| W (12) n ( a j ) | >r n } , {| W (34) n ( a j (cid:48) ) | >r n } (cid:105) + log (cid:12)(cid:12) E (cid:2) d ( a j , d ( a j , (cid:3) (cid:12)(cid:12) Cov (cid:104) {| W (12) n ( a j ) | >r n } , log | W (34) n ( a j (cid:48) ) | {| W (34) n ( a j (cid:48) ) | >r n } (cid:105) + log (cid:12)(cid:12) E (cid:2) d ( a j , d ( a j , (cid:3) (cid:12)(cid:12) log (cid:12)(cid:12) E (cid:104) d ( a j (cid:48) , d ( a j (cid:48) , (cid:105) (cid:12)(cid:12) Cov (cid:104) {| W (12) n ( a j ) | >r n } , {| W (34) n ( a j (cid:48) ) | >r n } (cid:105) = Cov (cid:104) log | W (12) n ( a j ) | {| W (12) n ( a j ) | >r n } , log | W (34) n ( a ( n )2 j (cid:48) ) | {| W (34) n ( a ( n )2 j (cid:48) ) | >r n } (cid:105) (C.58) + o (cid:16)(cid:16) a ( n ) h max − h min n ∗ (cid:17) (cid:17) . ULTIVARIATE HADAMARD SELF-SIMILARITY: TESTING FRACTAL CONNECTIVITY 49 The last two equalities in (C.58) are a consequence of Lemma C.1, ( ii ), as applied to E (cid:2) d ( a j , d ( a j , (cid:3) and E (cid:104) d ( a j (cid:48) , d ( a j (cid:48) , (cid:105) , and of the bound (C.54) (from LemmaC.7, ( iii )) under the condition (4.2).Therefore, it suffices to show that the main term on the right-hand side of (C.58) is equalto the main term on the right-hand side of (4.3). By accounting for absolute values, thecovariance term in the former can be broken up into a sum of four terms, namely,(C.59) Cov (cid:104) log W (12) n ( a j )1 { W (12) n ( a j ) >r n } , log W (34) n ( a j (cid:48) )1 { W (34) n ( a j (cid:48) ) >r n } (cid:105) plus the remainderCov (cid:104) log W (12) n ( a j ) 1 { W (12) n ( a j ) >r n } , log | W (34) n ( a j ) | { W (34) n ( a j ) < − r n } (cid:105) +Cov (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) < − r n } , log W (34) n ( a j ) 1 { W (34) n ( a j ) >r n } (cid:105) (C.60) +Cov (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) < − r n } , log | W (34) n ( a j ) | { W (34) n ( a j ) < − r n } (cid:105) By the Cauchy-Schwarz inequality, the bounds (C.52) and (C.53) (from Lemma C.7, ( ii ))and condition (4.2), the absolute value of the second term in the sum (C.60) is bounded by (cid:114) Var (cid:104) log | W (12) n ( a j ) | { W (12) n ( a j ) < − r n } (cid:105) Var (cid:104) log W (34) n ( a j )1 { W (34) n ( a j ) >r n } (cid:105) = o (cid:16)(cid:16) a ( n ) h max − h min n ∗ (cid:17) (cid:17) . By a similar argument, the same bound holds for the remaining terms in the sum (C.60).Thus, it suffices to focus on (C.59). In the following developments, expressions involvingindividual sample wavelet variance terms will be expressed in terms of W (12) n ( a j ), butanalogous expressions hold when substituting W (34) n ( a j (cid:48) ) for W (12) n ( a j ).Fix 0 < ξ < 1. For any given j , write out the almost sure Taylor expansionlog W (12) n ( a j ) 1 { W (12) n ( a j ) >r n } (C.61) = (cid:110) ( W (12) n ( a j ) − − (cid:16) W (12) n ( a j ) − θ ( W (12) n ( a j )) (cid:17) (cid:111) { W (12) n ( a j ) >r n } , where θ + ( W (12) n ( a j )) ∈ [min { W (12) n ( a j ) , } , max { W (12) n ( a j ) , } ]. Then, E (cid:104) log W (12) n ( a j ) log W (34) n ( a j (cid:48) ) 1 min {{ W (12) n ( a j ) ,W (34) n ( a j (cid:48) ) } >r n } (cid:105) = E (cid:104) ( W (12) n ( a j ) − W (34) n ( a j (cid:48) ) − { min { W (12) n ( a j ) ,W (34) n ( a j (cid:48) ) } >r n } (cid:105) − E (cid:34) ( W (12) n ( a j ) − (cid:16) W (34) n ( a j (cid:48) ) − θ + ( W (34) n ( a j (cid:48) )) (cid:17) { min { W (12) n ( a j ) ,W (34) n ( a j (cid:48) ) } >r n } (cid:35) − E (cid:34)(cid:16) W (12) n ( a j ) − θ + ( W (12) n ( a j )) (cid:17) ( W (34) n ( a j (cid:48) ) − 1) 1 { min { W (12) n ( a j ) ,W (34) n ( a j (cid:48) ) } >r n } (cid:35) (C.62) + 14 E (cid:34)(cid:16) W (12) n ( a j ) − θ + ( W (12) n ( a j )) (cid:17) (cid:16) W (34) n ( a j (cid:48) ) − θ + ( W (34) n ( a j (cid:48) )) (cid:17) { min { W (12) n ( a j ) ,W (34) n ( a j (cid:48) ) } >r n } (cid:35) . For 0 < r n < / 2, recast (cid:16) W (12) n ( a j ) − (cid:98) θ + ( W (12) n ( a j )) (cid:17) { W (12) n ( a j ) >r n } = (cid:16) W (12) n ( a j ) − (cid:98) θ + ( W (12) n ( a j )) (cid:17) (cid:16) { r n Remark C.2.