Patrice Abry

B Cooper Boniece

Gustavo Didier

Herwig Wendt

Wavelet eigenvalue regression in high dimensions * † ‡ §

Keywords: wavelets, operator self-similarity, random matrices

In this paper, we construct the wavelet eigenvalue regression methodology (Abry and Didier (2018b, 2018a)) in high dimensions. We assume that possibly non-Gaussian, finite-variance p-variate measurements are made of a low-dimensional r-variate (r p) fractional stochastic process with non-canonical scaling coordinates and in the presence of additive high-dimensional noise. The measurements are correlated both time-wise and between rows. Building upon the asymptotic and large scale properties of wavelet random matrices in high dimensions, the wavelet eigenvalue regression is shown to be consistent and, under additional assumptions, asymptotically Gaussian in the estimation of the fractal structure of the system. We further construct a consistent estimator of the effective dimension r of the system that significantly increases the robustness of the methodology. The estimation performance over finite samples is studied by means of simulations.

Introduction

A wavelet is a unit L 2 (R)-norm function that annihilates polynomials. For a fixed (octave) j ∈ N, a wavelet random matrix is given by

W(a(n)2 j ) = 1 n a,j n a,j k=1 D(a(n)2 j , k)D(a(n)2 j , k) * .
(1.1)

In (1.1), * denotes transposition, n a,j = n/a(n)2 j is the number of wavelet-domain observations for a sample size n, and each random vector D(a(n)2 j , k) is the wavelet transform of a multivariate stochastic process Y at the dyadic scale a(n)2 j and shift k ∈ Z. The entries of {D(a(n)2 j , k)} k∈N are generally correlated. The so-named wavelet eigenanalysis methodology consists in using the behavior across scales of the eigenvalues of wavelet random matrices to study the fractality of stochastic systems (Abry andDidier (2018b, 2018a)). In this paper, we build upon recent results on the properties of large wavelet random matrices [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]) to construct a wavelet eigenvalue regression methodology in high dimensions. The (possibly non-Gaussian) underlying stochastic process Y = {Y (t)} t∈Z is assumed to have the form Y (t) = P X(t) + Z(t), t ∈ Z.

(1.2)

In (1.2), both Y and the noise term Z = {Z(t)} t∈Z are (high-dimensional) p = p(n)-variate processes, P = P (n) is a rectangular coordinates matrix and, for fixed r ∈ N, X = {X(t)} t∈Z is a (low-dimensional) r-variate fractional process. In particular, the measurements Y are correlated time-wise and between rows. We show that, if the ratio a(n)p(n)/n converges to a positive constant, the wavelet eigenvalue regression provides a consistent and, under additional assumptions, asymptotically Gaussian estimator of the underlying low-dimensional fractal structure of the system. In addition, we construct a consistent estimator of the effective dimension r of the system that significantly increases the robustness of the statistical methodology. The performance of the statistical protocols over finite samples is further studied by means of simulations.

Since the 1950s, the spectral behavior of large dimensional random matrices has attracted considerable attention from the mathematical research community. In quantum mechanics, for example, random matrices are of great interest as statistical mechanical models of infinite dimensional and possibly unknown Hamiltonian operators (e.g., [START_REF] Mehta | On the density of eigenvalues of a random matrix[END_REF], [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF], [START_REF] Arous | Large deviations for Wigner's law and Voiculescu's noncommutative entropy[END_REF], [START_REF] Soshnikov | Universality at the edge of the spectrum in Wigner random matrices[END_REF], [START_REF] Mehta | Random Matrices[END_REF], [START_REF] Deift | Universality for mathematical and physical systems[END_REF], [START_REF] Anderson | An Introduction to Random Matrices[END_REF], [START_REF] Erdős | Rigidity of eigenvalues of generalized Wigner matrices[END_REF]). Random matrices have also naturally emerged as one essential mathematical framework for the modern era of "Big Data" [START_REF] Briody | Big data: Harnessing a game-changing asset[END_REF]), when hundreds to several tens of thousands of time series get recorded and stored on a daily basis. In coping with the data deluge, one is often interested in understanding the behavior of random constructs such as the spectral distribution of sample covariance matrices when the dimension p is comparable to the sample size n, including instances with dependent measurements (e.g., [START_REF] Tao | Random covariance matrices: universality of local statistics of eigenvalues[END_REF], [START_REF] Paul | Random matrix theory in statistics: a review[END_REF], [START_REF] Basu | Regularized estimation in sparse high-dimensional time series models[END_REF], [START_REF] Giraud | Introduction to High-Dimensional Statistics[END_REF], [START_REF] Yao | Large Sample Covariance Matrices and High-Dimensional Data Analysis[END_REF], [START_REF] Chakrabarty | From random matrices to long range dependence[END_REF], [START_REF] Merlevède | On the empirical spectral distribution for matrices with long memory and independent rows[END_REF], [START_REF] Che | Universality of random matrices with correlated entries[END_REF], [START_REF] Taylor | On partitioning multivariate self-affine time series[END_REF], [START_REF] Wang | Spectral analysis of sample autocovariance matrices of a class of linear time series in moderately high dimensions[END_REF], [START_REF] Zhang | Gaussian approximation for high dimensional time series[END_REF], [START_REF] Erdős | Random matrices with slow correlation decay[END_REF], [START_REF] Horváth | Asymptotics for empirical eigenvalue processes in high-dimensional linear factor models[END_REF], [START_REF] Merlevède | Unbounded largest eigenvalue of large sample covariance matrices: Asymptotics, fluctuations and applications[END_REF], [START_REF] Wainwright | High-Dimensional Statistics: a Non-Asymptotic Viewpoint[END_REF], [START_REF] Bourguin | Limiting behavior of large correlated Wishart matrices with chaotic entries[END_REF]).

In turn, scale invariance manifests itself in a wide range of natural and social phenomena such as in climate studies [START_REF] Isotta | The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine raingauge data[END_REF]), critical phenomena [START_REF] Sornette | Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforganization and Disorder: Concepts and Tools[END_REF]), dendrochronology [START_REF] Bai | How the instability of ranks under long memory affects largesample inference[END_REF]), hydrology [START_REF] Benson | Wavelet-based detection and estimation of fractional Lévy signals in high dimensions[END_REF]) and turbulence [START_REF] Kolmogorov | The local structure of turbulence in an incompressible fluid at very high Reynolds numbers[END_REF]). In a multidimensional setting, scaling behavior does not always appear along standard coordinate axes, and often involves multiple (scaling) relations. A R r -valued stochastic process X is called operator self-similar (o.s.s.; [START_REF] Laha | Operator self similar stochastic processes in rd[END_REF], [START_REF] Hudson | Operator-self-similar processes in a finite-dimensional space[END_REF]) if it exhibits the scaling property {X(ct)} t∈R f.d.d.

= {c H X(t)} t∈R , c > 0.

(1.3)

In (1.3), H is some (Hurst) matrix H whose eigenvalues have real parts lying in the interval (0, 1] and c H := exp{log(c)H} = ∞ k=0 (log(c)H) k k!

. A canonical model for multivariate fractional systems is operator fractional Brownian motion (ofBm), namely, a Gaussian, o.s.s., stationaryincrement stochastic process [START_REF] Maejima | Operator-self-similar stable processes[END_REF], [START_REF] Mason | Sample path properties of operator-self-similiar Gaussian random fields[END_REF], [START_REF] Didier | Exponents, symmetry groups and classification of operator fractional Brownian motions[END_REF]). In particular, ofBm is the natural multivariate generalization of the classical fBm [START_REF] Mandelbrot | Fractional Brownian motions, fractional noises and applications[END_REF]).

In the characterization of scaling properties, the use of eigenanalysis was first proposed in Meerschaert andScheffler (1999, 2003) and [START_REF] Becker-Kern | Parameter estimation of selfsimilarity exponents[END_REF]. It has also been used in the cointegration literature (e.g., [START_REF] Phillips | Testing for cointegration using principal components methods[END_REF], [START_REF] Li | On determination of cointegration ranks[END_REF]). In Abry andDidier (2018a, 2018b), wavelet eigenanalysis is proposed in the construction of a general methodology for the statistical identification of the scaling (Hurst) structure of ofBm in low dimensions.

Wavelet random matrices were used in high-dimensional modeling contexts first in Abry et al. (2018) and Boniece et al. (2019). In [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], the fundamental mathematical properties of the eigenvalues λ 1 W(a(n)2 j ) ≤ . . . ≤ λ p(n) W(a(n)2 j ) of large wavelet random matrices are, to the best of our knowledge, for the first time considered. In this paper, we build upon such fundamental properties to construct robust statistical methodology in high dimensions. The measurements are assumed to be of the form (1.2), where the fractional behavior of X is characterized by a scaling matrix of the Jordan form

H = P H diag(h 1 , . . . , h r ) P -1 H , (1.4)
with real eigenvalues. The measurements Y display correlation time-wise and between rows. In applications, models of the form (1.2) and related models appear, for example, in neuroscience, fMRI imaging and signal processing [START_REF] Ciuciu | Scale-free and multifractal properties of fMRI signals during rest and task[END_REF], [START_REF] Liu | On the Marčenko-Pastur law for linear time series[END_REF]cf. Chauduri et al. (2019), [START_REF] Stringer | Highdimensional geometry of population responses in visual cortex[END_REF]) and in econometrics (e.g., [START_REF] Brown | The number of factors in security returns[END_REF], [START_REF] Zhang | Identifying cointegration by eigenanalysis[END_REF]). The goal of statistical inference is to characterize the low-dimensional (r-variate) fractal structure in high-dimensional data. This requires robustness not only with respect to non-canonical scaling coordinates, as in fixed dimensions (the so-called amplitude and dominance effects; see [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]), but also to high-dimensional environmental noise. Starting from measurements Y , we propose a multiscale wavelet eigenvalue regression estimator

{ i } i=1,...,p := 1 2 j 2 j=j 1 w j log 2 λ i W(a(n)2 j ) -1 i=1,...,p
, for fixed wavelet octaves j 1 , j 2 and particular weights w j (see Definition 3.1). We show that the fixed-dimensional subvector

{ h q } q=1,...,r := { i } i=p(n)-r+1,...,p(n) (1.5)
is a consistent and, under additional assumptions, asymptotically Gaussian vector estimator of the scaling eigenvalue structure {h q } q=1,...,r of X (Theorems 3.1 and 3.2; see also Figure 1). The asymptotic properties of the estimator are grounded in the high-dimensional behavior of the eigenvalues of wavelet random matrices. In fact, as shown in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], as n and a(n) grow the r largest wavelet eigenvalues satisfy scaling relations of the form

λ p(n)-r+q W(a(n)2 j ) P ∼ ξ q (2 j ) a(n) 2hq+1 , q = 1, . . . , r, (1.6)
for deterministic ξ q (2 j ). Also, under additional assumptions, they display asymptotically Gaussian fluctuations (cf. Theorems B.1 and B.2). The high-dimensional statistical analysis requires taking the three-way limit

lim n→∞ p(n)a(n) n =: c ∈ [0, ∞). (1.7)
For the sake of illustration and comparison, note that traditional analysis of large sample covariance matrices usually involves the ratio lim n→∞ p(n)/n and the largest eigenvalue often exhibits universality in the form of Tracy-Widom fluctuations [START_REF] Bai | Spectral Analysis of Large Dimensional Random Matrices[END_REF], [START_REF] Lee | Tracy-Widom distribution for the largest eigenvalue of real sample covariance matrices with general population[END_REF]). Vis-à-vis [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], the main mathematical contribution of this paper stems from the fact that the wavelet eigenvalue regression demands careful analysis of the rate of convergence of the rescaled r largest eigenvalues of the deterministic wavelet matrices EW(a(n)2 j ). This involves studying the rate of angular convergence of the associated (deterministic) wavelet eigenvectors in terms of the high-dimensional coordinates P = P (n) (see Proposition C.1). Since (1.5) presupposes knowledge of the dimension r of the hidden fractional process X, we further put forth a wavelet eigenvalue regression-based estimator r of r (see Definition 3.2 and Theorem 3.3). The construction is based on the fact that, by contrast with (1.6), the lower eigenvalues λ W(a(n)2 j ) , = 1, . . . , p(n) -r, of wavelet random matrices are bounded in probability (cf. Theorem B.1). The use of the estimator r greatly increases the robustness of the overall statistical methodology over finite samples by providing a quantitative measure of separation between the scaling behaviors of large and small wavelet log-eigenvalues (cf. Figure 1; for effective dimension estimation in different contexts, see, for example, [START_REF] Nadakuditi | Sample eigenvalue based detection of high-dimensional signals in white noise using relatively few samples[END_REF], [START_REF] Little | Estimation of intrinsic dimensionality of samples from noisy low-dimensional manifolds in high dimensions with multiscale SVD[END_REF], [START_REF] Lam | Factor modeling for high-dimensional time series: inference for the number of factors[END_REF]).

As in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], for the sake of clarity and mathematical generality, the assumptions are stated directly in the wavelet domain (Section 2). In particular, the measurements Y are possibly non-Gaussian. A brief illustration of instances covered by the assumptions is provided in Examples 2.1-2.3.

We provide broad computational studies that demonstrate the convergence to Gaussianity of { h q } q=1,...,r in the high dimensional limit for various instances of n a,j , p and r. The experiments confirm that convergence takes place regardless of the (fixed) ratio p/n a,j . We further use computational experiments to study the optimal thresholding procedure involved in the use of r.

This paper is organized as follows. In Section 2, we provide the basic wavelet framework, definitions and wavelet-domain assumptions used throughout the paper. In Section 3, we state and discuss the main results on the wavelet eigenvalue regression. In Section 4, we display and discuss the broad simulation studies. In Section 5, we lay out conclusions and discuss several open problems that this work leads to. All proofs can be found in the appendix.

Framework and notation

In this section, we describe the notation, basic wavelet framework and wavelet-domain assumptions used throughout the paper. Note that this is the same wavelet framework put forth in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], Section 2.

For m ∈ N, let M (m, R) be the space of m × m real-valued matrices, and M (m 1 , m 2 , R) be the spaces of m 1 × m 2 real-valued matrices. Let S(m, R), S ≥0 (m, R) and S >0 (m, R), respectively, be the space of m × m symmetric matrices and the cones of m × m positive definite and positive semidefinite symmetric matrices. Throughout the manuscript, M denotes the operator norm of a matrix M ∈ M (p, R) in arbitrary dimension p, i.e., M = sup

u∈S p-1 u * M M * u = sup u∈S p-1 u * M * M u. For any M ∈ S(p, R), -∞ < λ 1 (M ) ≤ . . . ≤ λ q (M ) ≤ . . . ≤ λ p (M ) < ∞ (2.1)
denotes the set of ordered eigenvalues of the matrix M . For S = (s The high-dimensional asymptotic behavior of i . The plots above show the asymptotic behavior of local estimates q around each j as p and n j = n/2 j increase with their ratio fixed at p/n j = c, with c = 1/2 (left column) and c = 1/4 (right column). There are r = 6 scaling eigenvalues h i ∈ (0.1, 0.3, 0.5, 0.6, 0.8, 0.9). For each j, the reported values i are based on the octave range (j 1 , j 2 ) = (j -1, j + 1). Note the i associated with non-scaling eigenvalues stays below 0 as j increases in all instances, and the q associated with the r = 6 scaling eigenvalues approach their theoretical values. With c fixed (along each column), the larger sample size n (bottom row) results in near-convergence at smaller j (note the axes have been shifted to align the curves). Likewise, as the c decreases (between left and right columns), the near-convergence also occurs at smaller octaves j. Compared to the wavelet log-eigenvalues themselves (not shown), the multiscale estimator i displays reduced bias.

i 1 ,i 2 ) i 1 ,i 2 =1,...,n ∈ M (n,
to describe sequences of random vectors or matrices whose operator norms vanish and are bounded above, respectively, in probability.

Wavelet analysis

Throughout the paper, we make use of a wavelet multiresolution analysis (MRA; see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], chapter 7), which decomposes L 2 (R) into a sequence of approximation (low-frequency) and detail (high-frequency) subspaces V j and W j , respectively, associated with different scales of analysis 2 j . In almost all mathematical statements, we make assumptions (W 1 -W 4) on the underlying wavelet MRA. Such assumptions are accurately described in Section A. In particular, we make use of a compactly supported wavelet basis.

So, let φ and ψ be the scaling and wavelet functions, respectively, associated with the wavelet MRA. We further suppose the wavelet coefficients stem from Mallat's pyramidal algorithm [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], chapter 7). For expositional simplicity, in our description of the algorithm we use the R pvalued process Y , though analogous developments also hold for both X and Z. Initially, suppose an infinite time series {Y (k)} k∈Z , (2.4) associated with the starting scale 2 j = 1 (or octave j = 0), is available. Then, we can apply Mallat's algorithm to extract the so-named approximation (A(2 j+1 , •)) and detail (D(2 j+1 , •)) coefficients at coarser scales 2 j+1 by means of an iterative procedure. In fact, as commonly done in the wavelet literature, we initialize the algorithm with the process

R p Y (t) := k∈Z Y (k)φ(t -k), t ∈ R. (2.5)
By the orthogonality of the shifted scaling functions {φ(• -k)} k∈Z , [START_REF] Stoev | Estimation of the self-similarity parameter in linear fractional stable motion[END_REF], proof of Lemma 6.1, or Moulines et al. (2007b), p. 160;cf. Abry and Flandrin (1994), p. 33). In other words, the initial sequence, at octave j = 0, of approximation coefficients is given by the original time series. To obtain approximation and detail coefficients at coarser scales, we use Mallat's iterative procedure

R p A(2 0 , k) = R Y (t)φ(t -k)dt = Y (k), k ∈ Z (2.6) (see
A(2 j+1 , k) = k ∈Z u k -2k A(2 j , k ), D(2 j+1 , k) = k ∈Z v k -2k A(2 j , k ), j ∈ N ∪ {0}, k ∈ Z,
(2.7) where the (scalar) filter sequences

{u k := 2 -1/2 R φ(t/2)φ(t -k)dt} k∈Z , {v k := 2 -1/2
R ψ(t/2)φ(t -k)dt} k∈Z are called low-and high-pass MRA filters, respectively. Due to the assumed compactness of the supports of ψ and of the associated scaling function φ (see condition (A.2)), only a finite number of filter terms is nonzero, which is convenient for computational purposes [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]). Hereinafter, we assume without loss of generality that supp(φ) = supp(ψ) = [0, T ] (cf. [START_REF] Moulines | On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter[END_REF], p. 160). Moreover, the wavelet (detail) coefficients D(2 j , k) of Y can be expressed as

R p D(2 j , k) = ∈Z Y ( )h j,2 j k-, (2.8)
where the filter terms are defined by

R h j, = 2 -j/2 R φ(t + )ψ(2 -j t)dt.
(2.9)

If we replace (2.4) with the realistic assumption that only a finite length time series

{Y (k)} k=1,...,n (2.10) is available, writing Y (n) (t) := n k=1 Y (k)φ(t -k), we have Y (n) (t) = Y (t) for all t ∈ (T, n + 1) (cf. Moulines et al. (2007b)). Noting D(2 j , k) = R Y (t)2 -j/2 ψ(2 -j t -k)dt and D (n) (2 j , k) = R Y (n) (t)2 -j/2 ψ(2 -j t -k)dt, it follows that the finite-sample wavelet coefficients D (n) (2 j , k) of Y (n) (t) are equal to D(2 j , k) whenever supp ψ(2 -j • -k) = (2 j k, 2 j (k + T )) ⊆ (T, n + 1). In other words, D (n) (2 j , k) = D(2 j , k), ∀(j, k) ∈ {(j, k) : 2 -j T ≤ k ≤ 2 -j (n + 1) -T }. (2.11)
Equivalently, such subset of finite-sample wavelet coefficients is not affected by the so-named border effect (cf. [START_REF] Craigmile | Wavelet-based parameter estimation for polynomial contaminated fractionally differenced processes[END_REF], [START_REF] Percival | Wavelet Methods for Time Series Analysis[END_REF], [START_REF] Didier | Adaptive wavelet decompositions of stationary time series[END_REF]). Moreover, by (2.11) the number of such coefficients at octave j is given by n j = 2 -j (n+1-T )-T . Hence, n j ∼ 2 -j n for large n. Thus, for notational simplicity we suppose

n j = n 2 j (2.12)
holds exactly and only work with wavelet coefficients unaffected by the border effect.

Wavelet random matrices and assumptions

Throughout the paper, we assume observations stem from the model (1.2). The independent "signal" X = {X(t)} and the noise component Z = {Z(t)} are R r -valued and R p -valued, respectively, where r is fixed and p = p(n). The deterministic matrix P = P (n) can be expressed as

M (p, r, R) P (n) = p 1 (n), . . . , p r (n) , p q (n) = 1, q = 1, . . . , r. (2.13)
For the sake of clarity and mathematical generality, we state the appropriate conditions for the convergence in probability as well as for the asymptotic normality of wavelet log-eigenvalues directly in the wavelet domain. For j ∈ N, k ∈ Z and a dyadic sequence {a(n)} n∈N , the random vectors

D(a(n)2 j , k) ∈ R p , D X (a(n)2 j , k) ∈ R r and D Z (a(n)2 j , k) ∈ R p
denote the wavelet transform at scale a(n)2 j of the stochastic processes Y , X or Z, respectively. Whenever well-defined, the wavelet random matrix -or sample wavelet (co)variance -of Y at scale a(n)2 j is denoted by

S ≥0 (p, R) W(a(n)2 j )) = 1 n a,j n a,j k=1 D(a(n)2 j , k)D(a(n)2 j , k) * , n a,j = n a(n)2 j .
(2.14)

The remaining wavelet random matrix terms W X , W XZ , W Z are naturally defined as

W X (a(n)2 j )) = 1 n a,j n a,j k=1 D X (a(n)2 j , k)D X (a(n)2 j , k) * ∈ S ≥0 (r, R), W X,Z (a(n)2 j )) = 1 n a,j n a,j k=1 D X (a(n)2 j , k)D Z (a(n)2 j , k) * ∈ M (r, p, R), W Z (a(n)2 j )) = 1 n a,j n a,j k=1 D Z (a(n)2 j , k)D Z (a(n)2 j , k) * ∈ S ≥0 (p, R).
(2.15)

Further define the auxiliary random matrix

B a (2 j ) = P -1 H a(n) -H-(1/2)I W X (a(n)2 j )a(n) -H * -(1/2)I (P * H ) -1 ∈ S ≥0 (r, R), (2.16)
as well as its mean B a (2 j ) := E B a (2 j ). The matrix B a (2 j ) should be interpreted as a version of W X (a(n)2 j ) after compensating for scaling and non-canonical coordinates (cf. relation (2.26), which displays canonical scaling). In (2.16), we assume that the scaling matrix H has the Jordan form

H = P H diag(h 1 , . . . , h r )P -1 H , P H ∈ GL(r, R), -1/2 < h 1 ≤ . . . ≤ h r < ∞.
(2.17)

We make use of the following assumptions in the main results of this paper (Section 3). For expository purposes, we first state the assumptions, and then provide some interpretation.

Assumption (A1): The wavelet random matrix

W(a(n)2 j ) = P (n)W X (a(n)2 j )P * (n) + W Z (a(n)2 j ) +P (n)W X,Z (a(n)2 j ) + W * X,Z (a(n)2 j )P * (n) (2.18)
is well defined a.s.

Assumption (A2): In (2.18), for H as in (2.16), max EW Z (a(n)2 j ) , W Z (a(n)2 j ) , a(n) -H-(1/2)I W X,Z (a(n)2 j )) = O P (1). (2.19) Assumption (A3): The random matrix B a (2 j ) as in (2.16) satisfies √ n a,j (vec S B a (2 j ) -vec S B a (2 j )) j=j 1 ,...,jm d → N (0, Σ B (j 1 , . . . , j m )), n → ∞.
(2.20)

for some Σ B (j 1 , . . . , j m ) ∈ S ≥0 m • r(r+1) 2 , R . In addition, for some β > 0,

B a (2 j ) -B(2 j ) = O(a(n) -β ), n → ∞, j = j 1 , . . . , j m , (2.21)
where B(2 j ) ∈ S >0 (r, R).

(2.22)

Assumption (A4): The dimension p(n) and the scaling factor a(n) satisfy the relations

a(n) ≤ n 2 j 2 , a(n) n + n a(n) 2 +1 → 0, p(n) n/a(n) → c ∈ [0, ∞), n → ∞, (2.23)
where > 0 is defined as

= min min

{2≤i≤r :

h i -h i-1 >0} (h i -h i-1 ), β, h 1 2 + 1 4 . (2.24)
Assumption (A5): Let P (n) ∈ M (p, r, R) and P H ∈ GL(r, R) be as in (2.13) and (2.17), respectively. Let P (n

)P H = Q(n)R(n) be the QR decomposition of P (n)P H , where R(n) ∈ GL(r, R) and Q(n) ∈ M (p, r, R
) has orthonormal columns. Let be as in (2.24). Then, there exists a (deterministic

) matrix A ∈ S >0 (r, R) with Cholesky decomposition A = R * R such that R(n) -R = O(a(n) -).
(2.25)

Assumptions (A1 -A3) are stated in the wavelet domain. Assumption (A1) holds under very general conditions. In fact, since X and Z are assumed independent, it suffices that W X (a(n)2 j ) < ∞ and W Z (a(n)2 j ) < ∞ a.s., which holds as long as X and Z are welldefined discrete time stochastic processes. Assumption (A2) ensures that the influence of the random matrices W X,Z (a(n)2 j ) and W Z (a(n)2 j ) in the observed wavelet spectrum W(a(n)2 j ) is not too large. Assumption (A3) posits the asymptotic normality of the (wavelet domain) fractional component W X (a(n)2 j ) after compensating for scaling and non-canonical coordinates.

In turn, Assumption (A4) controls the divergence rates among n, a(n) and p(n). In particular, it states that the scaling factor a(n) must blow up slower than n, and that the threecomponent ratio p(n)a(n) n must converge to a constant (cf. the traditional ratio lim n→∞ p(n)/n for high-dimensional sample covariance matrices). Assumption (A5) ensures that, asymptotically speaking, the angles between the column vectors of the matrix P (n)P H converge in such a way that the matrix of asymptotic angles lim n→∞ P * H P * (n)P (n)P H = A has full rank. This entails that P (n) does not strongly perturb the scaling properties of the hidden random matrix W X (a(n)2 j ).

In the following examples, to fix ideas we briefly illustrate contexts where assumptions (A2) and (A3) hold.

Example 2.1 Suppose that, for any fixed n, X is an r-variate ofBm and the p(n)-variate process Z are independent. Further assume that, for fixed integers p and q, Z is made up of entry-wise independent, i.d. ARMA(p, q) processes. Then, under mild regularity assumptions it can be shown that condition (2.19) (i.e., (A2)) is satisfied. It can further be proved that relations (2.20) and (2.21) hold for the matrix sequences B a (2 j ), B a (2 j ) ∈ S ≥0 (R, r) and a matrix B(2 j ) ∈ S >0 (R, r) (i.e., (A3)) (see Example 2.1 and Proposition 4.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]; cf. Lemma C.2 in [START_REF] Abry | Wavelet estimation for operator fractional Brownian motion[END_REF]). Moreover, when H is diagonalizable with real eigenvalues, it can be easily shown that the matrix B(2 j ) satisfies the entrywise (i.e., along canonical axes) scaling relations

B(2 j ) = b(2 j ) ii i,i =1,...,r = 2 j(h i +h i +1) b(1) ii i,i =1,...,r . (2.26) Example 2.2 Let X = {X(t)} t∈Z = (X 1 (t), . . . , X r (t)) *
t∈Z be a possibly non-Gaussian, R r -valued stochastic process whose entry-wise components are independent linear fractional processes with finite fourth moments. Then, based on the framework constructed in [START_REF] Roueff | Asymptotic normality of wavelet estimators of the memory parameter for linear processes[END_REF], one can show that the associated random matrix B a (2 j ) ∈ S >0 (r, R) satisfies conditions (2.20) and (2.21) (i.e., (A3)) under conditions (W 1 -W 4) and mild additional assumptions on the wavelet ψ and on the process X (see Proposition C.2 in Abry et al. ( 2022)). For this instance of X, the model (1.2) is associated with the so-named blind source separation problems in the field of signal processing (e.g., [START_REF] Comon | Handbook of Blind Source Separation: Independent Component Analysis and Applications[END_REF]; see also [START_REF] Abry | Two-step wavelet-based estimation for Gaussian mixed fractional processes[END_REF] on fractional instances).

Example 2.3 Recall that a distribution is called sub-Gaussian when its tails are no heavier than those of the Gaussian distribution [START_REF] Vershynin | High-Dimensional Probability: an Introduction with Applications in Data Science[END_REF], Proposition 2.5.2). Sub-Gaussian distributions form a broad family that includes the Gaussian distribution itself, as well as compactly supported distributions, for example. Suppose the noise process {Z(t)} t∈Z consists of i.i. 2022), Example 4.3). On possible extensions beyond the sub-Gaussian case as well as on related results, see [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF] and [START_REF] Einmahl | Characterization of LIL behavior in Banach space[END_REF].

Main results

As mentioned in the Introduction, under the assumptions laid out in Section 2.2 the r largest eigenvalues of the random matrix W(a(n)2 j ) display asymptotic scaling relationships in high dimensions dictated by the eigenvalues of H (see (1.6)). Furthermore, under additional assumptions, the fluctuations of the rescaled wavelet log-eigenvalues around their limits are asymptotically Gaussian. For the reader's convenience, the precise statements are provided in Theorems B.1 and B.2. These facts indicate that the scaling eigenvalues (1.4) can be efficiently estimated in high dimensions through a linear regression based on the log-eigenvalues of W(a(n)2 j ) over a range of scales 2 j 1 , . . . , 2 j 2 , provided the scaling factor a(n) is sufficiently large. This motivates the following definition.

Definition 3.1 Let {W(a(n)2 j )} j=j 1 ,...,j 2 be the wavelet random matrices corresponding to scales {a(n)2 j 1 , . . . , a(n)2 j 2 }. Fix a range of octaves

j = j 1 , j 1 + 1, . . . , j 2 , m := j 2 -j 1 + 1. (3.1)
The wavelet eigenvalue regression is given by

{ i } i=1,...,p(n) := 1 2 j 2 j=j 1 w j log 2 λ i W(a(n)2 j ) -1 i=1,...,p(n) . (3.2)
In (3.2), w j , j = j 1 , . . . , j 2 , are weights satisfying the relations

j 2 j=j 1 w j = 0, j 2 j=j 1 jw j = 1, (3.3)
where

w j = 1 if m = 1 (see (3.1)).
If W(a(n)2 j ) ∈ S >0 (r, R), then the estimator (3.2) is well defined. Assuming for the moment that the dimension r of the hidden fractional process X is known, expression (3.2) based on the top r wavelet eigenvalues can be interpreted as wavelet eigenstructure estimators of the scaling eigenvalues, i.e., we can write

{ h q } q=1,...,r = { p(n)-r+q } q=1,...,r . (3.4)
Hence, the estimator (3.4) can be viewed as a high-dimensional extension of the wavelet eigenvalue regression estimator first developed in Abry andDidier (2018b, 2018a) in a noiseless, lowdimensional context for ofBm.

In the following theorems, we characterize the asymptotic and large scale behavior of the wavelet eigenvalue regression estimator (3.4) in high dimensions, namely, its consistency and, under assumptions, asymptotic normality. Consistency is, indeed, a direct consequence of the convergence in probability of the r largest rescaled wavelet eigenvalues (cf. Theorem B.1). By contrast, establishing the asymptotic normality of the wavelet eigenvalue regression estimator goes beyond the asymptotic distribution of wavelet eigenvalues over large scales in high dimensions (cf. Theorem B.2). In fact, it also requires establishing the rate of convergence of the rescaled r largest eigenvalues of the deterministic wavelet matrices EW(a(n)2 j ). In turn, this involves studying the rate of angular convergence of the associated (deterministic) wavelet eigenvectors in terms of the high-dimensional coordinates P = P (n).

Theorem 3.1 Assume (W 1 -W 4) and (A1 -A5) hold and fix j 1 , j 2 ∈ N, j 1 < j 2 . Then, h q -h q q=1,...,r P → 0, n → ∞.

(3.5)

Theorem 3.2 Suppose (W 1 -W 4) and (A1 -A5) hold, and fix j 1 , j 2 ∈ N, j 1 < j 2 . Further suppose one of the following conditions holds, namely,

(i) either -1/2 < h 1 < . . . < h r < ∞; (3.6)
(ii) or h 1 = . . . = h r and the functions ξ q (2 j ) in (B.1) satisfy

q 1 = q 2 ⇒ ξ q 1 (1) = ξ q 2 (1). (3.7) Then, n a(n) h q -h q q=1,...,r d → N (0, M Σ λ M * ), (3.8)
as n → ∞, for some weight matrix M (see (C.37)) and Σ λ as in Theorem B.2.

Remark 3.1 Under non-simple scaling matrix eigenvalues and without condition (3.7), the asymptotic distribution of { h q } q=1,...,r is generally expected to be non-Gaussian because, in this case, wavelet log-eigenvalues themselves are generally expected to be non-Gaussian (see the discussion in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], Remark 3.1). The broad characterization of the distribution of wavelet log-eigenvalues outside the framework of the assumptions of Theorem 3.2 is currently a topic of research.

In practice, the dimension r of the latent process X may not be known. In the following definition, we introduce an estimator of r based on the eigenvalues of W(a(n)2 j ).

Definition 3.2 Let j 1 < j 2 and consider any set of weights

v j 1 , . . . v j 2 ∈ R satisfying j 2 j=j 1 v j = 1. Let ∆ i (j 1 , j 2 ) := j 2 j=j 1 v j log λ i W(a(n)2 j ) log(a(n)2 j ) (3.9)
Given κ > 0, we define r(2 j 1 , 2 j 2 , κ) := #{i : ∆ i (j 1 , j 2 ) > κ}.

(3.10) Now note that, under the assumptions of Theorem 3.1, the lowest p -r eigenvalues stay bounded. Hence, for i ≤ p -r, the quantity ∆ i (j 1 , j 2 ) tends to zero in probability. On the other hand, for i > p -r, still under the assumptions of Theorem 3.1, ∆ i (j 1 , j 2 ) converges to 2h i-(p-r) + 1 > 0 in probability, thereby separating the non-scaling and scaling eigenvalues of W(a(n)2 j ). This phenomenon lies behind the following theorem, which establishes the consistency of the estimator r(2 j 1 , 2 j 2 , κ).

Theorem 3.3 Let r(2 j 1 , 2 j 2 , κ) be as in (3.10) with κ ∈ (0, 2h 1 + 1] and suppose the assumptions of Theorem 3.1 hold. Then,

lim n→∞ P r(2 j 1 , 2 j 2 , κ) = r = 1 (3.11)
Remark 3.2 In the context of the (univariate) wavelet regression, it is common practice to select the weights (3.3) to include information about the wavelet variance over the regression scales 2 j 1 , . . . 2 j 2 (cf. [START_REF] Veitch | A wavelet-based joint estimator of the parameters of long-range dependence[END_REF]). This can lead to improved estimation performance. For convenience, in all simulations we set v j = jw j , where w j are the regression weights in (3.2).

Monte Carlo studies

In this section, we describe our Monte Carlo-based studies of the asymptotic behavior of the wavelet eigenvalue regression estimator { i } i=1,...,p(n) (see (3.2)) and of the estimator r(j 1 , j 2 , κ) (see (3.10)) in high dimensions.

Figure 1 displays the asymptotic behavior of { i } i=1,...,p(n) over scales [j-1, j+1] in the specified ranges of j. As the effective sample size n j = n/2 j and dimension p = p(n) increase together at the specified ratio c := p/n j , the estimators display convergence to their theoretical values or boundedness. As expected, at more favorable (namely, smaller) values of c, the eigenvalues display near-convergence to theoretical values at smaller sample sizes. Computational studies not shown reveal that the asymptotic performance of i over multiple scales in Figure 1 displays reduced bias by comparison to using a single or a reduced number of scales.

Recall that the so-named Gamma plots are defined as plots of the empirical quantiles of squared Mahalanobis distance vs. the theoretical quantiles of a χ 2 r distribution [START_REF] Johnson | Applied Multivariate Statistical Analysis[END_REF], [START_REF] Oppong | Assessing univariate and multivariate normality: A guide for non-statisticians[END_REF]). Figure 2 displays Gamma plots of the sample distribution of the wavelet eigenvalue regression estimator { h q } q=1,...,6 = { p(n)-r+q } q=1,...,6 (r = 6) for various combinations of p and n j 2 = n/2 j 2 . In all instances, the hidden process X was simulated by means of CME as an ofBm with scaling (Hurst) parameter H = diag 0.1, 0.3, 0.5, 0.6, 0.8, 0.9 and (instantaneous) covariance matrix EB H (1)B H (1) = Toeplitz(1, 0.2, 0.2, 0.3, 0.2, 0.3). At each time t, the vector noise term Z(t) consisted of p i.i.d. N (0, 1) entries. At each run, the columns of the coordinates matrix P (n) were set to the first r = 6 canonical vectors. As n (and, hence, p) grows, in the high-dimensional limit the distribution of { h q } q=1,...,6 becomes increasingly closer to a multivariate Gaussian, as mathematically characterized in Theorem 3.2.

Figure 3 depicts the performance r(j 1 , j 2 , κ) as a function of the threshold κ > 0 for various combinations of p/n j 2 , n j 2 and r. To illustrate the asymptotic behavior of r(j 1 , j 2 , κ) while keeping the ratio c := p/n j 2 constant, the parameters j 1 , j 2 were chosen as follows: for n = 2 10 , (j 1 , j 2 ) = (3, 5) (hence, n j 2 = 32); for n = 2 12 , (j 1 , j 2 ) = (4, 6) (hence, n j 2 = 64); for n = 2 14 , (j 1 , j 2 ) = (5, 7) (hence, n j 2 = 128). The hidden process X was simulated by means of CME as an ofBm with scaling (Hurst) parameter H = diag h 1 , h 2 , ..., h r , where h q = q/(r + 1), q = 1, . . . , r. At each time t, the vector noise term Z(t) consisted of p i.i.d. N (0, 1) entries. At each run, the entries of the coordinates matrix P (n) were generated as i.i.d. standard normal entries, and then the norm of each column was renormalized to 1 so as to keep constant the signal-to-noise ratio. For each wavelet random matrix, eigenvalues that fell below 10 -10 in absolute value (stemming from deficient rank when p/n j < 1) were set to 0 manually to avoid spuriously high values of ∆ i (j 1 , j 2 ).

The computational studies reveal that, in all instances, too small or too large a value of κ ∈ (0, 1) leads to over-and underestimation, respectively, of X components. For any r ∈ {3, 5, 8}, and for each value of c, larger n j 2 (and, hence, larger dimension p) produces a larger range of κ where r(j 1 , j 2 , κ) is perfectly accurate or nearly so, thus reflecting the high-dimensional asymptotics in Theorem 3.3. Also, for each pair n j 2 and p, larger values of c lead to narrower ranges of κ for which r(j 1 , j 2 , κ) concentrates around the true value of r. This is expected since larger values of c amount to more extreme stochastic regimes.

Conclusion and open problems

In this paper, we construct the wavelet eigenvalue regression methodology (Abry andDidier (2018b, 2018a)) in high dimensions. We assume that possibly non-Gaussian, finite-variance pvariate measurements are made of a low-dimensional r-variate (r p) fractional stochastic process with non-canonical scaling coordinates and in the presence of additive high-dimensional noise. The measurements are correlated both time-wise and between rows. Due to the asymptotic and large scale properties of large wavelet random matrices, the wavelet eigenvalue regression is shown to be consistent and (under additional assumptions) asymptotically Gaussian in the estimation of the fractal structure of the underlying measurements. We further construct a consistent estimator of the effective dimension r of the system that significantly increases the robustness of the statistical methodology. The estimation performance over finite samples is further studied by means of simulations.

This research leads to many new research directions, some of which can be briefly described as follows. (i) In applications, the results in this paper naturally pave the way for the investigation of scaling behavior in high-dimensional data from fields such as physics, neuroscience and signal processing; (ii) Modeling requires a deeper study, in the wavelet domain, of the so-named eigenvalue repulsion effect (e.g., [START_REF] Tao | Topics in Random Matrix Theory[END_REF]), which may severely skew the observed scaling laws when the assumptions of Theorem 3.2 are violated. This is particularly important in the context of instances where all scaling eigenvalues are close to equal, with the same holding for the asymptotic rescaled eigenvalues ξ q (2 j ) (see [START_REF] Wendt | Bootstrap-based bias reduction for the estimation of the self-similarity exponents of multivariate time series[END_REF] on preliminary computational studies). In those cases, it is of great interest to develop efficient testing procedures for the statistical identification of distinct scaling eigenvalues in real data; (iii) An interesting direction of extension is the construction of statistical methodology for instances where the r largest eigenvalues of wavelet random matrices exhibit non-Gaussian fluctuations (cf. Remark 3.1), or other related instances where the hidden process X displays heavier tails. The asymptotic (joint) distribution of { h q } q=1,...,r = { p-r+q } q=1,...,r in high dimensions. Each plot displays Gamma plots (see text in Section 4 for a definition) based on 5000 independent realizations of { h q } q=1,...,6 for increasing sample size (from top to bottom, respectively) and p/n j2 = 1/2, 1, 2 (left to right column, respectively). The plots further contain Kolmogorov-Smirnov test decisions d KS for the null hypothesis that the Mahalanobis distance follows a χ 2 6 distribution (obtained as averages over 100 random subsets of size 1250 of the 5000 realizations). The plots show that, for any ratio p/n j2 considered, as the sample size increases the distribution of { h q } q=1,...,6 becomes statistically indistinguishable from a joint Gaussian distribution. p/n j 2 = 1/2: p/n j 2 = 1: p/n j 2 = 2: Figure 3: The asymptotic performance of r(j 1 , j 2 , κ). Each plot displays r versus κ based on 5000 observations of {∆ 1 (j 1 , j 2 ), . . . , ∆ p (j 1 , j 2 )} for p/n j2 ∈ {0.5, 1, 2} (top to bottom, respectively), over sample sizes n ∈ {2 10 , 2 12 , 2 14 } (left to right column, respectively). The solid lines represent the Monte Carlo average for the particular choice of r, n and p. The upper/lower limits of shaded regions correspond the empirical regions between quantiles 0.05 and 0.95. The marked lines (with small "x" markers) are the points at which the mean was exactly equal to the corresponding value of r.

exist, and the deterministic functions ξ q satisfy the scaling relation ξ q (2 j ) = 2 j (2hq+1) ξ q (1).

(B.2)

In addition,

λ 1 W(a(n)2 j ) ≤ . . . ≤ λ p-r W(a(n)2 j ) = O P (1). (B.3)
Moreover, the asymptotic normality of the r largest wavelet log-eigenvalues in high dimensions was also established. The statement covers the central cases of simple (h 1 < . . . < h r ) and of identical (h 1 = . . . = h r ) scaling eigenvalues. Note that, by comparison to Theorem B.1, the asymptotic normality of wavelet log-eigenvalues requires some additional condition so as to ensure the finite-sample differentiability of the eigenvalues with respect to scaling parameters (see also Remark 3.1).

Theorem B.2 (a consequence of Theorem 3.2 in Abry et al. ( 2022)) Under the assumptions of Theorem 3.2, for p = p(n),

R m•r √ n a,j log λ p-r+q W(a(n)2 j ) -log λ p-r+q EW(a(n)2 j ) q=1,...,r j=j 1 ,...,j 2 d → N (0, Σ λ ) (B.4) as n → ∞ for some Σ λ ∈ S ≥0 (m • r, R).
Remark B.1 Note that Assumption (A4) (see (2.23) and (2.24)) implies the analogous Assumption (A4) of [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], which in turn is assumed in Theorems 3.1 and 3.2 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. The slightly stronger conditions (2.23) and (2.24), as well as (2.25), are needed in Theorem 3.2 in this paper, i.e., in the context of the wavelet eigenvalue regression.

C Proofs

In this section, we provide the proofs the main results of the paper. We often use the following notation. For an arbitrary q ∈ {1, . . . , r}, I -, I 0 and I + are index sets given by the relations

I -:= { : h < h q }, I 0 := { : h = h q }, I + := { : h > h q }. (C.1)
Note that I -and I + are possibly empty. Also write

r 1 := card(I -), r 2 := card(I 0 ) ≥ 1, r 3 := card(I + ). (C.2)
Throughout this section, for notational simplicity, we write

P (n)P H ≡ P (n) ≡ P, (C.3)
whose column vectors are denoted by p (n) ≡ p , = 1, . . . , r. Also for notational convenience, we may write a ≡ a(n).

Proof of Theorem 3.1: By Theorem B.1, by the properties (3.3) of the regression weights w j , j = 1, . . . , m, and by relation (B.2),

h q = - 1 2 + 1 2 j 2 j=j 1 w j log 2 λ p-r+q (W(a(n)2 j )) = - 1 2 + 1 2 j 2 j=j 1 w j log 2 λ p-r+q (W(a(n)2 j ))/a(n) 2hq+1 P → - 1 2 + 1 2 j 2 j=j 1 w j log 2 ξ q (2 j ) = - 1 2 + 1 2 j 2 j=j 1 w j log 2 2 j (2hq+1) ξ q (1) = h q .
The proof of Theorem 3.2 mainly relies on Proposition C.1, stated and proved next. The proposition establishes not only that the r largest rescaled eigenvalues of the deterministic matrix EW(a2 j ) converge to their respective asymptotic rescaled eigenvalues ξ q (2 j ) (cf. expression (B.1)), but also it provides an upper bound on the associated convergence rate.

Proposition C.1 Fix j ∈ N and suppose conditions (W 1 -W 4) and (A1 -A5) hold. Suppose that either (i) 0 < h 1 < . . . < h r < 1; or (ii) h 1 = . . . = h r , and whenever q 1 = q 2 , the functions ξ q (2 j ) in (B.1) satisfy

ξ q 1 (1) = ξ q 2 (1). (C.4)
Then, for some C > 0 that does not depend on j,

λ p-r+q (EW(a(n)2 j )) a(n) 2hq+1 -ξ q (2 j ) ≤ C a(n)
, q = 1, . . . , r.

(C.5) Proof: Assume condition (i) holds. We prove only the statement for q ∈ {2, . . . , r -1} since the cases q = 1, r can be handled by a simplified version of the same argument. For

γ = γ 1 , . . . , γ r * ∈ R r , x ∈ R r 3 ,
consider the deterministic function

g n (γ; x) := y * n (γ; x)B a (2 j )y n (γ; x) ∈ R, (C.6) where y n (γ; x) = γ 1 a(n) h 1 -hq , . . . , γ r 1 a(n) hr 1 -hq r 1 , γ q , x q+1 , . . . , x r r 3 * . Now let g(γ; x) = y * (γ; x)B(2 j )y(γ; x) ∈ R, (C.7) where R r y(γ; x) = 0, . . . , 0 r 1 , γ q , x q+1 , . . . , x r r 3 * = lim n→∞ y n (γ; x).
(observe that y(γ; x) depends on γ only through its q-th entry γ q ). For = 1, . . . , p, let u (n) denote a unit eigenvector of EW(a2 j ) associated with its -th eigenvalue in nondecreasing order. Define γ q (n) := P * (n)u p-r+q (n).

(C.8) By Proposition D.1, u p-r+q (n) can be chosen so that the limit lim n→∞ γ q (n) =: γ q = (γ 1,q , . . . , γ r,q ) * , n → ∞, (C.9) exists. So, let x * (n) ∈ R r 3 and x * = (x * ,q+1 , . . . , x * ,r ) ∈ R r 3 (C.10) be the (unique) minimizers of the functions g n (γ q (n); •) and g(γ q ; •), respectively, where such functions are given in (C.6) and (C.7). Observe that, as n → ∞, x * (n) → x * , implying y n (γ q (n), x * (n)) → y(γ q , x * ). Hence,

g n (γ q (n), x * (n)) → g(γ q , x * ), n → ∞. (C.11)
By Lemma D.3, for large enough n we may take a sequence of unit vectors v(n) ∈ span{u p-r+q (n), . . . , u p (n)} such that

p (n), v(n) = x * , a h -hq , = q + 1, . . . , r, P * (n)v(n) → γ q . (C.12)
Therefore, as n → ∞,

a h-hqI P * (n)v(n) = y n P * (n)v(n), x * (n) → y(γ q ; x * ), which implies g n (P * (n)v(n), x * (n)) → g(γ q , x * ), n → ∞. (C.13) Moreover, let R r 3 x(n) := p q+1 (n), u p-r+q (n) a h q+1 -hq , . . . , p r (n), u p-r+q (n) a hr-hq , n ∈ N
(not to be confused with the minimizer x * (n) of g n (γ q (n), •) as in (C.10)). Then, we can express

λ p-r+q EW(a2 j ) a 2hq+1 = g n (γ q (n), x(n)) + u * p-r+q (n) EW Z (a2 j ) a 2hq+1 u p-r+q (n).
Thus, for all large n and for γ q (n) and v(n) as in (C.8) and (C.12), respectively,

g n (γ q (n), x * (n)) + u * p-r+q (n) EW Z (a2 j ) a 2hq+1 u p-r+q (n) ≤ g n (γ q (n), x(n)) + u * p-r+q (n) EW Z (a2 j ) a 2hq+1 u p-r+q (n) = u * p-r+q (n) EW(a2 j ) a 2hq+1 u p-r+q (n) = λ p-r+q (EW(a2 j )) a 2hq+1 ≤ v * (n) EW(a2 j ) a 2hq+1 v(n) = g n (v(n); x * ) + v * (n) EW Z (a2 j ) a 2hq+1 v(n), (C.14)
where in the second inequality we used the fact that v(n) ∈ span{u p-r+q (n), . . . , u p (n)}. However, note that EW Z (a2 j )/a 2hq+1 = o(a -) under assumption (2.19). Then, by Corollary D.1 with

M n = EW Z (a2 j )/a 2hq+1 , λ p-r+q (EW(a(n)2 j )) a(n) 2hq+1 → ξ q (2 j ), n → ∞. (C.15)
Therefore, in view of (C.11) and (C.13), by taking limits in (C.14), we see that

ξ q (2 j ) = g(γ q , x * ). (C.
λ p-r+q (EW(a(n)2 j )) a(n) 2hq+1 -ξ q (2 j ) ≤ R(n) -R + R Q * (n)u p-r+q (n) -e q = O(a -),
where in the last equality we make use of condition (2.25). Hence, γ q = lim n→∞ γ q (n) = Re q . Keeping in mind expressions (C.20) and (C.21), this implies that

|y q (n) -y q | = | p q (n), u p-r+q (n) -γ qq | = O(a -). (C.26)
In other words, (C.23) also holds for = q.

Turning to y (n) -y for ∈ I + = {q + 1, . . . , r}, recall that

x * (n) := (x * ,q+1 (n), . . . , x * ,r (n)) *
is the unique minimizer of g n (γ q (n); •) in (C.6), and that x * = lim n→∞ x * (n). Now, reexpress

B a (2 j ) = (B i (n)) i, =1,2,3 , B(2 j ) = (B i ) i, =1,2,3 , (C.27)
where B i (n) and B i denote blocks of size r i × r . The first order conditions for the minimization of g n (γ q (n); •) imply that

B 31 (n) B 32 (n) B 33 (n)   O(a -) p q (n), u p-r+q (n) x * (n)   = 0 (C.28)
where O(a -) ∈ R r 1 contains the first r 1 terms in y(n) (see (C.20)). Since B 33 (n) ∈ S >0 (r 3 , R) for all large n, by rearranging (C.28) and multiplying by B 33 (n) -1 on the left we have, for all large n,

x * (n) = -B 33 (n) -1 O(a -) + B 32 (n) p q (n), u p-r+q (n) , x * = -B -1 33 B 32 γ qq .
Together with (C.24) and (C.26), as well as assumption (2.21), this implies that (y (n) -y

) ∈I + = x * (n) -x * ≤ O(a -) + B 33 (n) -1 B 32 (n) p q (n), u p-r+q (n) -B -1 33 B 32 γ qq ≤ O(a -)+ (B 33 (n) -1 B 32 (n)-B -1 33 B 32 p q (n), u p-r+q (n) + B -1 33 B 32 p q (n), u p-r+q (n) -γ qq ≤ O(a -) + C (B 33 (n) -1 B 32 (n) -B -1 33 B 32 + p q (n), u p-r+q (n) -γ qq = O(a -).
Hence, (C.23) also holds for ∈ I + . Thus, (C.23) holds for = 1, . . . , r. This establishes (C.18) and, hence, (a). Now, we turn to (b). Consider the sequence of unit vectors v(n) given by (D.11) in Lemma D.3 for the deterministic matrix EW(a2 j ). For a sequence of scalars {c q

(n), c q+1 (n), . . . , c r (n)} satisfying r i=q c 2 i (n) = 1, we may write v(n) = r i=q c i (n)u p-r+i (n). By expression (D.11) of Lemma D.3, 1 -c q (n) 2 = r i=q+1 c i (n) 2 = O(a -2 ).
Moreover, by flipping the sign of u p-r+q (n) if necessary, we may suppose c q (n) → 1 as n → ∞. Hence, |1 -c q (n)| = O(a -2 ). Therefore,

| p q (n), v(n) -p q (n), u p-r+q (n) | = (c q (n) -1) p q (n), u p-r+q (n) + i∈I + c i (n) p q (n), u p-r+i (n) ≤ |1 -c q (n)| • | p q (n), u p-r+q (n) | + O(a -) = O(a -).
(C.29) Thus, (C.26) and (C.29) show that

p q (n), v(n) -γ qq ≤ p q (n), v(n) -p q (n), u p-r+q (n) + p q (n), u p-r+q (n) -γ qq = O(a -) + O(a -) = O(a -). (C.30) Now, define y(n) = a h-hqI P * (n)v(n) = y n P * (n)v(n), x * .
Thus, entry-wise, we can express

y(n) =      p (n), v(n) a h -hq = O(a -), ∈ I -; p (n), v(n) , = q; x * , , ∈ I + .
From (C.30) and for y as in (C.21), we obtain

y(n) -y 2 = p q (n), v(n) -γ qq 2 + O(a -2 ) = O(a -2 ). (C.31) Hence, g n (v(n); x * ) -ξ q (2 j ) = y * (n)B a (2 j ) y(n) -y * B(2 j )y = ( y(n) -y) * B a (2 j )( y(n) -y) + 2y * B(2 j )( y(n) -y) + y * B a (2 j ) -B(2 j ) y ≤ C B a (2 j ) -B(2 j ) + C y(n) -y = O(a -),
where the last equality is a consequence of condition (2.21) and of relation (C.31). This establishes (C.19) and, hence, (b). Thus, as anticipated, (C.5) holds under condition (i).

We now turn to (C.5) for the case (ii). Let W P X (a2 j ) = P W X (a2 j )P * .

Fix any q ∈ {1, . . . , r}. By Weyl's inequality (e.g., [START_REF] Vershynin | High-Dimensional Probability: an Introduction with Applications in Data Science[END_REF], Theorem 4.5.3),

λ p-r+q EW(a2 j ) a 2hq+1 - λ p-r+q EW P X (a2 j ) a 2hq+1
≤ EW(a2 j )/a 2hq+1 -EW P X (a2 j )/a 2hq+1 = O(a -).

So, it suffices to show (C.5) for the matrix EW P X (a2 j )/a 2hq+1 = P B a (2 j )P * in place of EW(a2 j ). Let u i (n), i = 1, . . . , p, be an orthonormal basis of eigenvectors of the (deterministic) matrix P B a (2 j )P * , where the associated ordering of eigenvalues is arbitrary in the case of ties. Observe that P * P B a (2 j )P * u i (n) = λ i (P B a (2 j )P * ) • P * u i (n) for i = 1, . . . , p. In particular, λ p-r+q (P B a (2 j )P * ) = λ q (P * P B a (2 j )), q = 1, . . . , r (C.32) (where P * u p-r+q (n) is the associated eigenvector of the non-symmetric matrix P * P B a (2 j ) ∈ M (r, R)). By a similar argument, and considering again the QR decomposition

P ≡ P (n) = Q(n)R(n), λ q R(n)B a (2 j )R * (n) = λ q R * (n)R(n)B a (2 j ) , q = 1, . . . , r. (C.33)
As a consequence of (C.32), (C.33) and of the fact that Q

* (n)Q(n) = I r , λ q P * (n)P (n)B a (2 j ) = λ q R(n)B a (2 j )R * (n) , q = 1, . . . , r.
So, by taking limits in (C.32), Lemma D.1 with M n ≡ 0 implies that, for q = 1, . . . , r, λ q (P * (n)P (n)B a (2 j )) → ξ q (2 j ) as n → ∞. Moreover, by a similar argument to the one leading to (C.33), ξ q (2 j ) = λ q (AB(2 j )) = λ q (RB(2 j )R * ) (C.34) due to assumptions (2.21) and (2.25). Therefore, again by Weyl's inequality,

λ q P * (n)P (n)B a (2 j ) -ξ q (2 j ) = λ q R * (n)R(n)B a (2 j ) -λ q (R * RB(2 j )) = λ q R(n)B a (2 j )R * (n) -λ q RB(2 j )R * ≤ R(n)B a (2 j )R * (n) -RB(2 j )R * ≤ (R(n) -R)B a (2 j )(R(n) -R) * + (R(n) -R)B a (2 j )R * + RB a (2 j )(R(n) -R) * + R(B a (2 j ) -B(2 j ))R * = O(a -), (C.35)
where the last equality is a consequence of conditions (2.21) and (2.25). This establishes (C.5) in the case (ii).

We are now in a position to prove Theorem 3.2 (cf. Corollary 2, (ii), in Abry and Didier (2018a)).

Proof of Theorem 3.2: For a fixed q ∈ {1, . . . , n}, the left-hand side of (3.8) can be recast in the form n a j 2 j=j 1 w j 2 log 2 λ p-r+q W(a2 j ) -log 2 λ p-r+q EW(a2 j ) + n a j 2 j=j 1 w j 2 log 2 λ p-r+q EW(a2 j ) -log 2 ξ q (a2 j ) + n a j 2 j=j 1 w j 2 log 2 ξ q (a2 j ) -h q . (C.36)

Note that by (B.2) in Theorem B.1, ξ q (a2 j ) = a 2hq+1 ξ q (2 j ). Therefore, by property (3.3), the third term in the sum (C.36) is zero. By Proposition C.1, the second term in the sum (C.36) is bounded by n a

j 2 j=j 1 |w j | 2 C a ≤ C n a 1+2 → 0, n → ∞,
where the limit is a consequence of condition (2.23). Therefore, we can rewrite the left-hand side of (3.8) as j 2 j=j 1 2 j/2-1 w j log 2 √ n a,j log λ p-r+q W(a2 j ) -log λ p-r+q EW(a2 j ) + o(1), and the weak limit (3.8) follows from Theorem B.2. In the limiting variance in (3.8), the weight matrix M ∈ M (r, mr, R) is given by

M = 2 j 1 /2 w j 1 log 2 I r ; 2 j 1 +1/2 w j 1 +1 log 2 I r ; . . . ; 2 j 2 /2 w j 2 log 2 I r , (C.37)
where I r ∈ M (r, R) is an identity matrix and m = j 2 -j 1 + 1.

Proof of Theorem 3.3: Theorem B.1 and the scaling relation (B.2) show that, as n → ∞, log 2 λ p-r+q W(a2 j ) log 2 (a2 j ) P → 2h q + 1 > 0, q = 1, . . . , r.

In addition, by relation (B.3) in Theorem B.1, log 2 λ p-r W(a2 j ) log 2 (a2 j ) P → 0 for every j. Thus, for q = 1, . . . , r and κ > 0 as in (3.10), ∆ p-r+q (j 1 , j 2 ) The following lemma is used in the proof of Proposition C.1.

P → j 2 j=j 1 v j (2h q + 1) = 2h q + 1 > κ. (C.38) Likewise, max i=1,...,p-r |∆ i (j 1 , j 2 )| ≤ j 2 j=j 1 |v j | • log 2 λ p-r W(a2 j ) log 2 (a2 j ) P → 0
Lemma C.1 Suppose the assumptions of Theorem 3.2 hold, as well as condition (i) in the same theorem. Let P

(n) = Q(n)R(n) be the QR decomposition of P (n), where R(n) ∈ GL(r, R) is upper-triangular, and Q(n) = ( p 1 (n), . . . , p r (n)) ∈ M (p,
(n), . . . , p r (n)} ⊥ span{p r+1 (n), . . . , p p (n)}.
Note that expression (D.3) of Lemma D.1 (with W(a2 j ) = EW(a2 j )) shows that for each fixed i ∈ {2, . . . , r},

p i (n), u p-r+ (n) 2 = O(a -2(h i -h ) ) = O(a -2 ), = 1, . . . , i -1. (C.42) Hence, from (C.41), r i=2 i-1 =1 p i (n), u p-r+ (n) 2 = O r i=2 i-1 =1 p i (n), u p-r+ (n) 2 = O(a -2 ). (C.43)
We now proceed inductively to show (C.25). Turning to the case k = 1, first note that, as a consequence of (C.43), r i=2 p i (n), u p-r+1 (n) 2 = O(a -2 ). In addition, Lemma D.2 shows that

p i=r+1 p i (n), u p-r+1 (n) 2 = O(a -(h 1 + 1 2 ) ) = O(a -2 ), implying 1 -p 1 (n), u p-r+1 (n) 2 = r i=2 p i (n), u p-r+1 (n) 2 + p i=r+1 p i (n), u p-r+1 (n) 2 = O(a -2 ) + O(a -2 ) = O(a -2 ),
i.e., (C.25) holds for k = 1. We now proceed by induction on k. It suffices to consider the range 2 ≤ k ≤ r. So, suppose we have 1

-p i (n), u p-r+i (n) 2 = O(a -2 ), i = 1, . . . , k. Therefore, by decomposing p i (n) in the eigenvector basis, 1 -p i (n), u p-r+i (n) 2 = =p-r+i p i (n), u (n) 2 = O(a -2 ), i = 1, . . . , k.
In particular,

p i (n), u p-r+(k+1) (n) 2 = O(a -2 ), i = 1, . . . , k. (C.44)
However, Lemma D.2 again shows that

p i=r+1 p i (n), u p-r+(k+1) (n) 2 = o(a -2 ). (C.45) Relations (C.44) and (C.45) imply that 1 -p k+1 (n), u p-r+(k+1) (n) 2 = k i=1 p i (n), u p-r+(k+1) (n) 2 + r i=k+2 p i (n), u p-r+(k+1) (n) 2 + p i=r+1 p i (n), u p-r+(k+1) (n) 2 = O(a -2 ) + O(a -2 ) + o(a -2 ) = O(a -2 ),
where we used that

r i=k+2 p i (n), u p-r+(k+1) (n) 2 = O r i=k+2 p i (n), u p-r+(k+1) (n) 2 = O(a -2
) by expressions (C.41) and (C.42). Thus, 1 -p k (n), u p-r+k (n) 2 = O(a -2 ), k = 1, . . . , r; i.e., (C.25) holds.

D Auxiliary lemmas

In this section, for the reader's convenience, we recap the statements of some useful lemmas established in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. For simplicity, we adjust the statements to the particular case involving deterministic matrices that is relevant for establishing the statements in Section C. In this section, we also use the notation (C.1)-(C.3), and write a instead of a(n) whenever convenient.

Lemma D.1 [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], Lemma B.3) Fix q ∈ {1, . . . , r} and consider the index sets I -, I 0 , and I + as in (C.1). Suppose I + = ∅, and let

W(a(n)2 j ) = P (n)a(n) h+ 1 2 I B a (2 j )a(n) h+ 1 2 I P * (n) + M n , (D.1)
where

M n = a(n) 2hq+1 O(1) a(n) 2hq+1 + P (n)a(n) h O(1) a(n) 2hq+1/2 + O(1) * a(n) h P * (n) a(n) 2hq+1/2 . (D.2)
In (D.2), the terms O(1) denote any sequence of matrices of appropriate dimension whose norms are bounded. Also, let u p-r+q (n) be a unit eigenvector associated with the (p -r + q)-th eigenvalue of W(a(n)2 j ). Then,

max ∈I + {| p (n), u p-r+q (n) | a(n) h -hq } = O(1), n → ∞. (D.3)
The following lemma is used in the proof of Lemma C.1. be an orthonormal basis for the nullspace of P * (n). Let be as in (2.24). Fix any q ∈ {1, . . . , r} and let u p-r+q (n) denote a unit eigenvector associated with the (p -r + q)-th eigenvalue of W(a(n)2 j ) as given in Lemma D.1. Then,

p i=r+1 p i (n), u p-r+q (n) 2 = O a(n) -(hq+ 1 2 ) . (D.5)
The following lemma is used in the proofs of Proposition C.1 and Lemma C.1. It is a slight refinement of Lemma B.6 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. The main difference lies in the second statement in (D.11), which provides the rate at which some of the coefficients of the constructed vector (namely, v(n )) shrink to zero as n → ∞. For the reader's convenience, we provide a proof of (D.11).

Lemma D.3 Fix q ∈ {1, . . . , r}, and for = 1, . . . , p, let u (n) be a unit eigenvector associated with the -th eigenvalue of the matrix W(a(n)2 j ) as defined in Lemma D.1. For each n, let 2022), Γ + is nonsingular. This implies that Γ + n has full rank for all large n, without loss of generality assumed so for all n. Therefore, for each fixed 0 ≤ s ≤ 1 and for the arbitrary scalars α q , ..., α r-r 3 as in (D.8), the system of equations In view of (D.14) and (D.19), the vector on the right-hand side of (D.17) tends to zero. Therefore, since Γ + n has full rank, the solution ς n (s) to the system (D.17 2022)). Further define c 1 (n) = . . . = c r 1 (n) = 0, and writing γ = α q γ q + . . . + α r γ r , define c (n) = α s * (n) for ∈ I 0 (cf. (B.49) in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]). In particular, the resulting vector c(n) = (c 1 (n), . . . , c r (n)) * of coefficients appearing in (D.10) is a unit vector. For the purpose of establishing the second statement in (D.11), we focus on studying the decay rate of the entries c (n), / ∈ I + . By construction, ϑ n = O(a -). In addition, for ∈ I 0 and i ∈ I + , p i (n), u p-r+ (n) = O(a h i -h ) = O(a -) due to Lemma D.1. Moreover, (D.20) further implies that the root (D.21) satisfies s * (n) → 1 as n → ∞. Thus, from (D.22), we obtain

Γ + n ς = ϑ n -s •
∈I + c (n) 2 ≤ (Γ + n ) -1 2 O(a -) + O(a -) 2 = O(a -2 ),
as was to be shown.

The following corollary and proposition are used in the proof of Proposition C.1. 2022)) Fix j ∈ N and suppose conditions (W 1 -W 4) and (A1 -A5) hold. Further assume that either (i) 0 < h 1 < . . . < h r < 1; or (ii) h 1 = . . . = h r and the functions ξ q (2 j ) in (B.1) satisfy q 1 = q 2 ⇒ ξ q 1 (1) = ξ q 2 (1). (D.24)

Let W(a(n)2 j ) be as given in Lemma D.1. Then, for each q ∈ {1, . . . , r}, there is a sequence of (p -r + q)-th unit eigenvectors {u p-r+q (n)} n∈N of W(a(n)2 j ) along which the limits lim n→∞ P * (n)u p-r+q (n) =: γ q (D.25)

exist.

E Conflicts of interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

F Data availability statement

The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

  Figure1: The high-dimensional asymptotic behavior of i . The plots above show the asymptotic behavior of local estimates q around each j as p and n j = n/2 j increase with their ratio fixed at p/n j = c, with c = 1/2 (left column) and c = 1/4 (right column). There are r = 6 scaling eigenvalues h i ∈ (0.1, 0.3, 0.5, 0.6, 0.8, 0.9). For each j, the reported values i are based on the octave range (j 1 , j 2 ) = (j -1, j + 1). Note the i associated with non-scaling eigenvalues stays below 0 as j increases in all instances, and the q associated with the r = 6 scaling eigenvalues approach their theoretical values. With c fixed (along each column), the larger sample size n (bottom row) results in near-convergence at smaller j (note the axes have been shifted to align the curves). Likewise, as the c decreases (between left and right columns), the near-convergence also occurs at smaller octaves j. Compared to the wavelet log-eigenvalues themselves (not shown), the multiscale estimator i displays reduced bias.

  Figure2: The asymptotic (joint) distribution of { h q } q=1,...,r = { p-r+q } q=1,...,r in high dimen-

  . (C.39) Expression (3.11) is now a consequence of (C.38) and (C.39).

Lemma D. 2 (

 2 [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF], Lemma B.4) For each n, let p (n), = r + 1, . . . , p, (D.4)

γ

  (n) := P * (n)u p-r+ (n), = 1, . . . , r.(D.6)Suppose n ∈ N ⊆ N is a subsequence along which the limitslim n →∞ γ (n ) =: γ , = 1, . . . , r, (D.7)exist. Fix any γ ∈ span{γ q , . . . , γ r-r 3 } (D.8)and any vector x * = (x * ,r-r 3 +1 , . . . , x * ,r ) * ∈ R r 3 . (D.9)Then, we can pick a sequence of unit vectorsv(n ) := r =q c (n)u p-r+ (n ) ∈ span{u p-r+q (n ), . . . , u p (n )} (D.10) satisfying P * (n )v(n ) → γ, n → ∞, ∈I + c 2 (n ) = O(a(n ) -2 ). (D.11) In (D.11), p (n ), v(n ) = x * , a(n ) h -hq , ∈ I + .(D.12)Proof: We seek only to establish the second statement in (D.11) since all other statements follow from the same argument (cf.Abry et al. (2022), Lemma B.6). For notational simplicity, write n = n . Also defineϑ n = x * , a(n) h -hq ∈I + ∈ R r 3 , (D.13)which is the vector that contains the right-hand terms in (D.12). Note thatϑ n → 0, n → ∞. (D.14)Consider the matrices Γ n ∈ M (r, R) and Γ + n ∈ M (r 3 , R) given byΓ n = p (n), u p-r+ (n) 1≤ , ≤r , Γ + n = p (n), u p-r+(n) , ∈I + . (D.15) By (D.7), as n → ∞, Γ n → Γ := (γ 1 , . . . , γ r ). (D.16) In particular, Γ + n → Γ + = Γ , , ∈I + . By Lemma B.5 in Abry et al. (

r-r 3 =qα

 3 p r-r 3 +1 (n), u p-r+ (n) , . . . , r-r 3 =q α p r (n), u p-r+ (n) * (D.17) has a unique solution ς = ς n (s) = ς r-r 3 +1,n (s), . . . , ς r,n (s) * ∈ R r 3 . (D.18) Note that, by Lemma D.1, as n → ∞, p (n), u p-r+ (n) → 0, ∈ I + , ∈ I 0 . (D.19)

  ς n (s) → 0, n → ∞.(D.20)Thus, for any small ε > 0, 0 ≤ sup s∈[0,1] ς n (s) < ε for large enough n. Now note that, for each n, the functions → f (s) = 1 -s 2 -ς n (s) 2 depends continuously on s. Moreover, f (0) > 1 -ε and f (1) = -ς n (1) 2 . Hence, f must have a root s * (n) ∈ (0, 1] (D.21)with probability tending to 1. So, for one such root s * (n) ∈ (0, 1], we can use expression (D.17) to define the coefficients c (n), ∈ I + , by means of   c r-r 3 +1 (n) . . . c r (n)    = ς n (s * (n)) = (Γ + n ) -1   ϑn -s * (n) •    r-r 3 =q α p r-r 3 +1 (n), u p-r+ (n) . . .r-r 3 =q α p r (n), u p-r+ (n) expression (B.51) in Abry et al. (

  Corollary D.1 (a consequence of Corollary B.1 in Abry et al. (2022)) Let W(a(n)2 j ) be as given in Lemma D.1. Then, λ p-r+W(a(n)2 j ) a(n) 2h +1→ ξ (2 j ), = 1, . . . , r. In particular,λ p-r+ W(a(n)2 j ) a(n) 2hq+1 → 23), ξ q (2 j) are the functions (B.1) appearing in Theorem B.1. Proposition D.1 (a consequence of Proposition B.1 in Abry et al. (

  R), let vec S (S) = (s 11 , s 21 , . . . , s n1 , s 22 , s 32 , . . . , s n2 , . . . , s nn ).

						i , p/n j 2 = 1/2							i , p/n j 2 = 1/4
		2															2
				n=256	n=1024	n=4096	n=16384	n=65536	n=262144	n=1048576		n=256	n=1024	n=4096	n=16384	n=65536	n=262144	n=1048576
				j=4	j=5	j=6		j=7		j=8		j=9		j=10		j=3	j=4	j=5	j=6	j=7	j=8	j=9
		1.5															1.5
		1															1
	i
	local	0 0.5															0 0.5
		-0.5															-0.5
		-1															-1
		p/n j -1.5	2	8/16	16/32	32/64		64/128	128/256	256/512	512/1024	p/n j -1.5	2	8/32	16/64	32/128	64/256	128/512	256/1024 512/2048
						i , p/n j 2 = 1/2							i , p/n j 2 = 1/4
		2
				n=64	n=256	n=1024	n=4096	n=16384	n=65536	n=262144	n=1048576
				j=2	j=3	j=4	j=5		j=6		j=7		j=8			j=9
		1.5
		1
	i
	local	0 0.5
		-0.5
		-1
		p/n j -1.5	2	8/16	16/32	32/64	64/128 128/256 256/512 512/1024 1024/2048	p/n j -1.5	2	8/32	16/64	32/128	64/256 128/512 256/1024 512/2048 1024/4096
																		(2.2)
	We use the asymptotic notation
														o P (1), O P (1)	(2.3)

* P.A. and H.W. were partially supported by ANR-16-CE33-0020 MultiFracs, France. H.W. was also partially supported by ANR-18-CE45-0007 MUTATION. G.D.'s long term visits to ENS de Lyon were supported by the school, the CNRS and the Simons Foundation collaboration grant #714014. The authors also gratefully acknowledge the support and resources from the Center for High Performance Computing at the University of Utah as well as the high performance computing (HPC) resources and services provided by Technology Services at Tulane University.

A Assumptions on the wavelet multiresolution analysis

In the main results of the paper, we make use of the following conditions on the underlying wavelet MRA.

Assumption (W 1): ψ ∈ L 2 (R) is a wavelet function, namely, it satisfies the relations R ψ 2 (t)dt = 1, R t p ψ(t)dt = 0, p = 0, 1, . . . , N ψ -1, R t N ψ ψ(t)dt = 0, (A.1) for some integer (number of vanishing moments) N ψ ≥ 1.

Assumption (W 2): the scaling and wavelet functions

and φ(0) = 1.

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

Conditions (A.1) and (A.2) imply that ψ(x) exists, is everywhere differentiable and its first N ψ -1 derivatives are zero at x = 0. Condition (A.3), in turn, implies that ψ is continuous (see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], Theorem 6.1) and, hence, bounded.

Note that assumptions (W 1 -W 4) are closely related to the broad wavelet framework for the analysis of κ-th order (κ ∈ N ∪ {0}) stationary increment stochastic processes laid out in Moulines et al. (2007a[START_REF] Moulines | On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter[END_REF][START_REF] Moulines | A wavelet Whittle estimator of the memory parameter of a nonstationary Gaussian time series[END_REF] and [START_REF] Roueff | Asymptotic normality of wavelet estimators of the memory parameter for linear processes[END_REF]. The Daubechies scaling and wavelet functions generally satisfy (W 1 -W 4) (see [START_REF] Moulines | A wavelet Whittle estimator of the memory parameter of a nonstationary Gaussian time series[END_REF]), p. 1927, or Mallat (1999), p. 253). Usually, the parameter α increases to infinity as N ψ goes to infinity (see [START_REF] Moulines | A wavelet Whittle estimator of the memory parameter of a nonstationary Gaussian time series[END_REF]), p. 1927[START_REF] Cohen | Numerical Analysis of Wavelet Methods[END_REF], Theorem 2.10.1).

B The asymptotic and large scale behavior of the eigenvalues of wavelet random matrices

In this section, we recap the main results in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]. It was established that, after proper rescaling, the r largest eigenvalues of a wavelet random matrix W(a(n)2 j ) in high dimensions converge in probability to deterministic functions ξ q (2 j ), q = 1, . . . , r. Thus, such functions can be interpreted as asymptotic rescaled eigenvalues. Notably, they display a scaling property. Moreover, the remaining p(n) -r eigenvalues of a wavelet random matrix are bounded in probability.

Theorem B.1 (a consequence of Theorem 3.1 in [START_REF] Abry | On high-dimensional wavelet eigenanalysis[END_REF]) Assume (W 1 -W 4) and (A1 -A5) hold, and fix j ∈ N. Then, for p = p(n), the limits p-lim

Consider the first and the second terms inside the max{•, •} on the right-hand side of expression (C.17). If (a)

then (C.5) holds under condition (i). So, we now establish (a) and (b). First, we show (a). Let γ q (n) and x * (n) be as in (C.8) and (C.10), respectively. For notational simplicity, write where lim n→∞ p q (n), u p-r+q (n) = γ qq as a consequence of (C.9). Recast

Therefore, in view of condition (2.21), if

Thus, (C.18) holds, which establishes (a). So, we now show (C.23). We establish relation (C.23) entry-wise for each of the ranges ∈ I -= {1, . . . , q -1}, = q and ∈ I + = {q + 1, . . . , r}.

In fact, for ∈ I -= {1, . . . , q -1}, expression (C.20) shows that

On the other hand, for = q, let P Consequently, after flipping the sign of u p-r+q (n) if necessary, Q * (n)u p-r+q (n) -e q = O(a -). Recalling (C.8) and the condition that R(n) → R as n → ∞ (see (2.25)), γ q (n) -Re q = R * (n)Q * (n)u p-r+q (n) -Re q