# On the empirical spectral distribution of large wavelet random matrices based on mixed-Gaussian fractional measurements in moderately high dimensions

Patrice Abry, Gustavo Didier, Oliver Orejola, Herwig Wendt

<figure><img src="image_1.png" /></figure>

## To cite this version:

Patrice Abry, Gustavo Didier, Oliver Orejola, Herwig Wendt. On the empirical spectral distribution of large wavelet random matrices based on mixed-Gaussian fractional measurements in moderately high dimensions. Electronic Journal of Probability, 2024, 30 (none), 10.1214/25-EJP1297. hal-05354835

HAL Id: hal-05354835

https://hal.science/hal-05354835v1

Submitted on 8 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Electron. J. Probab. 30 (2025), article no. 40, 1–48.
ISSN: 1083-6489 https://doi.org/10.1214/25-EJP1297

# On the empirical spectral distribution of large wavelet random matrices based on mixed-Gaussian fractional measurements in moderately high dimensions*

Patrice Abry† Gustavo Didier‡ Oliver Orejola§ Herwig Wendt¶

## Abstract

In this paper, we characterize the convergence of the (rescaled logarithmic) empirical spectral distribution of wavelet random matrices. We assume a moderately high-dimensional framework where the sample size $n$, the dimension $p(n)$ and, for a fixed integer $j$, the scale $a(n)2^j$ go to infinity in such a way that $\lim_{n \to \infty} p(n) \cdot a(n)/n = \lim_{n \to \infty} o(\sqrt{a(n)/n}) = 0$. We suppose the underlying measurement process is a random scrambling of a sample of size $n$ of a growing number $p(n)$ of fractional processes. Each of the latter processes is a fractional Brownian motion conditionally on a randomly chosen Hurst exponent. We show that the (rescaled logarithmic) empirical spectral distribution of the wavelet random matrices converges weakly, in probability, to the distribution of Hurst exponents.

**Keywords:** wavelets; random matrices; bulk; mixed-Gaussian distribution.

**MSC2020 subject classifications:** 60G18; 60B20; 42C40.

Submitted to EJP on January 5, 2024, final version accepted on February 17, 2025.

## 1 Introduction

A wavelet is a unit $L^2(\mathbb{R})$-norm function that annihilates polynomials (see (A.1)). For a fixed (octave) $j \in \mathbb{N} \cup \{0\}$, a wavelet random matrix is given by

$$ \mathbb{R}^{p^2} \ni \mathbf{W}_n \equiv \mathbf{W}(a(n)2^j) = \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} \mathcal{D}(a(n)2^j, k) \mathcal{D}(a(n)2^j, k)^*. \quad (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$. Each random vector $\mathcal{D}(a(n)2^j, k) \in \mathbb{R}^p$ is the wavelet

<footnote>*H.W. was partially supported by ANR-18-CE45-0007 MUTATION, France. G.D.'s long term visits to ENS de Lyon were supported by the school, the CNRS and the Simons Foundation collaboration grant #714014.</footnote>

<footnote>†CNRS, ENS de Lyon, Laboratoire de Physique, France. E-mail: patrice.abry@ens-lyon.fr</footnote>

<footnote>‡Tulane University, U.S.A. E-mail: gdidier@tulane.edu</footnote>

<footnote>§Tulane University, U.S.A. E-mail: oorejola@gmail.com</footnote>

<footnote>¶IRIT-ENSEEIHT, CNRS (UMR 5505), Université de Toulouse, France. E-mail: herwig.wendt@irit.fr</footnote>
transform of a p-variate stochastic process Y at dyadic scale $a(n)2^j$ and shift $k \in \mathbb{Z}$ (see (2.12) and (2.19) for the continuous- and discrete-time definitions, respectively). The entries of $\{D(a(n)2^j, k)\}_{k=1,...,n_{a,j}}$ are generally correlated. A fractal is an object or phenomenon that displays the property of self-similarity, in some sense, across a range of scales (Mandelbrot [62]). Due to its intrinsic multiscale character and fine-tuned mathematical properties, the wavelet transform has been widely used in the study of fractals (e.g., Wornell [106], Doukhan et al. [39], Massopust [65]). In this paper, we characterize the asymptotic behavior of the (rescaled logarithmic) empirical spectral distribution of large wavelet random matrices in a moderately high-dimensional framework where the sample size $n$, the dimension $p(n)$ and the scale $a(n)2^j$ go to infinity. We assume measurements of the form

$$ \mathbb{R}^p \ni Y(t) = \mathbf{P}(n)X(t), \qquad (1.2) $$

either in continuous or in discrete time $t$. In (1.2), the coordinates matrix $\mathbf{P}(n)$ is random and independent of $X(t)$. Moreover, each row of $X$ is, conditionally, an independent $H$-fractional Brownian motion, where each Hurst exponent $H$ is picked independently from a discrete probability distribution $\pi(dH)$. In the main results (Theorems 3.1 and 3.3), we show that the (rescaled logarithmic) empirical spectral distribution of $\mathcal{W}(a(n)2^j)$ converges weakly, in probability, to a known affine transformation of the distribution $\pi(dH)$ of the Hurst exponents.

In this paper, we combine two mathematical frameworks that are rarely considered jointly: (i) high-dimensional probability theory; (ii) fractal analysis. This is done by bringing together the study of large random matrices and scaling analysis in the wavelet domain.

Since the 1950s, the spectral behavior of large-dimensional random matrices has attracted considerable attention from the mathematical research community. For example, random matrices have proven to be prolific statistical mechanical models of Hamiltonian operators (e.g., Mehta and Gaudin [69], Dyson [40], Ben Arous and Guionnet [20], Soshnikov [94], Mehta [70], Anderson et al. [9], Erdős et al. [44]). They are also of great interest in combinatorics and numerical analysis (e.g., Baik et al. [14], Baik et al. [15], Menon and Trogdon [71]), as well as in the study of integrable systems and universality (Kuijlaars and McLaughlin [56], Baik et al. [16], Deift [36], Tao and Vu [98], Borodin and Petrov [23], Deift [37]). In particular, the literature on random matrices under dependence has been expanding at a fast pace (e.g., Xia et al. [107], Paul and Aue [82], Che [29], Steland and von Sachs [95], Wang et al. [103], Zhang and Wu [109], Erdős et al. [43], Merlevède et al. [74], Bourguin et al. [24]).

In turn, recall that the emergence of a fractal is typically the signature of a physical mechanism that generates scale invariance (e.g., Mandelbrot [62], West et al. [104], Zheng et al. [111], He [48], Shen et al. [90]). Unlike traditional statistical mechanical systems (e.g., Reif [86]), a scale-invariant system does not display a characteristic scale, namely, one that dominates its statistical behavior. Instead, the behavior of the system across scales is determined by specific parameters called scaling exponents. Scale invariance manifests itself in a wide range of natural and social phenomena such as in criticality (Sornette [93]), turbulence (Kolmogorov [55]), climate studies (Isotta et al. [52]), dendrochronology (Bai and Taqqu [11]) and hydrology (Benson et al. [21]). Mathematically, it is a topic of central importance in Markovian settings (e.g., diffusion, lattice models, universality classes) as well as in non-Markovian ones (e.g., anomalous diffusion, long-range dependence, non-central limit theorems).

In a multidimensional framework, scaling behavior does not always appear along standard coordinate axes, and often involves multiple scaling relations. A $\mathbb{R}^r$-valued stochastic process $X$ is called operator self-similar (o.s.s.; Laha and Rohatgi [57], Hudson
and Mason [51]) if it exhibits the scaling property

$$ \{X(ct)\}_{t \in \mathbb{R}} \stackrel{\text{f.d.d.}}{=} \{c^{\mathbf{H}}X(t)\}_{t \in \mathbb{R}}, \quad c > 0. \qquad (1.3) $$

In (1.3), we define $c^{\mathbf{H}} = \exp\{\log(c)\mathbf{H}\} = \sum_{k=0}^{\infty} \frac{(\log(c)\mathbf{H})^k}{k!}$, where, in Jordan form, the $r \times r$ (Hurst) matrix

$$ \mathbf{H} = \mathbf{P}\mathbf{J}_{\mathbf{H}}\mathbf{P}^{-1} \qquad (1.4) $$

only displays eigenvalues with strictly positive real parts. In particular, if $\mathbf{P}$ as in (1.4) is non-diagonal, the scaling behavior described in (1.3) involves non-standard coordinates. A canonical model for multivariate fractional systems is operator fractional Brownian motion (ofBm), namely, a Gaussian, o.s.s., stationary-increment stochastic process (Maejima and Mason [60], Mason and Xiao [64], Didier and Pipiras [38]). In particular, ofBm is the natural multivariate generalization of the classical fractional Brownian motion (fBm; Embrechts and Maejima [41]).

In applications, the importance of the role of multiple scaling laws along non-standard coordinates is now well established. For example, in econometrics, the detection of distinct scaling laws in multivariate fractional time series is indicative of the key property of cointegration – namely, the existence of meaningful and statistically useful long-run relationships among the individual series (e.g., Engle and Granger [42], NobelPrize.org [80], Hualde and Robinson [50], Shimotsu [91]). From a different perspective, it has been shown that ignoring the presence of multiple scaling laws in statistical inference may lead to severe biases (the so-called dominance and amplitude effects – see, for instance, Abry and Didier [4]).

In the study of scaling properties, eigenanalysis has proven to be a fecund analytic framework in a number of areas, such as in the characterization of heavy tails and correlation (e.g., Meerschaert and Scheffler [67, 68], Becker-Kern and Pap [19]) and in time series modeling (e.g., Phillips and Oullaris [84], Li et al. [59], Zhang et al. [110]). In Abry and Didier [3, 4], wavelet eigenanalysis is put forward in the construction of a general methodology for the statistical identification of the scaling (Hurst) structure of ofBm in low dimensions (cf. (1.4)).

To the best of our knowledge, the behavior of the eigenvalues of large-dimensional wavelet random matrices was mathematically studied for the first time in Abry et al. [1, 2]. This was done in the context of a high-dimensional signal-plus-noise model where measurements display a fixed number $r$ of scaling laws, each driven by a deterministic (and unknown) Hurst exponent. In those papers, the scaling, convergence and joint fluctuations of the top $r$ eigenvalues of the associated wavelet random matrix were established and characterized.

As opposed to the “strong” approach in Abry et al. [1, 2], where individual eigenvalues are tracked, in this paper we put forward a weak (or “Wigner-like”) approach. In other words, instead of assuming a fixed number of deterministic scaling laws (Hurst exponents), we suppose that, as the dimension grows, measurements display an increasing number of random scaling laws. This provides a natural model for multiparameter high-dimensional fractal systems that emerge in a wide range of fields of application. These fields include neuroscience and fMRI imaging (Li et al. [58], Gotts et al. [46]), network traffic (Abry and Didier [3]), climate science (Schmith et al. [89]), high-dimensional time series (Merlevède and Peligrad [75], Chan et al. [27], Alshammri and Pan [8]) as well as blind source separation problems and principal components analysis (Comon and Jutten [32], Chang et al. [28], Abry et al. [5]). In this context, it is of primary interest to understand the bulk behavior of the empirical spectral distribution (e.s.d.) of wavelet random matrices.
More specifically, following up on the computational study Orejola et al. [81], for each $n \in \mathbb{N}$ we assume a random vector

$$H_1, \dots, H_{p(n)} \in (0, 1)^{p(n)} \qquad (1.5)$$

is sampled from a discrete probability measure $\pi(dH)$ on $\mathcal{B}(0, 1)$. Then, measurements of the form (1.2) are made, where $X$ displays independent rows and the $\ell$-th row of $X$ is, conditionally on (1.5), a fBm with Hurst parameter $H_\ell$ (on the recent use of this class of processes in the modeling of anomalous diffusion, see Balcerek et al. [17]). The presence of the generally non-diagonal coordinates matrix $\mathbf{P}(n)$ implies that the scaling behavior of the measurement process (1.2) does not follow standard coordinate axes (cf. the matrix $\mathbf{P}$ in (1.4)). In addition, the marginal distributions of the measurement process (1.2) are random sums of mixed-Gaussian laws.

In the main results of this paper (Theorems 3.1 and 3.3), we describe the asymptotic behavior of the rescaled logarithmic e.s.d. of the wavelet random matrix

$$\mathcal{W}_n = \mathbf{W}(a(n)2^j) \text{ or } \mathcal{W}_n = \widetilde{\mathbf{W}}(a(n)2^j), \qquad (1.6)$$

corresponding to continuous- and discrete-time measurements, respectively. We consider both measurement frameworks because, in the former case, exact self-similarity-type relations hold (e.g., (C.65)), which is very mathematically convenient. Once results are obtained for continuous-time measurements, then we turn to the realistic situation where measurements are made in discrete time. When taking limits, we consider a moderately high-dimensional regime where

$$\frac{p(n) \cdot a(n)}{n} = o\left(\sqrt{\frac{a(n)}{n}}\right), \quad \lim_{n \to \infty} \frac{a(n)}{n} = 0 \qquad (1.7)$$

(see (2.7)). Namely, the sample size $n$, the dimension $p(n)$ and the scale $a(n)2^j$ go to infinity in a way that $n$ grows faster than $p(n) \cdot a(n)$. Over fixed scales (i.e., when $a(n)$ is assumed constant in (1.7)), moderately high-dimensional regimes have been studied in several works. These include Bai and Yin [12], Jiang [53] and Wang et al. [103] on Wigner, Jacobi and sample covariance matrices, respectively (on the use of condition (1.7) in this paper, see the discussion in Section 4).

It is by considering the three-way limit (1.7), which includes the scaling limit $a(n) \to \infty$, that large wavelet random matrices may be used in the characterization of low-frequency behavior in a (moderately) high-dimensional framework. In fact, let $F_{\mathcal{W}_n}(dv)$ be the rescaled logarithmic e.s.d. of $\mathcal{W}_n$ as in (1.6) (see (3.2) and (3.25) for the precise definition of $F_{\mathcal{W}_n}(dv)$ for continuous- and discrete-time measurements, respectively). In Theorems 3.1 and 3.3, we show that, for any appropriate test function $g$,

$$\int_{\mathbb{R}} g(v) F_{\mathcal{W}_n}(dv) \xrightarrow{P} \int_{\mathbb{R}} g(v) F_{2\mathcal{H}+1}(dv), \quad n \to \infty. \qquad (1.8)$$

In other words, $F_{\mathcal{W}_n}(dv)$ converges weakly, in probability, to the distribution given by the random variable $2\mathcal{H} + 1$. Namely, it converges to the distribution of Hurst exponents up to a known affine transformation (see Figure 1; see also Remark 3.2 and Figure 2 on a schematic representation of the “forces” acting on rescaled wavelet log-eigenvalues).

Note that (1.7) and (1.8) stand in sharp contrast with the traditional analysis of large random matrices, which involves what we call fixed scales – i.e., the scaling factor $a(n)$ as in (1.7) is held constant. With standard sample covariance matrices, for example, one considers the ratio $\lim_{n \to \infty} p(n)/n$ and the e.s.d. often converges to a Dirac mass at a constant or to a Marčenko-Pastur distribution in moderately high and in high dimensions, respectively (e.g., Chen and Pan [30], Bai and Silverstein [13], Tao and Vu [99]). In the
<figure><img src="image_4.png" /><figcaption>Figure 1: **The distribution of the (rescaled logarithmic) wavelet e.s.d. in the three-way limit (1.7).** A Monte Carlo study displays a tri-modal distribution emerging in the rescaled logarithmic wavelet e.s.d. in the three-way limit (1.7) (**n.b.:** after applying an affine transformation, the results are shown on the same scale as that of the distribution $\pi(dH)$). In the depicted simulation study based on 1000 realizations, $\pi(dH)$ is a discrete uniform distribution with support $\{0.2, 0.5, 0.8\}$. From left to right, respectively, (sample size, scale, dimension) = ($2^{10}$, $2^4$, $2^3$), ($2^{15}$, $2^5$, $2^5$) and ($2^{18}$, $2^6$, $2^6$). Wavelet log-eigenvalues weighted over multiple scales were used for enhanced ("debiased") finite-sample convergence (cf. Abry and Didier [3] and Abry et al. [1]).</figcaption></figure>

presence of strongly dependent measurements, important results on both bulk and edge behavior for fixed scales were established in several works (e.g., Chakrabarty et al. [26], Merlevède and Peligrad [75], Yaskov [108], Merlevède et al. [74], Tian [101]).

From the standpoint of probability theory, to the best of our knowledge this paper provides the first mathematical study of the asymptotic behavior of the e.s.d. of wavelet random matrices starting from measurements containing $p(n) \rightarrow \infty$ scaling laws. In addition, statement (1.8) suggests a new form of universality that is analogous to the law of large numbers (Tao [97]). In other words, regardless of the underlying $\pi(dH)$, the wavelet eigenspectrum reveals the true distribution $\pi(dH)$ of the Hurst modes in the moderately high-dimensional limit (1.7).

From the standpoint of fractal analysis, this paper takes a decisive step in the expansion, to the high-dimensional context, of the study of scale-invariant and non-Markovian phenomena started by Kolmogorov [54] and Mandelbrot and Van Ness [63], and later taken up by the likes of Flandrin [45], Wornell and Oppenheim [105], Meyer et al. [76], among many others (see Pipiras and Taqqu [85] and references therein).

This paper is organized as follows. In Section 2, we provide the basic wavelet framework, definitions and assumptions used throughout the paper. In Section 3, we state and establish the main results on the behavior of the rescaled logarithmic e.s.d. of wavelet random matrices in the three-way limit (1.7). The proofs of Theorems 3.1 and 3.3 are written so as to provide readers with immediate access to the core ideas behind the argument, whereas the auxiliary technical results can be found in the Appendices. In Section 4, we lay out conclusions. We also briefly discuss several open problems that this work leads to in the theory and applications of wavelet random matrices. With a view to universality claims, these include potential ways of expanding the mathematical framework of the paper.
<figure><img src="image_5.png" /><figcaption>Figure 2: Schematic representation of the “forces” acting on the rescaled log-eigenvalues in the three-way limit (1.7). We can interpret that two distinct “fractal forces” – one repulsive, one attractive – act on the rescaled wavelet log-eigenvalues. In the three-way limit (2.7), as a result of these forces only the large-scale features of the system remain. We are then left with the distribution of the Hurst exponents, as illustrated in Figure 1 (see Remark 3.2 for a discussion).</figcaption></figure>

# 2 Preliminaries and assumptions

## 2.1 Notation

Throughout the paper, we use the following notation. All with respect to the field $\mathbb{R}$, $M(m_1, m_2, \mathbb{R})$ and $M(n) = M(n, \mathbb{R})$ are the vector spaces, respectively, of all $m_1 \times m_2$ matrices and $n \times n$ matrices, $GL(n, \mathbb{R})$ is the general linear group (invertible matrices), and $O(n)$ is the orthogonal group (i.e., matrices $O \in M(n)$ such that $OO^* = I$). Also, $S(n, \mathbb{R})$ and $S_{\ge 0}(n, \mathbb{R})$ denote, respectively, the sets of $n \times n$ symmetric and symmetric positive semidefinite matrices, whereas $I = I_n$ denotes the $n \times n$ identity matrix. The $(n-1)$-dimensional unit sphere is represented by the symbol $S^{n-1}$. For any $A \in S(n, \mathbb{R})$, the ordered eigenvalues of $A$ are denoted $\lambda_1(A) \le \dots \le \lambda_n(A)$. More generally, for $M \in M(n)$, the ordered singular values of $M$ are written

$$ \sigma_1(\mathbf{M}) \le \dots \le \sigma_n(\mathbf{M}). \qquad (2.1) $$

For any $A \in M(n)$, $\|A\|_{op} := \sigma_n(A)$ is the operator norm of $A$. The symbol $\text{Haar}(O(n))$ represents the Haar probability measure on $O(n)$. The symbol $o_P(1)$ denotes a random variable that vanishes in probability as $n \to \infty$. Given a collection of scalars $x_1, \dots, x_n$,

$$ x_{(1)} \le \dots \le x_{(n)} \qquad (2.2) $$

denotes the associated ordered $n$-tuple.

## 2.2 Measurements

In regard to the underlying stochastic framework, first consider the following definition.

**Definition 2.1.** Let $\pi(dH)$ be a probability measure such that $\pi(0, 1) = 1$. The univariate stochastic process

$$ X_H = \{X_H(t)\}_{t \in \mathbb{R}} \qquad (2.3) $$
is called a random Hurst (exponent)-fractional Brownian motion (rH-fBm) when, conditionally on some value $\mathcal{H} = H$ picked from $\pi(dH)$, $X_{\mathcal{H}}$ is a standard fBm with Hurst exponent

$$H \in (0, 1). \qquad (2.4)$$

The basic properties of rH-fBm are provided in Lemma C.12.

So, throughout this manuscript, we make use of the following assumptions on the measurements. For expository purposes, we first state the assumptions, and then provide some interpretation. In the assumptions, $\mathcal{T}$ denotes either $\mathcal{T} = \mathbb{R}$ or $\mathcal{T} = \{1, \dots, n\}$, corresponding to the cases of continuous- or discrete-time measurements, respectively.

**Assumption (A1).** $\pi(dH)$ is a distribution such that

$$\mathrm{supp} \, \pi(dH) = \{\check{H}_1, \dots, \check{H}_r\}, \quad \varpi := \min\{\mathrm{supp} \, \pi(dH)\}, \qquad (2.5)$$

where $r \in \mathbb{N}$ and $0 < \check{H}_1 < \dots < \check{H}_r < 1$.

**Assumption (A2).** For each $n \in \mathbb{N}$, $p = p(n)$ independent rH-fBm sample paths are generated based on (independently) sampling from $\pi(dH)$ as in (2.5). In particular, when restricted to discrete time, a total of $p \times n$ entries is available. The associated $p$-variate (latent) stochastic process is denoted by $\{X(t)\}_{t \in \mathcal{T}}$.

**Assumption (A3).** For $p = p(n)$ and $\{X(t)\}_{t \in \mathcal{T}}$ as in assumption (A2), the measurements have the form

$$\mathbb{R}^{p(n)} \ni Y(t) = \mathbf{P}(n)X(t), \quad t \in \mathcal{T}. \qquad (2.6)$$

In (2.6), the so-named coordinates matrix $\mathbf{P}(n) \in GL(p(n), \mathbb{R})$ that is independent of $X(t)$.

**Assumption (A4).** Fix $j \in \mathbb{N} \cup \{0\}$. The dimension $p(n)$ and the scaling factor $a(n)$ (cf. $n_{a,j}$ as in (2.15) below) satisfy the relations

$$a(n) \le \frac{n}{2^j}, \quad \frac{a(n)}{n} + \frac{n}{a(n)^{1+2\varpi}} \to 0, \quad p(n) < \frac{n}{a(n)2^j}, \quad \infty \leftarrow p(n) = o\left(\sqrt{\frac{n}{a(n)}}\right), \qquad (2.7)$$

as $n \to \infty$, where $\varpi$ is as in (2.5).

**Assumption (A5).** For the random matrix $\mathbf{P}(n)$ as in (2.6),

$$\frac{\log \sigma_1(\mathbf{P}(n))}{\log a(n)} \xrightarrow{\mathbb{P}} 0, \quad \frac{\log \sigma_p(\mathbf{P}(n))}{\log a(n)} \xrightarrow{\mathbb{P}} 0, \quad n \to \infty. \qquad (2.8)$$

Assumption (A1) defines the distribution of Hurst exponents. Assumption (A2) postulates the latent $p(n)$-variate stochastic process $X$ as a collection of independent rH-fBms. Assumption (A3) describes the observed process $Y$, where the unknown random coordinates matrix $\mathbf{P}(n)$ determines the directions of the multiple (random) scaling relations stemming from the latent process $X$. In turn, assumption (A4) controls the divergence rates among $n$, $a(n)$ and $p(n)$ in the three-way limit. In particular, it states that the scaling factor $a(n)$ and the number of fractional processes, $p(n)$, must grow slower than $n$, and that the three-component ratio $\frac{p(n)a(n)}{n}$ must converge to 0 (see (1.7)). This establishes a moderately high-dimensional regime (cf. the traditional ratio $\lim_{n \to \infty} \frac{p(n)}{n}$ for sample covariance matrices). By assumption (A5), the minimum and maximum singular values of $\mathbf{P}(n)$ neither vanish nor diverge too fast, respectively. Hence, $\mathbf{P}(n)$ cannot interfere too heavily on the scaling relations emanating from the latent process.

Throughout this manuscript, for each $n \in \mathbb{N}$ it will be useful to write

$$\mathbb{H}_n = \mathrm{diag}(H_1, \dots, H_{p(n)}) \in \mathcal{M}(p(n)) \qquad (2.9)$$
<table><caption>Table 1: Different uses of the letter “H” throughout the paper.</caption><thead><tr><th>notation</th><th>domain</th><th>description</th><th>defined in</th></tr></thead><tbody><tr><td><b>H</b></td><td><i>M</i>(<i>p</i>)</td><td>Hurst matrix</td><td>(1.3)</td></tr><tr><td><b>H</b><sub>n</sub></td><td><i>M</i>(<i>p</i>)</td><td>diagonal matrix whose main diagonal entries are picked from &pi;(<i>dH</i>)</td><td>(2.9)</td></tr><tr><td><b>H&#x0302;</b><sub>n</sub></td><td><i>M</i>(<i>p</i>)</td><td>particular (deterministic) instance of <b>H</b><sub>n</sub></td><td>(B.37)</td></tr><tr><td><b>H</b></td><td>(0, 1)</td><td>random Hurst exponent</td><td>(2.3)</td></tr><tr><td><b>H</b></td><td>(0, 1)</td><td>particular (deterministic) instance of <b>H</b></td><td>(2.4)</td></tr><tr><td><b>H&#x0302;</b></td><td>(0, 1)</td><td>value in supp &pi;(<i>dH</i>)</td><td>(2.5)</td></tr></tbody></table>

to denote a diagonal random matrix whose (main diagonal) entries correspond to the sampled Hurst exponents as described in assumption (A2). Also, whenever convenient we write

$$a = a(n) \text{ or } p = p(n). \tag{2.10}$$

**Remark 2.2.** Due to the many uses of the letter “H” throughout the paper, for the readers’ convenience we provide Table 1 to help them keep track of the notation.

## 2.3 Wavelet multiresolution analysis

Recall that a wavelet $\psi$ is a unit $L^2(\mathbb{R})$-norm function that annihilates polynomials (see (A.1)). Throughout the paper, we make use of a wavelet multiresolution analysis (MRA; see Mallat [61], chapter 7). A wavelet MRA decomposes $L^2(\mathbb{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, j \in \mathbb{Z}$. In particular, given a wavelet $\psi$, there is a related scaling function $\phi \in L^2(\mathbb{R})$. Appropriate rescalings and shifts of $\phi$ and $\psi$ form bases for the subspaces $V_j$ and $W_j$, respectively (see Mallat [61], Theorems 7.1 and 7.3).

In almost all mathematical statements, we make assumptions (W1)-(W3) on the underlying wavelet MRA. Such assumptions are standard in the wavelet literature and are accurately described in Section A.1. In particular, we make use of a compactly supported wavelet basis.

## 2.4 Wavelet random matrices: continuous time

For any $p \in \mathbb{N}$, suppose $\mathcal{T} = \mathbb{R}$ in (2.6). Namely, assume we observe the continuous-time vector-valued stochastic process

$$\{Y(t)\}_{t \in \mathbb{R}} \subseteq \mathbb{R}^p. \tag{2.11}$$

For a wavelet function $\psi$, the wavelet transform vector of the stochastic process Y is defined as the convolution

$$\mathbb{R}^p \ni D(2^j, k) := 2^{-j/2} \int_{\mathbb{R}} \psi(2^{-j}t - k)Y(t)dt, \quad j \in \mathbb{N} \cup \{0\}, \quad k \in \mathbb{Z}, \tag{2.12}$$

whenever such expression is meaningful. Likewise, let

$$\mathbb{R}^p \ni D_X(2^j, k) := 2^{-j/2} \int_{\mathbb{R}} \psi(2^{-j}t - k)X(t)dt \tag{2.13}$$

be the wavelet transform of the latent process $\{X(t)\}_{t \in \mathbb{R}} \subseteq \mathbb{R}^p$ associated with (2.11). Note that, under (A1)-(A5), (2.13) is well defined in the mean squared sense by Lemma C.14. Therefore, under the same conditions, (2.12) is also well defined in the mean squared sense.
So, let

$$n \in \mathbb{N} \qquad (2.14)$$

be the wavelet-domain “sample size”, i.e., the number of wavelet coefficients available at scale $2^0 = 1$. Also, let $j \in \mathbb{N} \cup \{0\}$ and let $a(n)$ be a dyadic scaling factor. Starting from (2.12), for $p = p(n)$ the associated wavelet random matrices will be denoted throughout this manuscript by

$$S_{\ge 0}(p(n), \mathbb{R}) \ni \mathbf{W}(a(n)2^j) := \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} D(a(n)2^j, k)D(a(n)2^j, k)^*, \quad n_{a,j} := \frac{n}{a(n)2^j}. \qquad (2.15)$$

In (2.15), $n_{a,j}$ may be understood as the effective sample size, namely, the number of wavelet transform vectors $D(a(n)2^j, \cdot)$ available at scale $a(n)2^j$. For notational convenience, we assume $n_{a,j} \in \mathbb{N}$ (cf. the discussion around (A.9), for the case of discrete-time measurements).

Likewise, starting from (2.13), we can define the associated wavelet random matrix

$$S_{\ge 0}(p(n), \mathbb{R}) \ni \mathbf{W}_X(a(n)2^j) := \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} D_X(a(n)2^j, k)D_X(a(n)2^j, k)^*. \qquad (2.16)$$

Note that both (2.15) and (2.16) are well defined a.s. From expressions (2.12) and (2.16), we can conveniently recast the wavelet random matrix (2.15) for the measurements (2.11) in the form

$$\mathbf{W}(a(n)2^j) = \mathbf{P}(n)\mathbf{W}_X(a(n)2^j)\mathbf{P}(n)^*. \qquad (2.17)$$

## 2.5 Wavelet random matrices: discrete time

Now consider the realistic situation where $\mathcal{T} = \{1, \dots, n\}$ in (2.6). Namely, suppose we observe the discrete-time, vector-valued stochastic process

$$\{Y(k)\}_{k=1,\dots,n} \subseteq \mathbb{R}^p, \qquad (2.18)$$

associated with the starting scale $2^j = 1$ (or octave $j = 0$). In particular, $n \in \mathbb{N}$ in (2.18) denotes the sample size (cf. (2.14)). Starting from (2.18), we suppose the wavelet transform vector $\tilde{D}(2^j, k)$ of the high-dimensional process $Y$ stems from Mallat's pyramidal algorithm (Mallat [61], chapter 7). It is given by the convolution

$$\mathbb{R}^p \ni \tilde{D}(2^j, k) := \sum_{\ell \in \mathbb{Z}} Y(\ell)h_{j,2^j k-\ell}, \quad j \in \mathbb{N} \cup \{0\}, \qquad (2.19)$$

where we use the convention $Y(k) = 0$ for $k \notin \{1, \dots, n\}$ and the filter terms are defined by

$$\mathbb{R} \ni h_{j,\ell} := 2^{-j/2} \int_{\mathbb{R}} \phi(t+\ell)\psi(2^{-j}t)dt, \quad \ell \in \mathbb{Z}. \qquad (2.20)$$

Due to the assumed compactness of the supports of $\psi$ and of the associated scaling function $\phi$ (see condition (A.2)), only a finite number of filter terms is nonzero (Daubechies [35]). It follows that (2.19) is well defined a.s. A more detailed description of Mallat's algorithm is provided in Section A.2. For $n, j, a(n)$ and $n_{a,j}$ as in (2.12), the associated wavelet random matrices will be denoted throughout this manuscript by

$$S_{\ge 0}(p(n), \mathbb{R}) \ni \tilde{\mathbf{W}}(a(n)2^j) := \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} \tilde{D}(a(n)2^j, k)\tilde{D}(a(n)2^j, k)^*. \qquad (2.21)$$

Likewise, let

$$\mathbb{R}^{p(n)} \ni \tilde{D}_X(a(n)2^j, k) := \sum_{\ell \in \mathbb{Z}} X(\ell)h_{j,2^j k-\ell} \qquad (2.22)$$
be the discrete-time wavelet transformation of the latent process X. As with (2.19), (2.22) is well defined a.s. Thus, we can naturally define the associated wavelet random matrix

$$
\begin{align}
S_{\ge 0}(p(n), \mathbb{R}) \ni \tilde{\mathbf{W}}_X(a(n)2^j) &:= \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} \tilde{D}_X(a(n)2^j, k) \tilde{D}_X(a(n)2^j, k)^* \nonumber \\
&= \frac{1}{n_{a,j}} \tilde{D}_X(a(n)2^j) \tilde{D}_X(a(n)2^j)^*. \tag{2.23}
\end{align}
$$

In (2.23), each term $\tilde{D}_X(a(n)2^j) \in \mathcal{M}(p(n), n_{a,j}, \mathbb{R})$ is a random matrix with $k$-th column given by (2.22). From expressions (2.19) and (2.23), we can conveniently recast the wavelet random matrix (2.21) for the measurements (2.18) in the form

$$
\tilde{\mathbf{W}}(a(n)2^j) = \mathbf{P}(n)\tilde{\mathbf{W}}_X(a(n)2^j)\mathbf{P}(n)^*. \quad (2.24)
$$

3 Main results

As anticipated in the Introduction, in this section we state and prove the main
results. Theorems 3.1 and 3.3, pertaining to continuous- and discrete-time measurements,
respectively. For ease of understanding, the proofs of the theorems contain the bulk of
the arguments. The auxiliary technical results can be found in the Appendices.

**3.1 Continuous time**

In this section, we consider the framework in which measurements are made in
continuous time (see (2.11)). In particular, recall that, in this context, the wavelet
transform is given by (2.12).

Fix $j \in \mathbb{N} \cup \{0\}$ as in assumption (A4). For notational simplicity, write

$$
\mathbf{W}_n = \mathbf{W}(a(n)2^j). \tag{3.1}
$$

In the first main result of this paper, we describe the asymptotic behavior of the rescaled
log-eigenspectrum of the wavelet random matrix $\mathbf{W}_n$. So, to state the theorem, let

$$
\lambda_1(\mathbf{W}_n) \le \dots \le \lambda_{p(n)}(\mathbf{W}_n)
$$

be the ordered eigenvalues of the wavelet random matrix $\mathbf{W}_n$. Also, let

$$
F_{\mathbf{W}_n}(dv) := \frac{1}{p} \sum_{\ell=1}^{p} \delta \left\{ \frac{\log \lambda_{\ell}(\mathbf{w}_n)}{\log a(n)} \le \upsilon \right\} (dv) \quad (3.2)
$$

be the empirical spectral distribution (e.s.d.) of the rescaled log-eigenvalues of $\mathbf{W}_n$. In Theorem 3.1, stated and proven next, we show that $F_{\mathbf{W}_n}(dv)$ converges weakly, in probability, to the distribution $\pi(dH)$ of the Hurst exponents up to the affine transformation $x \mapsto 2x + 1$ of the latter.

**Theorem 3.1.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Then, the e.s.d. $F_{\mathbf{W}_n}$ as in (3.2) converges weakly, in probability, to $F_{2\mathcal{H}+1}$. In other words, for any $g \in C_b(\mathbb{R})$ and for any $\varepsilon > 0$,

$$
\lim_{n \to \infty} \mathbb{P} \left( \left| \int_{\mathbb{R}} g(v) F_{\mathbf{W}_n}(dv) - \int_{\mathbb{R}} g(v) F_{2\mathcal{H}+1}(dv) \right| > \varepsilon \right) = 0. \quad (3.3)
$$

Proof. Let

$$
F_{\mathbf{W}_n}(v) := \frac{1}{p} \sum_{\ell=1}^{p} \mathbf{1} \left\{ \frac{\log \lambda_{\ell}(\mathbf{w}_n)}{\log a(n)} \leq v \right\}, \quad v \in \mathbb{R}, \qquad (3.4)
$$
be the empirical cumulative distribution function (e.c.d.f.) associated with $F_{W_n}(dv)$ as in (3.2). By a natural adaptation of the proof of the Helly-Bray theorem (Shiryaev [92], pp. 377–378), it suffices to show that, for any point of continuity $v_0$ of $F_{2\mathcal{H}+1}$,

$$F_{W_n}(v_0) \xrightarrow{P} F_{2\mathcal{H}+1}(v_0), \quad n \to \infty. \qquad (3.5)$$

So, towards establishing (3.5), recall expressions (2.16) and (2.17) for $W_X(a(n)2^j)$ and $W(a(n)2^j)$, respectively. Now define

$$S_{\ge 0}(p(n), \mathbb{R}) \ni W_X(2^j) := \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} D_X(2^j, k) D_X(2^j, k)^* = \frac{1}{n_{a,j}} D_X(2^j) D_X(2^j)^*, \quad (3.6)$$

where $D_X(2^j) \in \mathcal{M}(p(n), n_{a,j}, \mathbb{R})$ is a random matrix with k-th column given by the vector $D_X(2^j, k) \in \mathbb{R}^{p(n)}$ (n.b.: the wavelet random matrix $W_X(2^j)$ involves $n_{a,j} = \frac{n}{a(n)2^j}$ terms, which depends on the scaling factor $a(n)$). Starting from the wavelet random matrix $W_X(2^j)$ as in (3.6), let

$$\check{W}(a(n)2^j) := P(n)a(n)^{H_n + (1/2)I} W_X(2^j)a(n)^{H_n + (1/2)I} P(n)^*, \quad (3.7)$$

where $H_n$ is given by (2.9). Then, as a consequence of the random operator self-similarity relation (C.95) (see Lemma C.16, (ii)), as well as of expressions (2.17) and (3.7),

$$\{\lambda_{\ell}(\mathbf{W}(a(n)2^j))\}_{\ell=1,\dots,p(n)} \stackrel{d}{=} \{\lambda_{\ell}(\check{\mathbf{W}}(a(n)2^j))\}_{\ell=1,\dots,p(n)}. \quad (3.8)$$

For convenience of notation, write $\check{W}_n = \check{W}(a(n)2^j)$. So, consider the e.c.d.f. of the rescaled log-eigenvalues coming from the right-hand side of (3.8), i.e.,

$$F_{\check{W}_n}(v) := \frac{1}{p} \sum_{\ell=1}^{p} \mathbb{1}_{\left\{\frac{\log \lambda_{\ell}(\check{W}_n)}{\log a(n)} \le v\right\}}, \quad v \in \mathbb{R}. \qquad (3.9)$$

Note that, as a consequence of Lemma B.2, (i), and of the fact that $P(n) \in GL(p(n), \mathbb{R})$ (assumption (A3)), (3.9) is well defined a.s. Let $v_0$ be a continuity point of $F_{2\mathcal{H}+1}$ and suppose, for the moment, that

$$F_{\check{W}_n}(v_0) \xrightarrow{P} F_{2\mathcal{H}+1}(v_0), \quad n \to \infty. \qquad (3.10)$$

Also, consider $F_{W_n}$ as defined in (3.4). We obtain

$$F_{W_n}(v_0) \stackrel{d}{=} F_{\check{W}_n}(v_0) \stackrel{d}{=} F_{2\mathcal{H}+1}(v_0), \quad n \to \infty, \qquad (3.11)$$

where equality in distribution in (3.11) is a consequence of (3.8). Since, in addition, $F_{2\mathcal{H}+1}(v_0)$ is a constant, then $F_{W_n}(v_0) \xrightarrow{P} F_{2\mathcal{H}+1}(v_0)$ as $n \to \infty$. This shows (3.5).

So, we need to establish the convergence (3.10). For this purpose, we first construct almost sure bounds for the eigenvalues $\lambda_{\ell}(\check{W}(a(n)2^j))$ as in (3.8). We do this by utilizing Weyl's inequalities for matrix products (B.2), and also by applying twice the cyclic property (B.5) of the eigenvalues of square matrix products. In fact, recall the notation (2.1) for singular values. Starting from expression (3.7), first we obtain

$$\lambda_{\ell}(a(n)^{H_n+(1/2)I} W_X(2^j)a(n)^{H_n+(1/2)I}) \cdot \sigma_1^2(\mathbf{P}(n)) \le \lambda_{\ell}(\check{W}(a(n)2^j)) \\ \le \lambda_{\ell}(a(n)^{H_n+(1/2)I} W_X(2^j)a(n)^{H_n+(1/2)I}) \cdot \sigma_p^2(\mathbf{P}(n)), \quad \ell = 1, \dots, p(n). \qquad (3.12)$$

Second, based on (3.12), we arrive at

$$\lambda_{\ell}(a(n)^{2H_n+I}) \cdot \lambda_1(W_X(2^j)) \cdot \sigma_1^2(P(n)) \le \lambda_{\ell}(\check{W}(a(n)2^j))$$
On the e.s.d. of large wavelet random matrices

$$ \leq \lambda_{\ell}(a(n)^{2\mathbb{H}_n+I}) \cdot \lambda_p(\mathbf{W}_X(2^j)) \cdot \sigma_p^2(\mathbf{P}(n)), \quad \ell = 1, \dots, p(n). \quad (3.13) $$

Recall that, by Lemma B.2, (i), $\log \lambda_{\ell}(\mathbf{W}_X(2^j))$ is well defined a.s. for $\ell = 1, \dots, p$. Therefore, from expression (3.13), taking logs while dividing by $\log a(n)$ yields

$$ \lambda_{\ell}(2\mathbb{H}_n + I) + G_L(n) \leq \frac{\log \lambda_{\ell}(\tilde{\mathbf{W}}(a(n)2^j))}{\log a(n)} \leq \lambda_{\ell}(2\mathbb{H}_n + I) + G_U(n), \quad (3.14) $$

for $\ell = 1, \dots, p$. In (3.14), we define

$$ G_L(n) := \frac{\log \lambda_1(\mathbf{W}_X(2^j))}{\log a(n)} + \frac{2 \log \sigma_1(\mathbf{P}(n))}{\log a(n)} \quad (3.15) $$

and

$$ G_U(n) := \frac{\log \lambda_p(\mathbf{W}_X(2^j))}{\log a(n)} + \frac{2 \log \sigma_p(\mathbf{P}(n))}{\log a(n)}. \quad (3.16) $$

Now observe that by making use of expression (2.8) (see assumption (A5)) and of Lemma B.2, (ii), it follows that

$$ G_L(n) \xrightarrow{P} 0 \quad \text{and} \quad G_U(n) \xrightarrow{P} 0, \quad (3.17) $$

as $n \to \infty$. Further note that

$$ \lambda_{\ell}(2\mathbb{H}_n + I) = 2\mathcal{H}_{(\ell)} + 1, \quad \ell = 1, \dots, p(n) \quad (3.18) $$

(cf. (2.2)), which stems from the fact that $\mathbb{H}_n$ as in (2.9) is a diagonal matrix.

So, we can now rewrite the bounds (3.14) in terms of the limits (3.17) and the equalities (3.18). That is, for every $n$ and each $\ell = 1, \dots, p(n)$, recast

$$ (2\mathcal{H}_{(\ell)} + 1) + o_{\mathbb{P}}(1)_L \leq \frac{\log \lambda_{\ell}(\tilde{\mathbf{W}}(a(n)2^j))}{\log a(n)} \leq (2\mathcal{H}_{(\ell)} + 1) + o_{\mathbb{P}}(1)_U \quad (3.19) $$

(n.b.: the subscripts L and U are used simply to make a distinction between the vanishing terms in the lower and upper bounds, respectively). Now fix $v \in \mathbb{R}$ and consider the e.c.d.f.s associated with the middle, right and left terms in (3.19), respectively. In other words, for $v \in \mathbb{R}$, consider $F_{\mathbf{W}_n}(v)$ as in (3.9), as well as

$$ F_{L_n}(v) := \frac{1}{p} \sum_{\ell=1}^{p} 1_{\{(2\mathcal{H}_{(\ell)}+1)+o_{\mathbb{P}}(1)_L \le v\}} \quad \text{and} \quad F_{U_n}(v) := \frac{1}{p} \sum_{\ell=1}^{p} 1_{\{(2\mathcal{H}_{(\ell)}+1)+o_{\mathbb{P}}(1)_U \le v\}}. \quad (3.20) $$

Then, for each $v \in \mathbb{R}$, the bounds (3.19) imply the (random) stochastic dominance relation

$$ F_{U_n}(v) \le F_{\tilde{\mathbf{W}}_n}(v) \le F_{L_n}(v) \quad \text{a.s.} \quad (3.21) $$

Now let $v_0$ be a point of continuity of $F_{2\mathbb{H}+1}$ as in (3.10). By setting $E_n = o_{\mathbb{P}}(1)_U$ and $E_n = o_{\mathbb{P}}(1)_L$ in Lemma C.10, (ii), we obtain, respectively,

$$ F_{U_n}(v_0) \xrightarrow{P} F_{2\mathbb{H}+1}(v_0) \quad \text{and} \quad F_{L_n}(v_0) \xrightarrow{P} F_{2\mathbb{H}+1}(v_0), \quad n \to \infty. \quad (3.22) $$

As a consequence, from the bounds (3.21) it follows that (3.10) holds, as claimed. □

**Remark 3.2.** It is informative to revisit the proof of Theorem 3.1 in light of some basic concepts from classical random matrix theory, pertaining to fixed-scale behavior.

So, consider the simplest framework of a single Hurst mode (i.e., $r = 1$ in (2.5)). Now recall the fact that the wavelet transform quasi-decorrelates the original measurements in the sense that the wavelet coefficients display fast-decaying correlations
(e.g., Flandrin [45], Bardet [18]). As a consequence of this property, over fixed scales the wavelet random matrix $n_{a,j}W(a(n)2^j)$ can be approximated by a classical Wishart matrix up to some constant $C_1 = C_1(p, a(n), j, n) > 0$. In other words, we can write $C_1 \cdot n_{a,j}W(a(n)2^j) \stackrel{d}{\approx} M^*M$, where the entries of the random matrix $M \in \mathcal{M}(n_{a,j}, p, \mathbb{R})$ are i.i.d. $\mathcal{N}(0, 1)$ random variables. In turn, for some constant $C_2 = C_2(p, n_{a,j}) > 0$, the joint density $f$ of the ordered eigenvalues of $M^*M$ satisfies

$$f(\lambda_1, \dots, \lambda_p) = C_2 \cdot 1_{\{0 \le \lambda_1 \le \dots \le \lambda_p\}} \cdot \prod_{\ell=1}^{p} \lambda_{\ell}^{(n_{a,j}-p-1)/2} \cdot \exp\left\{-\sum_{\ell=1}^{p} \lambda_{\ell}/2\right\} \cdot \underbrace{\prod_{i<\ell} (\lambda_{\ell} - \lambda_i)}_{\text{eigenvalue repulsion}} \quad (3.23)$$

(e.g., Anderson [10], p. 539). Based on the right-hand side of (3.23), we can interpret that, over fixed scales, there is a “force” based on which (wavelet) eigenvalues mutually repel (cf. Menon and Trogdon [72], chapter 8, or Tao [97], pp. 48–49). These conclusions naturally extend to log-eigenvalues.

The situation is markedly different over large scales. The proof of Theorem 3.1 (cf. relation (3.19)) allows us to interpret that two distinct “fractal forces” – i.e., both stemming from scaling effects – act on the rescaled wavelet log-eigenvalues. One “force” is repulsive, the other one, attractive (see Figure 2 for a schematic representation).

The repulsive force appears when $\pi(dH)$ displays more than one Hurst mode (i.e., $r > 1$ in (2.5)). It splits log-eigenvalues according to the different scaling laws trickling into the system over high dimensions. It mainly results from a combination of: (i) the scaling property of wavelet random matrices (see (3.8)); (ii) the variational characterization of eigenvalues, which appears implicitly in (3.12); (iii) the strong law of large numbers, which ensures that Hurst exponents are sampled asymptotically with the proportions given by the distribution $\pi$ (see (3.22) and its implicit connection with (B.13)).

On the other hand, the attractive force eventually overcomes fixed-scale repulsion (cf. (3.23)) and causes rescaled wavelet log-eigenvalues to collapse around each Hurst mode. This force stems from the vanishing small-scale features of the system under rescaling (see (3.17)).

As a consequence of the action of these two “forces”, in the three-way limit (2.7) (in particular, in the large scale limit), we are left with the distribution of Hurst (scaling) exponents in the rescaled logarithmic e.s.d. of wavelet random matrices.

## 3.2 Discrete time

Vis-à-vis Section 3.1, in this section we consider the more realistic situation in which measurements are made in discrete time (see (2.18)). In particular, recall that, in this context, the wavelet transform is given by (2.19).

The following theorem, Theorem 3.3, is the second main result of the paper. It is the discrete-time counterpart of Theorem 3.1. The fundamental conclusions in both theorems are analogous, as well as their proofs. However, the proof of Theorem 3.3 relies on developments for the continuous-time framework, since it involves controlling the difference between the spectra of wavelet random matrices for continuous- and discrete-time measurements (see Remark 3.4 for a short description).

Fix $j \in \mathbb{N} \cup \{0\}$ as in assumption (A4). As in Section 3.1, and, for notational simplicity, write

$$\tilde{\mathbf{W}}_n = \tilde{\mathbf{W}}(a(n)2^j). \quad (3.24)$$

To state the theorem, let

$$\lambda_1(\tilde{\mathbf{W}}_n) \le \dots \le \lambda_{p(n)}(\tilde{\mathbf{W}}_n)$$
be the ordered eigenvalues of the wavelet random matrix $\tilde{\mathbf{W}}_n$. Also, let

$$F_{\tilde{\mathbf{W}}_n}(dv) := \frac{1}{p} \sum_{\ell=1}^{p} \delta \left\{ \frac{\log \lambda_{\ell}(\tilde{\mathbf{W}}_n)}{\log a(n)} \leq \upsilon \right\} (dv) \quad (3.25)$$

be the e.s.d. of the rescaled log-eigenvalues of $\tilde{\mathbf{W}}_n$. In Theorem 3.3, stated and proven next, we show that $F_{\tilde{\mathbf{W}}_n}(dv)$ converges weakly, in probability, to the distribution $\pi(dH)$ of the Hurst exponents up to an affine transformation of the latter.

**Theorem 3.3.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Then, the e.s.d. $F_{\tilde{\mathbf{W}}_n}$ as in (3.25) converges weakly, in probability, to $F_{2\mathcal{H}+1}$. In other words, for any $g \in C_b(\mathbb{R})$ and for any $\varepsilon > 0$,

$$\lim_{n \to \infty} \mathbb{P}\left(\left|\int_{\mathbb{R}} g(v) F_{\tilde{\mathbf{W}}_n}(dv) - \int_{\mathbb{R}} g(v) F_{2\mathcal{H}+1}(dv)\right| > \varepsilon\right) = 0. \quad (3.26)$$

*Proof.* The proof follows a similar argument to the one for showing Theorem 3.1. Let

$$F_{\tilde{\mathbf{W}}_n}(v) := \frac{1}{p} \sum_{\ell=1}^{p} 1 \left\{ \frac{\log \lambda_{\ell}(\tilde{\mathbf{W}}_n)}{\log a(n)} \leq \upsilon \right\}, \quad v \in \mathbb{R}, \quad (3.27)$$

be the e.c.d.f. associated with $F_{\tilde{\mathbf{W}}_n}(dv)$ as in (3.25). As in the proof of Theorem 3.1, it suffices to show that, for any point of continuity $v_0$ of $F_{2\mathcal{H}+1}$,

$$F_{\tilde{\mathbf{W}}_n}(v_0) \xrightarrow{p} F_{2\mathcal{H}+1}(v_0), \quad n \to \infty. \quad (3.28)$$

So, recall relation (2.23). We define the auxiliary wavelet random matrix

$$\begin{align}
S_{\ge 0}(p(n), \mathbb{R}) &\ni \tilde{\mathbf{B}}_a(2^j) := a(n)^{-(H_n+(1/2)I)} \tilde{\mathbf{W}}_X(a(n)2^j)a(n)^{-(H_n+(1/2)I)} \nonumber \\
&= a(n)^{-(H_n+(1/2)I)} \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} \tilde{D}_X(a(n)2^j, k) \tilde{D}_X(a(n)2^j, k)^* a(n)^{-(H_n+(1/2)I)} \nonumber \\
&= a(n)^{-(H_n+(1/2)I)} \frac{1}{n_{a,j}} \tilde{D}_X(a(n)2^j) \tilde{D}_X(a(n)2^j)^* a(n)^{-(H_n+(1/2)I)}. \tag{3.29}
\end{align}$$

Now observe that, by (2.24), we can express the discrete-time wavelet random matrix (3.24) as

$$\tilde{\mathbf{W}}_n = \mathbf{P}(n)a(n)^{H_n+(1/2)I}\tilde{\mathbf{B}}_a(2^j)a(n)^{H_n+(1/2)I}\mathbf{P}(n)^*, \quad (3.30)$$

where $\tilde{\mathbf{B}}_a(2^j)$ is as in (3.29).

As in the proof of Theorem 3.1, by (B.5) and by Weyl's inequalities for matrix products (B.2), we obtain

$$\begin{align*}
\lambda_{\ell} (a(n)^{H_n+(1/2)I} \tilde{\mathbf{B}}_a(2^j) a(n)^{H_n+(1/2)I}) \cdot \sigma_1^2(\mathbf{P}(n)) &\leq \lambda_{\ell} (\tilde{\mathbf{W}}_n) \\
&\leq \lambda_{\ell} (a(n)^{H_n+(1/2)I} \tilde{\mathbf{B}}_a(2^j) a(n)^{H_n+(1/2)I}) \cdot \sigma_p^2(\mathbf{P}(n)), \quad \ell = 1, \dots, p(n).
\end{align*}$$

Then, again an application of (B.5) and (B.2) yields upper and lower bounds

$$\begin{align}
\lambda_{\ell} (a(n)^{2H_n+I}) \cdot \lambda_1 (\tilde{\mathbf{B}}_a(2^j)) \cdot \sigma_1^2 (\mathbf{P}(n)) &\leq \lambda_{\ell} (\tilde{\mathbf{W}}_n) \nonumber \\
&\leq \lambda_{\ell} (a(n)^{2H_n+I}) \cdot \lambda_p (\tilde{\mathbf{B}}_a(2^j)) \cdot \sigma_p^2 (\mathbf{P}(n)), \quad \ell = 1, \dots, p(n). \tag{3.31}
\end{align}$$

By Lemma B.5, (i), $\log \lambda_{\ell}(\tilde{\mathbf{B}}_a(2^j))$ is well defined a.s. for $\ell = 1, \dots, p(n)$. Therefore, from expression (3.31), taking logs while dividing by $\log a(n)$ yields

$$\lambda_{\ell}(2\mathcal{H}_n + I) + \tilde{G}_L(n) \leq \frac{\log \lambda_{\ell}(\tilde{\mathbf{W}}_n)}{\log a(n)} \leq \lambda_{\ell}(2\mathcal{H}_n + I) + \tilde{G}_U(n), \quad (3.32)$$
for $\ell = 1, \dots, p$. In (3.32), we define

$$ \tilde{G}_L(n) := \frac{\log \lambda_1(\tilde{\mathbf{B}}_a(2^j))}{\log a(n)} + \frac{2 \log \sigma_1(\mathbf{P}(n))}{\log a(n)}, \quad \tilde{G}_U(n) := \frac{\log \lambda_p(\tilde{\mathbf{B}}_a(2^j))}{\log a(n)} + \frac{2 \log \sigma_p(\mathbf{P}(n))}{\log a(n)}. $$

Moreover, by Lemma B.5, (ii), and (2.8) (see assumption (A5)), $\tilde{G}_L(n) \xrightarrow{P} 0$ and $\tilde{G}_U(n) \xrightarrow{P} 0$, as $n \to \infty$. The remainder of the proof is identical to that of Theorem 3.1, mutatis mutandis. This establishes (3.26). □

**Remark 3.4.** Proving Theorem 3.3 involves developments for the continuous-time framework since it requires controlling the deviation of the auxiliary wavelet random matrix (3.29) from its continuous-time analogue (B.31). In fact, this control in terms of the operator norm is provided in Lemma B.7. In turn, this lemma is used in the proof of Lemma B.5, (ii), which the proof of Theorem 3.3 relies on.

# 4 Conclusion, discussion and open problems

In this paper, we characterize the convergence of the (rescaled logarithmic) e.s.d. of wavelet random matrices. We assume a moderately high-dimensional framework where the sample size $n$, the dimension $p(n)$ and the scale $a(n)2^j$ go to infinity. We suppose the underlying measurement process is a random scrambling of a sample of size $n$ of a growing number $p(n)$ of fractional processes. Each of the latter processes is a fractional Brownian motion conditionally on a randomly chosen Hurst exponent. In our main results (Theorems 3.1 and 3.3), we show that the (rescaled logarithmic) e.s.d. of the wavelet random matrices converges weakly, in probability, to the distribution of Hurst exponents.

This work leads to a number of open questions, which we now briefly discuss.

With a view to universality claims, it is of great interest to better understand the extent to which the main results of this paper hold while some assumptions get lifted.

The moderately high-dimensional condition

$$ p(n) = o\left(\sqrt{\frac{n}{a(n)}}\right) \qquad (4.1) $$

(see (2.7)) is crucially used only once in this paper, namely, in the convergence of the truncated Haar random matrix $||U_{\text{trunc}}||_{\text{op}} \xrightarrow{P} 0$ (see (C.17) and (C.19)). This convergence undergirds Proposition C.2 (see the bound (C.23)), which in turn is used in the construction of the fixed-scale bounds (B.13) in Lemma B.2 and related auxiliary results. This leads to two open problems, namely, can the main conclusions of this paper be generalized from (4.1) to (i) the most general moderately high-dimensional framework $\lim_{n \to \infty} p(n) \cdot a(n)/n = 0$?; or (ii) to the high-dimensional framework $\lim_{n \to \infty} p(n) \cdot a(n)/n = c > 0$?

Likewise, it is in Lemma B.2 (and related auxiliary results) that we use the assumptions of a discrete distribution $\pi(dH)$ and of the independence among rH-fBms making up the latent process in (2.6). This leads to two other open problems, namely, can the main conclusions of this paper be generalized assuming (iii) a general Borel measure $\pi(dH)$ on (0,1)?; or (iv) a latent process $X$ with dependent rows? (here, preliminary simulation studies suggest a positive answer).

The condition $N_{\psi} \ge 2$ (see assumption (W1)) is used in Lemma C.13, (ii). In turn, this lemma is used in Lemma B.2, based on the bounds (B.7). However, the case where

$$ N_{\psi} = 1 \qquad (4.2) $$
is still of interest because of the following reason. Heuristically, under the choice (4.2), the wavelet random matrix $W(a(n)2^j)$ behaves, approximately, like a multiresolution sample covariance matrix (cf. Abry et al. [1], Remark 3.4). This is so because taking the wavelet transform based on a Haar basis, for example – for which (4.2) holds –, roughly corresponds to a differencing procedure (e.g., Mallat [61]). Hence, whether or not the conclusions of Theorems 3.1 and 3.3 remain true for the case (4.2) remains an open question.

Still related to the number of vanishing moments $N_\psi$, wavelet analysis in low dimensions further indicates that the assumption of a latent process $X$ as in (1.2) made up of (independent) stationary-increment processes is inessential. Namely, as long as $N_\psi$ is large enough, it is expected that, mutatis mutandis, the main conclusions of Theorems 3.1 and 3.3 hold for large (non-mixed-Gaussian) classes of stationary processes or stationary-increment processes of arbitrary order (cf. Moulines et al. [79] and Roueff and Taqqu [87]).

The research on the theory of wavelet random matrices is not only of mathematical interest. Rather, it naturally connects with a number of applications in the modeling of high-dimensional fractal systems, such as those emerging in neuroscience and fMRI imaging, network traffic, climate science, cointegration and high-dimensional time series. For those systems, statistically identifying Hurst multimodality from the available data demands the characterization of the fluctuations of the e.s.d. of wavelet random matrices, both in the bulk and around the edges (cf. Orejola et al. [81]). So far, though, this has only been done for measurements of the form of a high-dimensional signal-plus-noise model, with a fixed number of scaling exponents (Abry et al. [1, 2]).

# A Wavelet framework

## A.1 Assumptions on the wavelet multiresolution analysis

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

**Assumption (W1).** $\psi \in L^2(\mathbb{R})$ is a wavelet function, namely, it satisfies the relations

$$ \int_{\mathbb{R}} \psi^2(t) dt = 1, \quad \int_{\mathbb{R}} t^p \psi(t) dt = 0, \quad p = 0, 1, \dots, N_\psi - 1, \quad \int_{\mathbb{R}} t^{N_\psi} \psi(t) dt \neq 0, \quad (A.1) $$

for some integer (number of vanishing moments) $N_\psi \ge 2$.

**Assumption (W2).** The scaling and wavelet functions

$$ \phi \in L^1(\mathbb{R}) \text{ and } \psi \in L^1(\mathbb{R}) \text{ are compactly supported} \quad (A.2) $$

and $\hat{\phi}(0) = 1$.

**Assumption (W3).** There exist positive constants $A < \infty$ and $\alpha > 1$ such that

$$ \sup_{x \in \mathbb{R}} |\hat{\psi}(x)|(1 + |x|)^\alpha < A. \quad (A.3) $$

Conditions (A.1) and (A.2) imply that $\hat{\psi}(x)$ exists, is everywhere infinitely differentiable and its first $N_\psi - 1$ derivatives are zero at $x = 0$. Condition (A.3), in turn, implies that $\psi$ is continuous (see Mallat [61], Theorem 6.1) and, hence, bounded.

Note that assumptions (W1)-(W3) are closely related to the broad wavelet framework for the analysis of $\kappa$-th order ($\kappa \in \mathbb{N} \cup \{0\}$) stationary-increment stochastic processes laid out in Moulines et al. [77, 78, 79] and Roueff and Taqqu [87]. The Daubechies scaling and wavelet functions generally satisfy (W1)-(W3) (see Moulines et al. [79], p. 1927, or Mallat [61], p. 253). Usually, the parameter $\alpha$ increases to infinity as $N_\psi$ goes to
infinity (see Moulines et al. [79], p. 1927, or Cohen [31], Theorem 2.10.1). Also, under the orthogonality of the underlying wavelet and scaling function basis, (W1)-(W3) imply the so-called Strang-Fix condition (see Mallat [61], Theorem 7.4, and Moulines et al. [78], p. 159, condition (W-4)).

## A.2 Mallat's algorithm and discrete time measurements

Initially, suppose an infinite sequence of vector-valued measurements

$$ \{Y(k)\}_{k \in \mathbb{Z}} \subseteq \mathbb{R}^p \qquad (A.4) $$

is available. Then, we can apply Mallat's algorithm to extract the so-named approximation ($A(2^{j+1},)$) and detail ($\tilde{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

$$ \mathbb{R}^p \ni \tilde{Y}(t) := \sum_{k \in \mathbb{Z}} Y(k)\phi(t-k), \quad t \in \mathbb{R}. \qquad (A.5) $$

By the orthogonality of the shifted scaling functions {$\phi(\cdot - k)$}$_{k \in \mathbb{Z}}$,

$$ \mathbb{R}^p \ni A(2^0, k) = \int_{\mathbb{R}} \tilde{Y}(t)\phi(t-k)dt = Y(k), \quad k \in \mathbb{Z} \qquad (A.6) $$

(see Stoev et al. [96], proof of Lemma 6.1, or Moulines et al. [78], p. 160; cf. Abry and Flandrin [6], p. 33). In other words, the initial sequence, at octave $j = 0$, of approximation coefficients is given by the original sequence of random vectors. To obtain approximation and detail coefficients at coarser scales, we use Mallat's iterative procedure

$$ A(2^{j+1}, k) = \sum_{k' \in \mathbb{Z}} u_{k'-2k} A(2^j, k'), \quad \tilde{D}(2^{j+1}, k) = \sum_{k' \in \mathbb{Z}} v_{k'-2k} A(2^j, k'), \quad j \in \mathbb{N} \cup \{0\}, \quad k \in \mathbb{Z}, \qquad (A.7) $$

where the (scalar) filter sequences {$u_k := 2^{-1/2} \int_{\mathbb{R}} \phi(t/2)\phi(t-k)dt$}$_{k \in \mathbb{Z}}$, {$v_k := 2^{-1/2} \int_{\mathbb{R}} \psi(t/2)\phi(t-k)dt$}$_{k \in \mathbb{Z}}$ are called low- and high-pass MRA filters, respectively. Due to the assumed compactness of the supports of $\psi$ and of the associated scaling function $\phi$ (see condition (A.2)), only a finite number of filter terms is nonzero, which is convenient for computational purposes (Daubechies [35]). So, assume without loss of generality that $\text{supp}(\phi) = \text{supp}(\psi) = [0, T]$ (cf. Moulines et al. [78], p. 160). Then, the wavelet (detail) coefficients $\tilde{D}(2^j, k)$ of $Y$ can be expressed as (2.19), where the filter terms are defined by (2.20).

If we replace (A.4) with the realistic assumption that only a finite length series (2.18) is available, writing $\tilde{Y}^{(n)}(t) := \sum_{k=1}^n Y(k)\phi(t-k)$, we have $\tilde{Y}^{(n)}(t) = \tilde{Y}(t)$ for all $t \in (T, n+1)$ (cf. Moulines et al. [78]). Noting that $\tilde{D}(2^j, k) = \int_{\mathbb{R}} \tilde{Y}(t)2^{-j/2}\psi(2^{-j}t-k)dt$ and $\tilde{D}^{(n)}(2^j, k) = \int_{\mathbb{R}} \tilde{Y}^{(n)}(t)2^{-j/2}\psi(2^{-j}t-k)dt$, it follows that the finite-sample wavelet coefficients $\tilde{D}^{(n)}(2^j, k)$ of $\tilde{Y}^{(n)}(t)$ are equal to $\tilde{D}(2^j, k)$ whenever $\text{supp} \psi(2^{-j} \cdot -k) = (2^j k, 2^j(k+T)) \subseteq (T, n+1)$. In other words,

$$ \tilde{D}^{(n)}(2^j, k) = \tilde{D}(2^j, k), \quad \forall (j, k) \in \{(j, k) : 2^{-j}T \le k \le 2^{-j}(n+1) - T\}. \qquad (A.8) $$

Equivalently, such subset of finite-sample wavelet coefficients is not affected by the so-named border effect (cf. Craigmile et al. [33], Percival and Walden [83]). Moreover, by (A.8) the number of such coefficients at octave $j$ is given by $n_j = \lfloor 2^{-j}(n+1-T) - T \rfloor$. Hence, $n_j \sim 2^{-jn}$ for large $n$. Thus, for notational simplicity we suppose

$$ n_j = \frac{n}{2^j} \in \mathbb{N} \qquad (A.9) $$

holds exactly and only work with wavelet coefficients unaffected by the border effect.
## B Auxiliary results

In this section, we state and prove a number of lemmas and propositions used in the proofs of Theorems 3.1 and 3.3. In several arguments, we make use of the so-called Weyl's inequalities for both sums and products of Hermitian matrices. In particular, for $A, B \in \mathcal{M}(p, \mathbb{R})$, and for $\ell = 1, \dots, p$, we have the following Weyl's inequalities

$$ \sigma_{\ell}(A + B) \leq \sigma_{\ell}(A) + \sigma_{p}(B) \quad (B.1) $$

and

$$ \sigma_{\ell}(A) \cdot \sigma_{1}(B) \leq \sigma_{\ell}(AB) \leq \sigma_{\ell}(A) \cdot \sigma_{p}(B). \quad (B.2) $$

Moreover, if $A$ and $B$ are Hermitian, then

$$ \lambda_{\ell}(A) + \lambda_{1}(B) \leq \lambda_{\ell}(A + B) \leq \lambda_{\ell}(A) + \lambda_{p}(B), \quad \ell = 1, \dots, p, \quad (B.3) $$

and

$$ \max_{1 \leq \ell \leq p} |\lambda_{\ell}(A) - \lambda_{\ell}(B)| \leq \|A - B\|_{\text{op}}. \quad (B.4) $$

Proofs of the Weyl's inequalities (B.1)-(B.4) can be found in Merikoski and Kumar [73], Gray [47] and Vershynin [102]. In addition, we also make use of the classical cyclic property of eigenvalues. Namely, for $A, B \in \mathcal{M}(p, \mathbb{R})$,

$$ AB \text{ and } BA \text{ have the same eigenvalues.} \quad (B.5) $$

Throughout this and the next sections, further recall that we write $a = a(n)$ or $p = p(n)$ as in (2.10) whenever convenient.

### B.1 Continuous time

In this section, we provide some lemmas and propositions used in the proof of Theorem 3.1. First, we recap some notions associated with Toeplitz matrices.

**Definition B.1.** A Toeplitz matrix $T \in \mathcal{M}(n)$ is a matrix of the form

$$ T = \begin{bmatrix} \tau(0) & \tau(-1) & \tau(-2) & \dots & \tau(-(n-1)) \\ \tau(1) & \tau(0) & \tau(-1) & & \\ \tau(2) & \tau(1) & \tau(0) & & \\ \vdots & & \ddots & & \\ \tau(n-1) & & & \tau(0) \end{bmatrix}. $$

We say a sequence of Toeplitz matrices $T_n = \{\tau(\ell_1 - \ell_2)\}_{\ell_1, \ell_2=1, \dots, n}$, $n \in \mathbb{N}$, is generated by an integrable function $f: [-\pi, \pi] \to (0, \infty)$ when

$$ \tau(\ell) = \frac{1}{2\pi} \int_{-\pi}^{\pi} f(\lambda)e^{-i\ell\lambda}d\lambda, \quad \ell \in \mathbb{Z}. \quad (B.6) $$

Moreover, if $\{\tau(\ell)\}_{\ell \in \mathbb{Z}}$ is absolutely summable, then the sequence of Toeplitz matrices $\{T_n\}_{n \in \mathbb{N}}$ is said to be in the Wiener class.

Let $\{\Sigma_n\}_{n \in \mathbb{N}}$ be a sequence of symmetric Toeplitz matrices in the Wiener class, generated by a function $f$. By Gray [47], Lemma 4.1, for every $n \in \mathbb{N}$ and $\ell = 1, \dots, n$, the eigenvalues of $\Sigma_n$ satisfy

$$ m_f \leq \lambda_{\ell}(\Sigma_n) \leq M_f, \quad (B.7) $$

where

$$ m_f := \underset{x \in [-\pi, \pi]}{\operatorname{ess\,inf}} f(x) \quad \text{and} \quad M_f := \underset{x \in [-\pi, \pi]}{\operatorname{ess\,sup}} f(x). \quad (B.8) $$
We are now in a position to state and prove Lemma B.2. In this lemma, which is used in the proof of Theorem 3.1, we show that $W_X(2^j)$ is a full rank matrix a.s. Consequently, $\log \lambda_\ell(W_X(2^j))$ is well defined a.s. for $\ell = 1, ..., p(n)$. In addition, we establish the limiting behavior of

$$ \log \lambda_\ell(W_X(2^j)) / \log a(n). $$

To do so, we construct almost sure bounds on the eigenvalues of the fixed-scale wavelet random matrix $W_X(2^j)$ as in (3.6).

**Lemma B.2.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Let $W_X(2^j)$ be the auxiliary wavelet random matrix as in (3.6). Then, for $\ell = 1, ..., p(n)$,

(i) $\log \lambda_\ell(W_X(2^j))$ is well defined a.s.;

(ii) and

$$ \frac{\log \lambda_\ell(W_X(2^j))}{\log a(n)} \xrightarrow{\mathbb{P}} 0, \quad n \to \infty. \qquad (B.9) $$

*Proof.* First, we show (i). From expression (3.6), $W_X(2^j) \in S_{\ge 0}(n, \mathbb{R})$ a.s. Hence, $\lambda_1(W_X(2^j)) \ge 0$ a.s. Now recall that, by assumption (A4),

$$ p = p(n) < n_{a,j}. \qquad (B.10) $$

Since $\mathrm{rank}(W_X(2^j)) = \mathrm{rank}(D_X(2^j))$, then in light of (B.10) it suffices to show that $\mathrm{rank}(D_X(2^j)) = p(n)$ a.s. Indeed, by expression (C.92), $D_X(2^j)^* \stackrel{d}{=} (\Gamma_{\mathcal{H}_1} Z_1, \dots, \Gamma_{\mathcal{H}_p} Z_p)$, where $\mathbb{R}^{n_{a,j}} \supseteq Z_1, \dots, Z_p \stackrel{\mathrm{i.i.d.}}{\sim} \mathcal{N}(0, I)$ and $\Gamma_{\mathcal{H}_\ell}, \ell = 1, \dots, p$, are as in (C.83). Consequently, for any Hurst (matrix) exponent $\tilde{H}_n$ (see (B.37)), again by (B.10) and by Lemma C.9,

$$ \mathbb{P}(\mathrm{dim}(\mathrm{span}\{\Gamma_{\mathcal{H}_1} Z_1, \dots, \Gamma_{\mathcal{H}_p} Z_p\}) = p | \mathbb{H}_n = \tilde{\mathbb{H}}_n) = 1. $$

Since this holds for every $\tilde{\mathbb{H}}_n$, it follows that the rows of $D_X(2^j)$ are linearly independent a.s. Therefore, $\mathrm{rank}(D_X(2^j)) = p(n)$ a.s. This establishes (i).

We now turn to (ii). Recall the simple fact that, for any matrix $M \in \mathcal{M}(p, n, \mathbb{R})$ with $p = p(n) \le n$,

$$ \lambda_\ell(\mathbf{M}\mathbf{M}^*) = \lambda_{n+\ell-p}(\mathbf{M}^*\mathbf{M}), \quad \ell = 1, \dots, p. \qquad (B.11) $$

Therefore, by expression (3.6),

$$ \lambda_\ell(W_X(2^j)) = \lambda_\ell\left(\frac{1}{n_{a,j}} D_X(2^j) D_X(2^j)^*\right) = \lambda_{n_{a,j}-p+\ell}\left(\frac{1}{n_{a,j}} D_X(2^j)^* D_X(2^j)\right) \\ \stackrel{d}{=} \lambda_{n_{a,j}-p+\ell}\left(\frac{1}{n_{a,j}} \sum_{i=1}^r \Gamma_{\tilde{H}_i} Z_i Z_i^* \Gamma_{\tilde{H}_i}\right), \quad \ell = 1, \dots, p, \qquad (B.12) $$

where the equality in distribution follows from Lemma C.15, (iv), and the random matrices $Z_i$, $i = 1, \dots, r$, are defined in (C.89). Consider the random variables $p_i = p_i(n)$, $i = 1, \dots, r$, as in (C.88). By Kolmogorov's strong law of large numbers, $\lim_{n \to \infty} p_i(n)/p(n) = \mathbb{P}(\mathcal{H} = \tilde{\mathbb{H}}_i)$ a.s., for $i = 1, \dots, r$. Thus, from assumption (A4) and Lemma C.13, (ii), an application of Proposition B.4 to the right-hand side of (B.12) establishes the existence of strictly positive constants $m$ and $M$ such that, for $\ell = 1, \dots, p$,

$$ m + o_{\mathbb{P}}(1) \le \lambda_{n_{a,j}-p+\ell} \left( \frac{1}{n_{a,j}} \sum_{i=1}^{r} \Gamma_{\tilde{H}_i} Z_i Z_i^* \Gamma_{\tilde{H}_i} \right) \le M + o_{\mathbb{P}}(1). \qquad (B.13) $$

Note that the eigenvalues in (B.13) are strictly positive a.s. Therefore, by Lemma C.5, again for each $\ell = 1, \dots, p$,

$$ \frac{\log \lambda_{n_{a,j}-p+\ell} \left( \frac{1}{n_{a,j}} \sum_{i=1}^{r} \Gamma_{\tilde{H}_i} Z_i Z_i^* \Gamma_{\tilde{H}_i} \right)}{\log a(n)} \xrightarrow{\mathbb{P}} 0. \qquad (B.14) $$
As a consequence of expressions (B.12) and (B.14), the limit (B.9) follows. This completes the proof of (ii). □

Proposition B.4 is used in the proof of Lemma B.2. However, Proposition B.3 is needed in establishing Proposition B.4. We state and prove these two propositions next.

**Proposition B.3.** Let {$p_i$}$_{n \in \mathbb{N}}$ = {$p_i(n)$}$_{n \in \mathbb{N}}$ ⊆ $\mathbb{N}$ ∪ {0}, i = 1, 2, be two sequences of random variables such that

$$p_1 + p_2 = p \le n \quad \text{and} \quad p_1, p_2 \stackrel{\text{a.s.}}{=} o(\sqrt{n}) \qquad (\text{B.15})$$

(n.b.: in (B.15), $p = p(n)$ is deterministic). For $i = 1, 2$, if $p_i \ge 1$ let $\mathcal{Z}_i \in \mathcal{M}(n, p_i, \mathbb{R})$ be independent random matrices with independent, standard normal entries, and let $\mathcal{Z}_i = 0$ if $p_i = 0$. Also, for $i = 1, 2$, let $\Sigma_i \in \mathcal{M}(n)$ be symmetric Toeplitz matrices generated by functions $f_i$ such that, for $M_{f_i}$ and $m_{f_i}$ defined as in (B.8), $M_{f_i} \ge m_{f_i} > 0$. Let $\Gamma_i := \Sigma_i^{1/2}$, $i = 1, 2$. Then, almost surely,

$$\lambda_{n-p+1} \left( \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1 + \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \right) \\ \ge \max \left\{ \min \left\{ \frac{m_{f_2} m_{f_1}}{M_{f_2}}, m_{f_2} \right\}, \min \left\{ \frac{m_{f_1} m_{f_2}}{M_{f_1}}, m_{f_1} \right\} \right\} + o_P(1). \qquad (\text{B.16})$$

*Proof.* Observe that, by Weyl's inequalities (B.2), we obtain the two bounds

$$\lambda_{n-p+1} \left( \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1 + \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \right) \ge \lambda_1(\Sigma_2) \cdot \lambda_{n-p+1} \left( \Gamma_2^{-1} \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1 \Gamma_2^{-1} + \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \right) \quad (\text{B.17})$$

and

$$\lambda_{n-p+1} \left( \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1 + \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \right) \ge \lambda_1(\Sigma_1) \cdot \lambda_{n-p+1} \left( \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} + \Gamma_1^{-1} \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \Gamma_1^{-1} \right). \quad (\text{B.18})$$

Note that the right-hand side of inequalities (B.17) and (B.18) are of identical form. Moreover, the method by which we bound (B.18) will analogously apply to (B.17). So, to establish (B.16) we simply take the max of the two lower bounds. Therefore, without loss of generality we focus on (B.18). That is, we want to show that

$$\lambda_{n-p+1} \left( \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} + \Gamma_1^{-1} \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \Gamma_1^{-1} \right) \ge \min \left\{ \frac{m_{f_2}}{M_{f_1}}, 1 \right\} + o_P(1) \quad \text{a.s.} \qquad (\text{B.19})$$

To do so, consider the following decomposition of the two matrices in the sum

$$\frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} + \Gamma_1^{-1} \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \Gamma_1^{-1} = \mathbf{U}_1 \mathbf{D}_1 \mathbf{U}_1^* + \mathbf{U}_2 \mathbf{D}_2 \mathbf{U}_2^* \quad \text{a.s.} \qquad (\text{B.20})$$

In (B.20), $\mathbf{U}_1, \mathbf{U}_2 \in O(n)$ and $\mathbf{D}_1, \mathbf{D}_2$ are diagonal matrices satisfying

$$\lambda_{\ell}(\mathbf{D}_1) = \lambda_{\ell}\left(\frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n}\right) \quad \text{and} \quad \lambda_{\ell}(\mathbf{D}_2) = \lambda_{\ell}\left(\Gamma_1^{-1} \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \Gamma_1^{-1}\right)$$

for $\ell = 1, \dots, n$. Thus,

$$\lambda_{n-p+1} \left( \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} + \Gamma_1^{-1} \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 \Gamma_1^{-1} \right) = \lambda_{n-p+1} (\mathbf{D}_1 + \mathbf{U}_1^* \mathbf{U}_2 \mathbf{D}_2 \mathbf{U}_2^* \mathbf{U}_1).$$

Note that $\mathbf{U}_1$ and $\mathbf{U}_2$ are independent random matrices. Moreover, by Lemma C.4 we can assume that $\mathbf{U}_1 \sim \text{Haar}(O(n))$. Hence, by Lemma C.3, $\mathbf{U}_1^* \mathbf{U}_2 = \mathbf{U} \sim \text{Haar}(O(n))$. Thus,

$$\lambda_{n-p+1}(\mathbf{D}_1 + \mathbf{U}_1^* \mathbf{U}_2 \mathbf{D}_2 \mathbf{U}_2^* \mathbf{U}_1) \stackrel{d}{=} \lambda_{n-p+1}(\mathbf{D}_1 + \mathbf{U} \mathbf{D}_2 \mathbf{U}^*). \qquad (\text{B.21})$$
For now, suppose that, almost surely, the non-zero eigenvalues of $D_1$ and $D_2$ satisfy the boundedness relations

$$1 + o_{\mathbb{P}}(1) \le \lambda_{n-p_1+1}(\mathbf{D}_1) \le \lambda_n(\mathbf{D}_1) \le 1 + o_{\mathbb{P}}(1) \quad (B.22)$$

and

$$\frac{m_{f_2}}{M_{f_1}} + o_{\mathbb{P}}(1) \le \lambda_{n-p_2+1}(\mathbf{D}_2) \le \lambda_n(\mathbf{D}_2) \le \frac{M_{f_2}}{m_{f_1}} + o_{\mathbb{P}}(1). \quad (B.23)$$

Namely, $D_1$ and $D_2$ satisfy the boundedness relations (C.21) and (C.22). Then, by Proposition C.2,

$$\lambda_{n-p+1}(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*) \ge \min\left\{\frac{m_{f_2}}{M_{f_1}}, 1\right\} + o_{\mathbb{P}}(1) \quad \text{a.s.}$$

Moreover, under conditions (B.15), by Lemma C.6 and by the equality in distribution (B.21), we obtain

$$\lambda_{n-p+1}\left(\frac{\mathcal{Z}_1\mathcal{Z}_1^*}{n} + \Gamma_1^{-1}\Gamma_2\frac{\mathcal{Z}_2\mathcal{Z}_2^*}{n}\Gamma_1^{-1}\Gamma_2\right) \ge \min\left\{\frac{m_{f_2}}{M_{f_1}}, 1\right\} + o_{\mathbb{P}}(1).$$

This establishes (B.16).

So, we now need to establish relations (B.22) and (B.23). In fact, as a consequence of condition (B.15), by Proposition C.11,

$$\lambda_1\left(\frac{\mathcal{Z}_1^*\mathcal{Z}_1}{n}\right), \lambda_{p_1}\left(\frac{\mathcal{Z}_1^*\mathcal{Z}_1}{n}\right), \lambda_1\left(\frac{\mathcal{Z}_2^*\mathcal{Z}_2}{n}\right), \lambda_{p_2}\left(\frac{\mathcal{Z}_2^*\mathcal{Z}_2}{n}\right) \xrightarrow{\mathbb{P}} 1. \quad (B.24)$$

From relation (B.11), $\lambda_\ell(\mathcal{Z}_1^*\mathcal{Z}_1/n) = \lambda_{n-p_1+\ell}(\mathbf{D}_1)$ for $\ell = 1, \dots, p_1$. It follows from expression (B.24) that $D_1$ satisfies conditions (B.22), as claimed.

Similarly, from relation (B.11) and by Weyl's inequalities (B.2), for $\ell = 1, \dots, p_2$,

$$\lambda_1(\Sigma_2) \cdot \lambda_1(\Sigma_1^{-1}) \cdot \lambda_1\left(\frac{\mathcal{Z}_2^*\mathcal{Z}_2}{n}\right) = \lambda_1(\Sigma_2) \cdot \lambda_1(\Sigma_1^{-1}) \cdot \lambda_{n-p_2+1}\left(\frac{\mathcal{Z}_2\mathcal{Z}_2^*}{n}\right)$$

$$\le \lambda_{n-p_2+\ell}\left(\Gamma_1^{-1}\Gamma_2\frac{\mathcal{Z}_2\mathcal{Z}_2^*}{n}\Gamma_2^*\Gamma_1^{*-1}\right)$$

$$= \lambda_{n-p_2+\ell}(\mathbf{D}_2) \le \lambda_n(\Sigma_2) \cdot \lambda_n(\Sigma_1^{-1}) \cdot \lambda_{n-p_2+\ell}\left(\frac{\mathcal{Z}_2\mathcal{Z}_2^*}{n}\right) = \lambda_n(\Sigma_2) \cdot \lambda_n(\Sigma_1^{-1}) \cdot \lambda_{p_2}\left(\frac{\mathcal{Z}_2^*\mathcal{Z}_2}{n}\right).$$

Then, by the facts that $\lambda_1(\Sigma_1^{-1}) = \lambda_n(\Sigma_1)^{-1}$, $\lambda_n(\Sigma_1^{-1}) = \lambda_1(\Sigma_1)^{-1}$, and by expression (B.7), we arrive at the bounds

$$\frac{m_{f_2}}{M_{f_1}} \cdot \lambda_1\left(\frac{\mathcal{Z}_2^*\mathcal{Z}_2}{n}\right) \le \lambda_{n-p_2+\ell}(\mathbf{D}_2) \le \frac{M_{f_2}}{m_{f_1}} \cdot \lambda_{p_2}\left(\frac{\mathcal{Z}_2^*\mathcal{Z}_2}{n}\right), \quad (B.25)$$

for $\ell = 1, \dots, p_2$. Hence, from expression (B.24) and inequalities (B.25), it follows that $D_2$ also satisfies (B.23), as claimed. $\square$

The following proposition is used in the proof of Lemma B.2. It establishes upper and lower bounds on the eigenvalues appearing on the right-hand side of (B.12). Note that the lower bound (B.28) is not sharp. Instead, it is one possible bound we can obtain based on the max bound in (B.16).

**Proposition B.4.** Fix $r \in \mathbb{N}$. For $i = 1, \dots, r$, let $p_i = p_i(n) \in \mathbb{N} \cup \{0\}$, $i = 1, \dots, r$, be a nondecreasing sequence of random variables such that

$$\sum_{i=1}^{r} p_i = p = o(\sqrt{n}) \quad \text{and} \quad \frac{p_i}{p} \xrightarrow{\text{a.s.}} c_i \in (0, 1) \quad (B.26)$$
(n.b.: in (B.26), p = p(n) is deterministic). Moreover, for i = 1, ..., r, let $\mathcal{Z}_i$ be indepen-
dent random matrices as in Proposition B.3. In addition, for i = 1, ..., r, let $\Sigma_i \in \mathcal{M}(n)$
be symmetric Toeplitz matrices generated by a strictly positive function $f_i$ such that for
$M_{f_i}$ and $m_{f_i}$ defined as in (B.8), $M_{f_i} \geq m_{f_i} > 0$. Also, define $\Gamma_i = \Sigma_i^{1/2}$. Then,

$$
\lambda_n \left( \sum_{i=1}^{r} \Gamma_i \frac{\mathcal{Z}_i \mathcal{Z}_i^*}{n} \Gamma_i \right) \le \sum_{i=1}^{r} M_{f_i} + o_{\mathbb{P}}(1), \quad (B.27)
$$

whereas

$$
\lambda_{n-p+1} \left( \sum_{i=1}^{r} \Gamma_i \frac{\mathcal{Z}_i \mathcal{Z}_i^*}{n} \Gamma_i^* \right) \ge m_{f_r} \cdot \min \left\{ 1, \frac{m_{f_{r-1}}}{M_{f_r}}, \frac{m_{f_{r-1}} m_{f_{r-2}}}{M_{f_r} M_{f_{r-1}}}, \dots, \frac{\prod_{i=1}^{r-1} m_{f_i}}{\prod_{i=2}^{r} M_{f_i}} \right\} + o_{\mathbb{P}}(1). \tag{B.28}
$$

Proof. We first establish (B.27). By Weyl’s inequalities (B.2)-(B.3) as well as the bounds (B.7),

$$
\lambda_n \left( \sum_{i=1}^{r} \Gamma_i \frac{\mathcal{Z}_i \mathcal{Z}_i^*}{n} \Gamma_i \right) \le \sum_{i=1}^{r} \lambda_n \left( \Gamma_i \frac{\mathcal{Z}_i \mathcal{Z}_i^*}{n} \Gamma_i \right) \le \sum_{i=1}^{r} \lambda_n(\Sigma_i) \cdot \lambda_n \left( \frac{\mathcal{Z}_i \mathcal{Z}_i^*}{n} \right) \le \sum_{i=1}^{r} M_{f_i} \cdot \lambda_n \left( \frac{\mathcal{Z}_i \mathcal{Z}_i^*}{n} \right).
$$

Then, under condition (B.26), by Proposition C.11 we obtain (B.27).

In turn, to establish (B.28), we proceed in an inductive fashion. We only explicitly
derive a bound for the case $r = 3$, since the general case of $r \in \mathbb{N}$ can be established by
a simple adaptation of the argument.

Let $\mathbf{U}_2\mathbf{D}_2\mathbf{U}_2^* = \Gamma_3^{-1}\left(\Gamma_2\frac{\mathcal{Z}_2\mathcal{Z}_2^*}{n}\Gamma_2 + \Gamma_1\frac{\mathcal{Z}_1\mathcal{Z}_1^*}{n}\Gamma_1\right)\Gamma_3^{-1}$, where $\mathbf{U}_2 \in O(n)$ and $\lambda_\ell(\mathbf{D}_2) = \lambda_\ell\left(\Gamma_3^{-1}\left(\Gamma_2\frac{\mathcal{Z}_2\mathcal{Z}_2^*}{n}\Gamma_2 + \Gamma_1\frac{\mathcal{Z}_1\mathcal{Z}_1^*}{n}\Gamma_1\right)\Gamma_3^{-1}\right)$ for $\ell = 1, \dots, n$. Then,

$$
\begin{align}
\lambda_{n-(p_1+p_2)+1}(\mathbf{D}_2) &\ge \lambda_1(\Sigma_3^{-1}) \cdot \lambda_{n-(p_1+p_2)+1} \left(\Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 + \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1\right) \nonumber \\
&\ge \frac{m_{f_2}}{M_{f_3}} \min \left\{ \frac{m_{f_1}}{M_{f_2}}, 1 \right\} + o_{\mathbb{P}}(1) \quad \text{a.s.} \tag{B.29}
\end{align}
$$

where the first inequality is a consequence of Weyl's inequalities (B.2), and the second
follows from Proposition B.3, and expression (B.7). Moreover,

$$
\lambda_n(\mathbf{D}_2) \le \frac{M_{f_1}}{m_{f_3}} \lambda_n\left(\frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n}\right) + \frac{M_{f_2}}{m_{f_3}} \lambda_n\left(\frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n}\right) \le \frac{M_{f_1}}{m_{f_3}} + \frac{M_{f_2}}{m_{f_3}} + o_{\mathbb{P}}(1) \quad \text{a.s.,} \quad (B.30)
$$

where the first inequality follows from Weyl's inequalities (B.2)-(B.3) and bounds (B.7),
and the second inequality is a consequence of condition (B.26) and Proposition C.11.
Hence, $\mathbf{D}_2$ satisfies the bounds (C.21) and (C.22). In addition, by Lemma C.4, we can
recast $\frac{\mathcal{Z}_3 \mathcal{Z}_3^*}{n} \stackrel{d}{=} \mathbf{U}_1 \mathbf{D}_1 \mathbf{U}_1^*$ where $\mathbf{U}_1 \sim \text{Haar}(O(n))$ and $\lambda_\ell(\mathbf{D}_1) = \lambda_\ell(\frac{\mathcal{Z}_3 \mathcal{Z}_3^*}{n})$ for $\ell = 1, \dots, n$.
Then, as in the proof of Proposition B.3, by Proposition C.2 and Lemma C.6, it follows
that

$$
\lambda_{n-(p_1+p_2+p_3)+1} \left( \frac{\mathcal{Z}_3 \mathcal{Z}_3^*}{n} + \Gamma_3^{-1} \left( \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 + \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1 \right) \Gamma_3^{-1} \right) \\
\ge \min \left\{ 1, \frac{m_{f_2}}{M_{f_3}} \min \left\{ \frac{m_{f_1}}{M_{f_2}}, 1 \right\} \right\} + o_{\mathrm{P}}(1) = \min \left\{ 1, \frac{m_{f_2} m_{f_1}}{M_{f_3} M_{f_2}}, \frac{m_{f_2}}{M_{f_3}} \right\} + o_{\mathrm{P}}(1) \quad \text{a.s.}
$$

Therefore,

$$
\lambda_{n-(p_1+p_2+p_3)+1} \left( \Gamma_3 \frac{\mathcal{Z}_3 \mathcal{Z}_3^*}{n} \Gamma_3 + \Gamma_2 \frac{\mathcal{Z}_2 \mathcal{Z}_2^*}{n} \Gamma_2 + \Gamma_1 \frac{\mathcal{Z}_1 \mathcal{Z}_1^*}{n} \Gamma_1 \right) \\
\geq m_{f_3} \cdot \min \left\{ 1, \frac{m_{f_2} m_{f_1}}{M_{f_3} M_{f_2}}, \frac{m_{f_2}}{M_{f_3}} \right\} + o_{\mathbb{P}}(1).
$$

This shows (B.28).

<figure><img src="image_6.png" /></figure>
## B.2 Discrete time

The following lemma is the discrete-time analogue of Lemma B.2. The lemma is used in the proof of Theorem 3.3. In the lemma, we make use of the continuous-time analogue of the auxiliary wavelet random matrix (3.29). So, assuming continuous-time measurements, let

$$ 
\begin{align} 
S_{\ge 0}(p(n), \mathbb{R}) &\ni \hat{\mathbf{B}}_a(2^j) \nonumber \\ 
&:= a(n)^{-(\mathbb{H}_n+(1/2)I)} \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} D_X(a(n)2^j, k) D_X(a(n)2^j, k)^* a(n)^{-(\mathbb{H}_n+(1/2)I)} \nonumber \\ 
&\stackrel{d}{=} \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} D_X(2^j, k) D_X(2^j, k)^* = \mathbf{W}_X(2^j), \tag{B.31}
\end{align} 
$$

where the equality in distribution follows from Lemma C.16, (ii) (n.b.: in discrete time, there is no exact self-similarity-like relation such as (B.31)).

**Lemma B.5.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Let $\tilde{\mathbf{B}}_a(2^j)$ be a fixed-scale wavelet random matrix for discrete-time measurements as in (3.29). Then, for $\ell = 1, \dots, p(n)$,

(i) $\log \lambda_\ell(\tilde{\mathbf{B}}_a(2^j))$ is well defined a.s.;

(ii) and

$$ \frac{\log \lambda_{\ell}(\tilde{\mathbf{B}}_{a}(2^{j}))}{\log a(n)} \xrightarrow{P} 0, \quad n \to \infty. \tag{B.32} $$

*Proof.* The proof of (i) is analogous to the proof of Lemma B.2, (i). In particular, from expression (3.29) and under condition (B.10), it suffices to show that $\text{rank}(a(n)^{-\mathbb{H}-(1/2)I} \tilde{D}_X(a(n)2^j)) = p(n)$ a.s. Moreover, since $a(n)^{-\mathbb{H}-(1/2)I}$ has full rank, we only need to show that $\text{rank}(\tilde{D}_X(a(n)2^j)) = p(n)$ a.s. Indeed, by the independence of the $p$-rH-fBms that make up $\{X(t)\}_{t \in \mathcal{T}}$ (see assumption (A2) and expression (2.6)), and by Lemma C.18, conditionally on $\mathbb{H}_n = \tilde{\mathbb{H}}_n$, the rows of $\tilde{D}_X(a(n)2^j)$ are equal in distribution to independent Gaussian random vectors. That is, for $i = 1, \dots, p$, there exist a full rank matrix $\tilde{\Gamma}_{H_i}$ and an independent standard normal vector $\mathbf{Z}_\ell \in \mathbb{R}^{n_{a,j}}$ such that $\tilde{\Gamma}_{H_i} \mathbf{Z}_\ell \in \mathbb{R}^{n_{a,j}}$ is equal, in law, to the $i$-th row of $\tilde{D}_X(a(n)2^j)$. Therefore, by Lemma C.9,

$$ P(\dim(\text{span}\{\tilde{\Gamma}_{H_1} \mathbf{Z}_1, \dots, \tilde{\Gamma}_{H_p} \mathbf{Z}_p\}) = p | \mathbb{H}_n = \tilde{\mathbb{H}}_n) = 1. $$

Since this holds for every $\tilde{\mathbb{H}}_n$, it follows that the rows of $D_X(2^j)$ are linearly independent a.s. Therefore, $\text{rank}(D_X(2^j)) = p(n)$ a.s. This establishes (i).

We now show (ii). Observe that, by an application of Weyl's inequalities (B.3),

$$ 
\begin{align} 
\lambda_1(\tilde{\mathbf{B}}_a(2^j) - \hat{\mathbf{B}}_a(2^j)) + \lambda_1(\hat{\mathbf{B}}_a(2^j)) &\le \lambda_\ell(\tilde{\mathbf{B}}_a(2^j)) \nonumber \\ 
&\le \lambda_p(\tilde{\mathbf{B}}_a(2^j) - \hat{\mathbf{B}}_a(2^j)) + \lambda_p(\hat{\mathbf{B}}_a(2^j)), \quad \ell = 1, \dots, p(n), \tag{B.33}
\end{align} 
$$

with $\hat{\mathbf{B}}_a(2^j)$ as in (B.31). In addition, as a consequence of Lemma B.7 and Weyl's inequality (B.4), it follows that

$$ |\lambda_1(\tilde{\mathbf{B}}_a(2^j) - \hat{\mathbf{B}}_a(2^j))| \xrightarrow{P} 0, \quad |\lambda_p(\tilde{\mathbf{B}}_a(2^j) - \hat{\mathbf{B}}_a(2^j))| \xrightarrow{P} 0. \tag{B.34} $$

Moreover, $\hat{\mathbf{B}}_a(2^j) \stackrel{d}{=} \mathbf{W}_X(2^j)$ (see expression (B.31)) implies that $\lambda_\ell(\tilde{\mathbf{B}}_a(2^j)) \stackrel{d}{=} \lambda_\ell(\mathbf{W}_X(2^j))$ for $\ell = 1, \dots, p(n)$. Then, by expressions (B.12), (B.13), and by Lemma C.6, there exist strictly positive constants $m, M$ such that

$$ o_P(1) + m \le \lambda_\ell(\tilde{\mathbf{B}}_a(2^j)) \le M + o_P(1) \quad \text{a.s., } \ell = 1, \dots, p(n). \tag{B.35} $$
Consequently, making use of (B.33), (B.34) and (B.35), we obtain

$$o_P(1) + m \le \lambda_\ell(\tilde{\mathbf{B}}_a(2^j)) \le M + o_P(1) \quad \text{a.s., } \ell = 1, \dots, p(n).$$

Hence, by part (i) of this lemma and by Lemma C.5, (B.32) holds. This completes the proof of (ii). □

The next proposition is analogous to Proposition C.1 in Abry and Didier [4]. It establishes that, at a given scale $2^j$, the root-mean squared difference between the discrete- and continuous-time wavelet coefficients for rH-fBm is bounded, up to a constant factor, by $2^{j/2}$ times the square root of the dimension. The proposition is used in Lemma B.7.

**Proposition B.6.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. For the fixed $j \in \mathbb{N} \cup \{0\}$ and for any $k \in \mathbb{Z}$, let $\tilde{D}_X(2^j, k)$ be the discrete-time wavelet coefficient given by Mallat's algorithm (see (A.7)). Then, there is a constant $C(\psi, \phi) \ge 0$, not depending on $j$ or $k$, such that

$$(E\|D_X(2^j, k) - \tilde{D}_X(2^j, k)\|^2)^{1/2} \le C(\psi, \phi) 2^{j/2} p(n)^{1/2}. \quad (B.36)$$

*Proof.* Conditionally on

$$\mathbb{H}_n = \tilde{\mathbb{H}}_n, \quad (B.37)$$

$X$ is an ofBm. So, let $H_\ell$ denote $\ell$-th entry of the diagonal of the matrix $\tilde{\mathbb{H}}_n$. Then,

$$E[\|D_X(2^j, k) - \tilde{D}_X(2^j, k)\|^2 | \mathbb{H}_n = \tilde{\mathbb{H}}_n] = \sum_{\ell=1}^{p(n)} E[|d_{\mathcal{H}_\ell}(2^j, k) - \tilde{d}_{\mathcal{H}_\ell}(2^j, k)|^2 | \mathcal{H}_\ell = H_\ell], \quad (B.38)$$

where $\tilde{d}_{\mathcal{H}_\ell}(2^j, k)$ is the univariate discrete wavelet transformation corresponding to the $\ell$-th entry of the vector $\tilde{D}_X(2^j, k)$. Suppose for the moment that, for each $\tilde{H} \in \text{supp } \pi(dH)$,

$$E|d_{\tilde{H}}(2^j, k) - \tilde{d}_{\tilde{H}}(2^j, k)|^2 \le C(\psi, \phi)2^j, \quad (B.39)$$

where $C(\psi, \phi)$ is a constant only depending on $\psi$ and $\phi$. Then, (B.36) is an immediate consequence of expressions (B.38) and (B.39).

So, it remains to show that (B.39) holds. Indeed, by Taqqu [100], the univariate standard fBm admits the time-domain representation

$$\{B_H(t)\}_{t \in \mathbb{R}} \stackrel{\text{f.d.d.}}{=} \left\{ \int_{\mathbb{R}} f_H(t, u) B(du) \right\}_{t \in \mathbb{R}}, \quad (B.40)$$

where

$$f_H(t, u) = \begin{cases} \frac{\sqrt{\Gamma(2H+1) \sin(\pi H)}}{\Gamma(H+\frac{1}{2})} \{(t-u)_+^{H-1/2} - (-u)_+^{H-1/2}\}, & H \in (0, 1/2) \cup (1/2, 1) \\ 1_{\{0 \le u \le t\}} - 1_{\{t \le u < 0\}}, & H = 1/2 \end{cases}. \quad (B.41)$$

In (B.40), $B(du)$ is a real-valued, orthogonal-increment Gaussian random measure, and in (B.41), $x_+ = x$ if $x \ge 0$ and $x_+ = 0$ otherwise. In light of (B.40), of the compact support of $\psi$ (see assumption (W2)) as well as of relation (C.81), an adaptation of the proof of Theorem 11.4.1 in Samorodnitsky and Taqqu [88] (a type of stochastic Fubini's theorem) yields

$$d_H(2^j, k) \stackrel{d}{=} 2^{-j/2} \int_{\mathbb{R}} \left( \int_{\mathbb{R}} \psi(2^{-j}t - k) f_H(t, s) dt \right) B(ds). \quad (B.42)$$

Then, under assumptions (W1)-(W3), the proofs of Lemma 6.1, Lemma 6.2 and Proposition 2.4 in Stoev et al. [96] can be directly adapted for the time-domain kernel (B.41).
However, for $H = 1/2$, we note that in order to adapt the proofs of Lemma 6.1 and Proposition 2.4 in Stoev et al. [96], the following identities are required, i.e.,

$$ \int_{\mathbb{R}} |f_{H=1/2}(t_1, u) - f_{H=1/2}(t_2, u)|^2 du = |t_1 - t_2| $$

and

$$ f_{H=1/2}(i + 2^j k, u) = f_{H=1/2}(i/2^j, u/2^j - k) + 1_{\{0 \le u/2^j < k\}} - 1_{\{k < u/2^j \le 0\}} \quad (B.43) $$

(see Stoev et al. [96], p. 1895, first and second columns, respectively). This yields expression (B.39). □

The following lemma is analogous to Lemma C.2 in Abry and Didier [4]. It is used in the proof of Lemma B.5.

**Lemma B.7.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Let $\tilde{\mathbf{B}}_a(2^j)$, $\hat{\mathbf{B}}_a(2^j)$ be as in (3.29) and (B.31), respectively. Then,

$$ \mathbb{E} \| \tilde{\mathbf{B}}_a(2^j) - \hat{\mathbf{B}}_a(2^j) \|_{\text{op}} \to 0, \quad n \to \infty. \quad (B.44) $$

*Proof.* By Minkowsky's inequality, the left-hand side of (B.44) is bounded by

$$ \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} \mathbb{E} \| a(n)^{-(\mathbb{H}_n+(1/2)I)} \tilde{D}_X(a(n)2^j, k) \tilde{D}_X(a(n)2^j, k)^* a(n)^{-(\mathbb{H}_n+(1/2)I)} \\ - a(n)^{-(\mathbb{H}_n+(1/2)I)} D_X(a(n)2^j, k) D_X(a(n)2^j, k)^* a(n)^{-(\mathbb{H}_n+(1/2)I)} \|_{\text{op}}. \quad (B.45) $$

However, for $k = 1, \dots, n_{a,j}$, the deviation between the discrete- and continuous-time (non-normalized) wavelet variance terms can be recast in the form

$$ \tilde{D}_X(a(n)2^j, k) \tilde{D}_X(a(n)2^j, k)^* - D_X(a(n)2^j, k) D_X(a(n)2^j, k)^* \\ = [\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)][\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)]^* \\ + [\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)] D_X(a(n)2^j, k)^* \\ + D_X(a(n)2^j, k)[\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)]^*. $$

For notational simplicity, let

$$ \tilde{D}_X(a(n)2^j, k) := a(n)^{-(\mathbb{H}_n+(1/2)I)} D_X(a(n)2^j, k). \quad (B.46) $$

Then, from the fact that the operator norm is sub-multiplicative, it follows that each term under the summation sign in (B.45) can be bounded by

$$ \begin{aligned} & \mathbb{E} \left( \|a(n)^{-(\mathbb{H}_n+(1/2)I)} \left\{ [\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)][\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)]^* \right. \right. \\ & \qquad \left. + [\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)] D_X(a(n)2^j, k)^* \right. \\ & \qquad \left. + D_X(a(n)2^j, k)[\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)]^* \right\} a(n)^{-(\mathbb{H}_n+(1/2)I)} \|_{\text{op}} \right) \\ & \le \mathbb{E} \left( \|a(n)^{-(\mathbb{H}_n+(1/2)I)} \|_{\text{op}}^2 \cdot \| [\tilde{D}_X(a(n)2^j, k) \right. \\ & \qquad \left. - D_X(a(n)2^j, k)][\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)]^* \|_{\text{op}} \right) \\ & \qquad + \mathbb{E} \left( \|a(n)^{-(\mathbb{H}_n+(1/2)I)} \|_{\text{op}} \cdot \| [\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)][\tilde{D}_X(a(n)2^j, k)^* \|_{\text{op}} \right) \end{aligned} $$
On the e.s.d. of large wavelet random matrices

$$
\begin{align}
&+\mathrm{E}\left(\|\tilde{D}_X(a(n)2^j, k)[\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)]^*\|_{\mathrm{op}} \cdot \|a(n)^{-(H_n+(1/2)I)}\|_{\mathrm{op}}\right) \nonumber \\
&= \mathrm{E}\left(\|a(n)^{-(H_n+(1/2)I)}\|_{\mathrm{op}}^2 \cdot \|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|^2\right) \nonumber \\
&+2\mathrm{E}\left(\|a(n)^{-(H_n+(1/2)I)}\|_{\mathrm{op}} \cdot \|\tilde{D}_X(a(n)2^j, k)\| \cdot \|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|\right). \tag{B.47}
\end{align}
$$

In (B.47), the last equality follows from the fact that, for a rank one matrix, $\mathbf{x}_1\mathbf{x}_2^{\top}$, with $\mathbf{x}_1, \mathbf{x}_2 \in \mathbb{R}^{p(n)}$, $\|\mathbf{x}_1\mathbf{x}_2^{\top}\|_{\mathrm{op}} = \|\mathbf{x}_1\|\|\mathbf{x}_2\|$. We now develop upper bounds for each term in the equality (B.47). Since $a(n)^{-(H_n+(1/2)I)}$ is a diagonal matrix, it follows that, for $\varpi$ as in (2.5) (see assumption (A1)),

$$
\|a(n)^{-(H_n+(1/2)I)}\|_{\mathrm{op}} \leq a(n)^{-(\varpi+1/2)}. \qquad (\mathrm{B}.48)
$$

Hence, by Proposition B.6 and expression (B.48),

$$
\begin{align}
& \mathrm{E}\left(\|a(n)^{-(H_n+(1/2)I)}\|_{\mathrm{op}}^2 \cdot \|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|^2\right) \nonumber \\
& \le a(n)^{-2\varpi-1} \mathrm{E}\left(\|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|^2\right) \nonumber \\
& \le a(n)^{-2\varpi} C(\psi, \phi)^2 2^j p(n). \tag{B.49}
\end{align}
$$

In addition, again by (B.48),

$$
\begin{align}
& \mathrm{E}\left(\|a(n)^{-(H_n+(1/2)I)}\|_{\mathrm{op}} \cdot \|\tilde{D}_X(a(n)2^j, k)\| \cdot \|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|\right) \nonumber \\
& \le a(n)^{-(\varpi+1/2)} \mathrm{E}\left(\|\tilde{D}_X(a(n)2^j, k)\| \cdot \|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|\right). \tag{B.50}
\end{align}
$$

Then, by the Cauchy-Schwarz inequality,

$$
\begin{align*}
& \mathrm{E}\left(\|\tilde{D}_X(a(n)2^j, k)\| \cdot \|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|\right) \\
& \le \left(\mathrm{E}\|\tilde{D}_X(a(n)2^j, k)\|^2 \cdot \mathrm{E}\|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|^2\right)^{1/2} \\
& = \left(\mathrm{E}\|D_X(2^j, k)\|^2 \cdot \mathrm{E}\|\tilde{D}_X(a(n)2^j, k) - D_X(a(n)2^j, k)\|^2\right)^{1/2},
\end{align*}
$$

where the equality is a consequence of (B.46) and Lemma C.16, (i). Observe that, since
the $p(n)$ entries of $D_X(2^j, k)$ are i.i.d.,

$$
\mathrm{E}\|D_X(2^j, k)\|^2 = \sum_{\ell=1}^{p(n)} \mathrm{E}d_{\mathcal{H}}(2^j, k)_{\ell}^2 = p(n) \cdot \mathrm{E}d_{\mathcal{H}}(2^j, k) = p(n) \cdot \mathrm{E}d_{\mathcal{H}}(2^j, 0)^2,
$$

where the last equality follows from the conditional stationarity of $d_{\mathcal{H}}(2^j, k)$ (see Lemma
C.13, (i)). Thus, by Proposition B.6, the right-hand side of expression (B.50) can be
bounded by

$$
a(n)^{-\varpi} p(n) \left( \mathrm{E} d_{\mathcal{H}}(2^j, 0)_{1}^{2} \right)^{1/2} C(\psi, \phi) (2^j)^{1/2}. \quad (B.51)
$$

Therefore, applying the bounds (B.49) and (B.51) to the terms in the sum (B.47) yields

$$
\mathrm{E}\|\tilde{\mathbf{B}}_a(2^j) - \hat{\mathbf{B}}_a(2^j)\|_{\mathrm{op}} \le A \cdot a(n)^{-2\varpi} \cdot p(n) + B \cdot a(n)^{-\varpi} \cdot p(n), \quad (\mathrm{B}.52)
$$

where A and B are positive constants that do not depend on n. By condition (2.7) (see
assumption (A4)), (B.44) holds. $\square$
C Additional results

In this section, we establish auxiliary results used throughout Sections B.1 and B.2. In regard to the notation, throughout this section recall that $O(n)$ and Haar$(O(n))$ denote, respectively, the orthogonal group on $\mathbb{R}^n$ and the Haar probability measure on $O(n)$.

The following lemma is used to establish Propositions B.3 and B.4. The lemma shows that a subset of $p_1 = p_1(n)$ Euclidean vectors and of $p_2 = p_2(n)$ columns obtained from a Haar-distributed random matrix are approximately orthogonal as long as $p_1 + p_2$ is sufficiently small with respect to the overall dimension $n \to \infty$.

**Lemma C.1.** For any $n \in \mathbb{N}$ and $1 \le p = p(n) \le n$, suppose that $p_i = p_i(n, \omega) \in \mathbb{N}$, $i = 1, 2$, are two sequences of random variables such that

$$p_1 + p_2 = p = o(\sqrt{n}) \qquad (C.1)$$

(n.b.: in (C.1), $p = p(n)$ is deterministic). Let

$$\mathcal{E} = \operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{p_1}\} \subseteq \mathbb{R}^n, \qquad (C.2)$$

where $\mathbf{e}_\ell$ denotes the $\ell$-th Euclidean vector. Also, let

$$\mathcal{U} \equiv \mathcal{U}_n = \operatorname{span}\{\mathbf{u}_1, \dots, \mathbf{u}_{p_2}\} \subseteq \mathbb{R}^n, \qquad (C.3)$$

where $\mathbf{u}_\ell = \mathbf{u}_{\ell,n}(\omega)$ denotes the $\ell$-th column of the random matrix

$$\mathbf{U} \equiv \mathbf{U}_n \sim \operatorname{Haar}(O(n)). \qquad (C.4)$$

The following claims hold.

$$(i) \qquad \dim(\operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{p_1}, \mathbf{u}_1, \dots, \mathbf{u}_{p_2}\}) = p \quad \text{a.s.} \qquad (C.5)$$

(ii) For any (random) sequence of unit vectors $\{v_n\}_{n \in \mathbb{N}} \subseteq \operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{p_1}, \mathbf{u}_1, \dots, \mathbf{u}_{p_2}\}$, almost surely there exist sequences of vectors $v_{\mathcal{E},n} \in \mathcal{E}$, $v_{\mathcal{U},n} \in \mathcal{U}$ based on which we can write the direct sum

$$v_n = v_{\mathcal{E},n} + v_{\mathcal{U},n}, \qquad (C.6)$$

such that

$$\|v_{\mathcal{E},n}\|^2 + \|v_{\mathcal{U},n}\|^2 \ge 1 - o_P(1) \qquad (C.7)$$

and

$$\max\{\|v_{\mathcal{E},n}\|^2, \|v_{\mathcal{U},n}\|^2\} \le 1 + o_P(1). \qquad (C.8)$$

**Proof.** We start by establishing (i). First, we prove (C.5) conditionally on fixed values of $p_1 = p_1(n)$ and $p_2 = p_2(n)$. For simplicity, we use the notation $\tilde{p}_1 = \tilde{p}_1(n)$ and $\tilde{p}_2 = \tilde{p}_2(n)$ for these deterministic values. Note that choosing any (deterministic number) $\tilde{p}_2$ of columns of $\mathbf{U}$ is equal, in law, to the following procedure: (i) pick $\tilde{p}_2$ vectors $\mathbb{R}^n \supseteq \{\mathbf{Z}_i\}_{i=1,\dots,\tilde{p}_2}^{\text{i.i.d.}} \sim \mathcal{N}(0, I_n)$ (ii) next, orthonormalize these $\tilde{p}_2$ vectors by means of the Gram-Schmidt algorithm (cf. Meckes [66], p. 8). In particular, the resulting vectors $\tilde{\mathbf{u}}_1, \dots, \tilde{\mathbf{u}}_{\tilde{p}_2}$ satisfy

$$\tilde{\mathbf{u}}_1, \dots, \tilde{\mathbf{u}}_{\tilde{p}_2} \stackrel{d}{=} \mathbf{u}_1, \dots, \mathbf{u}_{\tilde{p}_2}. \qquad (C.9)$$

Let $\mathbf{e}_1, \dots, \mathbf{e}_{\tilde{p}_1}$ be the first $\tilde{p}_1$ Euclidean vectors (n.b.: $\tilde{p}_1$ is deterministic). Then,

$$\mathbb{P}(\dim(\operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{\tilde{p}_1}, \mathbf{Z}_1, \dots, \mathbf{Z}_{\tilde{p}_2}\}) = p) = \mathbb{P}\left(\bigcap_{j=1}^{\tilde{p}_1} \{\mathbf{e}_j \notin \operatorname{span}\{\mathbf{Z}_1, \dots, \mathbf{Z}_{\tilde{p}_2}\}\}\right)$$
On the e.s.d. of large wavelet random matrices

$$
\begin{align*}
&= \mathbb{P}\left(\bigcap_{j=1}^{\tilde{p}_1} \{\mathbf{e}_j \notin \operatorname{span}\{\tilde{\mathbf{u}}_1, \dots, \tilde{\mathbf{u}}_{\tilde{p}_2}\}\right) = \mathbb{P}\left(\bigcap_{j=1}^{\tilde{p}_1} \{\mathbf{e}_j \notin \operatorname{span}\{\mathbf{u}_1, \dots, \mathbf{u}_{\tilde{p}_2}\}\right) \\
&= \mathbb{P}(\operatorname{dim}(\operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{\tilde{p}_1}, \mathbf{u}_1, \dots, \mathbf{u}_{\tilde{p}_2}\}) = p),
\end{align*}
$$

where the first and fourth equalities follow from Lemma C.8, the second equality follows
from the Gram-Schmidt algorithm, and the third equality is a consequence of (C.9). Since
$\dim(\operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{\tilde{p}_1}\}) = \tilde{p}_1$ and by Lemma C.9, $\dim(\operatorname{span}\{\mathbf{Z}_1, \dots, \mathbf{Z}_{\tilde{p}_2}\}) = \tilde{p}_2$ a.s., we now
consider the case where $\tilde{p}_2 \neq 0$. In particular, for $\tilde{p}_1, \tilde{p}_2 \neq 0$, a simple adaptation of the
proof of Lemma C.9 yields $\mathbb{P}(\dim(\operatorname{span}\{\mathbf{e}_1, \dots, \mathbf{e}_{\tilde{p}_1}, \mathbf{Z}_1, \dots, \mathbf{Z}_{\tilde{p}_2}\}) = \tilde{p}_1 + \tilde{p}_2 = p) = 1$. Since
this is true for any fixed pair $\tilde{p}_1, \tilde{p}_2$, then (C.5) also holds unconditionally. Thus, (i) is
established.

We now turn to (ii). The almost sure existence of the decomposition (C.6) is a
consequence of (C.5). Now, for any $v_{\mathcal{E},n}$ and $v_{\mathcal{U},n}$ as in (C.6), observe that $1 = \|v_n\|^2 =$
$\|v_{\mathcal{E},n}\|^2 + \|v_{\mathcal{U},n}\|^2 + 2v_{\mathcal{E},n}^* v_{\mathcal{U},n}$. Hereinafter, we write $v_{\mathcal{E}} = v_{\mathcal{E},n}$ and $v_{\mathcal{U}} = v_{\mathcal{U},n}$ for
notational simplicity. From this, we obtain the almost sure upper and lower bounds

$$
\|v_{\mathcal{E}}\|^2 + \|v_{\mathcal{U}}\|^2 + 2|v_{\mathcal{E}}^* v_{\mathcal{U}}| \ge 1 \ge \|v_{\mathcal{E}}\|^2 + \|v_{\mathcal{U}}\|^2 - 2|v_{\mathcal{E}}^* v_{\mathcal{U}}|.
$$

Hence,

$$
\|v_{\mathcal{E}}\|^2 + \|v_{\mathcal{U}}\|^2 \geq 1 - 2|v_{\mathcal{E}}^* v_{\mathcal{U}}| \quad \text{and} \quad \max\{\|v_{\mathcal{E}}\|^2, \|v_{\mathcal{U}}\|^2\} \leq 1 + 2|v_{\mathcal{E}}^* v_{\mathcal{U}}|. \quad (C.10)
$$

Therefore, from (C.10), to establish (C.7) and (C.8), it suffices to show that

$$
v_{\mathcal{E}}^{*} v_{\mathcal{U}} = o_{\mathbb{P}}(1). \qquad (C.11)
$$

So, let

$$
\mathbf{U}_{\text{trunc}} \in \mathcal{M}(p_1, p_2, \mathbb{R}) \qquad (\text{C.12})
$$

be the upper $p_1 \times p_2$ corner of the random matrix $\mathbf{U}_n$ as in (C.4), and consider the basis expansions

$$
v_{\mathcal{E}} = \sum_{i=1}^{p_1} \alpha_i \mathbf{e}_i, \quad v_{\mathcal{U}} = \sum_{\ell=1}^{p_2} \beta_{\ell} \mathbf{u}_{\ell} \qquad (C.13)
$$

(n.b.: the coefficients in (C.13) are random and depend on n). Then,

$$
v_{\mathcal{E}}^{*} v_{\mathcal{U}} = \sum_{i=1}^{p_1} \sum_{\ell=1}^{p_2} \alpha_i \mathbf{e}_i^* \mathbf{u}_{\ell} \beta_{\ell} = [\alpha_1 \ \dots \ \alpha_{p_1}] \begin{bmatrix} \mathbf{e}_1^* \mathbf{u}_1 & \dots & \mathbf{e}_1^* \mathbf{u}_{p_2} \\ \vdots & & \vdots \\ \mathbf{e}_{p_1}^* \mathbf{u}_1 & \dots & \mathbf{e}_{p_1}^* \mathbf{u}_{p_2} \end{bmatrix} \begin{bmatrix} \beta_1 \\ \vdots \\ \beta_{p_2} \end{bmatrix} = \alpha^* \mathbf{U}_{\text{trunc}} \beta, \tag{C.14}
$$

where $\alpha \in \mathbb{R}^{p_1}$ and $\beta \in \mathbb{R}^{p_2}$ are vectors of coefficients. Note that, by (C.13),

$$
\|\alpha\| = \|v_{\mathcal{E}}\|, \quad \|\beta\| = \|v_{\mathcal{U}}\|. \tag{C.15}
$$

Now let

$$
\mathbf{U}_{\text{trunc},p} \in \mathcal{M}(p, \mathbb{R}) \qquad (\text{C.16})
$$

be the upper $p \times p$ corner ($p \le n$) of the random matrix $\mathbf{U}_n$ as in (C.4). In particular, since $\max\{p_1, p_2\} \le p$, the matrix $\mathbf{U}_{\text{trunc}} \in \mathcal{M}(p_1, p_2, \mathbb{R})$ is embedded in the matrix $\mathbf{U}_{\text{trunc},p} \in \mathcal{M}(p, \mathbb{R})$. Hence,

$$
\| \mathbf{U}_{\text{trunc}} \|_{\text{op}} \le \| \mathbf{U}_{\text{trunc},p} \|_{\text{op}}. \qquad (\text{C.17})
$$

Then, by relation (C.14) and the Cauchy-Schwarz inequality, as well as by relations (C.10) and (C.17),

$$
|v_{\mathcal{E}}^{*} v_{\mathcal{U}}| \le \|v_{\mathcal{E}}\| \cdot \|v_{\mathcal{U}}\| \cdot \|\mathbf{U}_{\text{trunc}}\|_{\text{op}}
$$
$$ \leq \|\mathbf{U}_{\text{trunc}}\|_{\text{op}} \cdot (1 + 2|v_{\varepsilon}^{*}v_{\mathcal{U}}|) \leq \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} \cdot (1 + 2|v_{\varepsilon}^{*}v_{\mathcal{U}}|). \quad (C.18) $$

Consequently, we can write $|v_{\varepsilon}^{*}v_{\mathcal{U}}|(1 - 2\|\mathbf{U}_{\text{trunc},p}\|_{\text{op}}) \leq \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}}$. However, by Jiang [53], Theorem 3, under $p = o(\sqrt{n})$ (see condition (C.1)),

$$ \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} = o_{\mathbb{P}}(1). \quad (C.19) $$

So, fix $0 < \delta < 1$ and let $\Omega_{n,\delta} = \{\omega : \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} < \delta/2\}$. Then, for any $\varepsilon > 0$, relations (C.18) and (C.19) imply that

$$ \begin{aligned} \mathbb{P}(|v_{\varepsilon}^{*}v_{\mathcal{U}}| > \varepsilon) &= \mathbb{P}(|v_{\varepsilon}^{*}v_{\mathcal{U}}| > \varepsilon \cap \Omega_{n,\delta}) + \mathbb{P}(|v_{\varepsilon}^{*}v_{\mathcal{U}}| > \varepsilon \cap \Omega_{n,\delta}^{c}) \\ &\leq \mathbb{P}(\|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} > \varepsilon(1 - \delta)) + o(1) = o(1). \end{aligned} $$

This shows (C.11). Hence, (ii) is established. □

The next proposition is used in the proofs of Propositions B.3 and B.4. Note that, in the statement and proof, we make use of expressions appearing in Lemma C.1. The proposition shows that, for random $p_1, p_2 \in \mathbb{N}$, a rank $p_1$ positive diagonal matrix plus a rank $p_2$ symmetric positive definite perturbation tends to have rank $p_1 + p_2$.

**Proposition C.2.** For any $n \in \mathbb{N}$ and $1 \leq p = p(n) \leq n$, suppose that $p_i = p_i(n, \omega) \in \mathbb{N}$, $i = 1, 2$, are two sequences of random variables satisfying (C.1). Let

$$ \mathbf{D}_1 \equiv \mathbf{D}_{1,n} = \operatorname{diag}(\zeta_{1,n}, \dots, \zeta_{p_1,n}, 0, \dots, 0), \quad \mathbf{D}_2 \equiv \mathbf{D}_{2,n} = \operatorname{diag}(\eta_{1,n}, \dots, \eta_{p_2,n}, 0, \dots, 0) \quad (C.20) $$

be two independent sequences of random matrices in $S_{\geq 0}(n, \mathbb{R})$ with

$$ 0 < \zeta_{1,n} \leq \dots \leq \zeta_{p_1,n} \quad \text{and} \quad 0 < \eta_{1,n} \leq \dots \leq \eta_{p_2,n} \quad \text{a.s.} $$

Further assume that there exist positive constants $\zeta_{\min}$, $\zeta_{\max}$, $\eta_{\min}$ and $\eta_{\max}$ such that

$$ \zeta_{1,n} \geq \zeta_{\min} + o_{\mathbb{P}}(1) \quad \text{and} \quad \eta_{1,n} \geq \eta_{\min} + o_{\mathbb{P}}(1) \quad (C.21) $$

and

$$ \zeta_{p_1,n} \leq \zeta_{\max} + o_{\mathbb{P}}(1) \quad \text{and} \quad \eta_{p_2,n} \leq \eta_{\max} + o_{\mathbb{P}}(1). \quad (C.22) $$

Now let $\mathbf{U} \equiv \mathbf{U}_n \sim \text{Haar}(O(n))$ be a sequence of random matrices which are independent of $\mathbf{D}_1 = \mathbf{D}_{1,n}$ and $\mathbf{D}_2 = \mathbf{D}_{2,n}$. Then,

$$ \lambda_{n-p+1}(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*) \geq \min\{\eta_{\min}, \zeta_{\min}\} + o_{\mathbb{P}}(1). \quad (C.23) $$

*Proof.* Starting from the left-hand side of (C.23), note that $\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^* \in S_{\geq 0}(n, \mathbb{R})$. Now let $\mathcal{V}_p \subseteq \mathbb{R}^n$ be any $p$-dimensional subspace. By the Courant-Fischer min-max theorem (e.g., Horn and Johnson [49], p. 236),

$$ \lambda_{n-p+1}(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*) \geq \inf_{v \in \mathcal{V}_p \cap \mathbb{S}^{n-1}} v^*(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*)v. \quad (C.24) $$

Let $\mathcal{E}$ and $\mathcal{U}$ be as in (C.2) and (C.3), respectively. Define $\mathcal{V} := \mathcal{E} \oplus \mathcal{U}$. By Lemma C.1, $\dim \mathcal{V} = p$ a.s. Also, by the same lemma, for any $v \in \mathcal{V}$ we can almost surely find $v_{\mathcal{E}} \in \mathcal{E}$ and $v_{\mathcal{U}} \in \mathcal{U}$ such that $v = v_{\mathcal{E}} + v_{\mathcal{U}}$. Thus,

$$ \begin{aligned} v^*(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*)v &= (v_{\mathcal{E}}^* + v_{\mathcal{U}}^*)(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*)(v_{\mathcal{E}} + v_{\mathcal{U}}) \\ &= v_{\mathcal{E}}^*\mathbf{D}_1v_{\mathcal{E}} + 2v_{\mathcal{E}}^*\mathbf{D}_1v_{\mathcal{U}} + v_{\mathcal{U}}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_{\mathcal{E}} + 2v_{\mathcal{E}}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_{\mathcal{U}} + v_{\mathcal{U}}^*\mathbf{D}_1v_{\mathcal{U}} + v_{\mathcal{U}}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_{\mathcal{U}}. \end{aligned} \quad (C.25) $$

In light of (C.20), since $\mathbf{D}_1, \mathbf{U}\mathbf{D}_2\mathbf{U}^* \in S_{\geq 0}(n, \mathbb{R})$, then

$$ v_{\mathcal{E}}^*\mathbf{D}_1v_{\mathcal{E}} \geq \|v_{\mathcal{E}}\|^2\zeta_{1,n}, \quad v_{\mathcal{U}}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_{\mathcal{U}} \geq \|v_{\mathcal{U}}\|^2\eta_{1,n}. \quad (C.26) $$
By (C.25) and (C.26), we obtain the bound

$$ \begin{aligned} v^*(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*)v &\ge \|v_\mathcal{E}\|^2\zeta_{1,n} + \|v_\mathcal{U}\|^2\eta_{1,n} + 2v_{\mathcal{U}}^*\mathbf{D}_1v_\mathcal{E} + 2v_\mathcal{E}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_\mathcal{U}. \\ &\ge \min\{\zeta_{1,n}, \eta_{1,n}\}(1 - o_\mathbb{P}(1)) + 2v_{\mathcal{U}}^*\mathbf{D}_1v_\mathcal{E} + 2v_\mathcal{E}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_\mathcal{U}, \end{aligned} \quad (C.27) $$

where the last inequality follows from relation (C.7) (see Lemma C.1). Suppose, for the moment, that

$$ v_{\mathcal{U}}^*\mathbf{D}_1v_\mathcal{E} = o_\mathbb{P}(1), \quad v_\mathcal{E}^*\mathbf{U}\mathbf{D}_2\mathbf{U}^*v_\mathcal{U} = o_\mathbb{P}(1). \quad (C.28) $$

Then, applying (C.28) to expression (C.24) with bound (C.27) yields

$$ \lambda_{n-p+1}(\mathbf{D}_1 + \mathbf{U}\mathbf{D}_2\mathbf{U}^*) \ge \min\{\zeta_{1,n}, \eta_{1,n}\}(1 - o_\mathbb{P}(1)) + o_\mathbb{P}(1). $$

Thus, by condition (C.21), we establish (C.23).

So, we need to show the two asymptotic relations in (C.28). We start off with the left one. Consider the basis expansions (C.13), in the proof of Lemma C.1. Let $\mathbf{D}_{1,\text{trunc}} := \text{diag}(\zeta_{1,n}, \dots, \zeta_{p_1,n}) \in S_{\ge 0}(p_1, \mathbb{R})$. Then,

$$ \begin{aligned} |v_{\mathcal{U}}^* \mathbf{D}_1 v_{\mathcal{E}}| &= \left| \sum_{\ell=1}^{p_2} \sum_{i=1}^{p_1} \beta_{\ell} \alpha_i \zeta_i e_i^* \mathbf{u}_{\ell} \right| = \left| \begin{bmatrix} \zeta_{1,n} \alpha_1 & \dots & \zeta_{p_1,n} \alpha_{p_1} \end{bmatrix} \begin{bmatrix} e_1^* \mathbf{u}_1 & \dots & e_1^* \mathbf{u}_{p_1} \\ \vdots & \ddots & \vdots \\ e_{p_1}^* \mathbf{u}_1 & \dots & e_{p_1}^* \mathbf{u}_{p_1} \end{bmatrix} \begin{bmatrix} \beta_1 \\ \vdots \\ \beta_{p_2} \end{bmatrix} \right| \\ &= |(\mathbf{D}_{1,\text{trunc}} \alpha)^* \mathbf{U}_{\text{trunc}} \beta| \le \|\mathbf{D}_{1,\text{trunc}} \alpha\| \cdot \|\mathbf{U}_{\text{trunc}} \beta\| \le \zeta_{p_1,n} \|\alpha\| \cdot \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} \|\beta\|, \end{aligned} \quad (C.29) $$

where $\mathbf{U}_{\text{trunc}}$ and $\mathbf{U}_{\text{trunc},p}$ are as in (C.12) and (C.16), respectively. In (C.29), the first and second inequalities follow, respectively, from an application of the Cauchy-Schwarz inequality and of the bound (C.17). Therefore, by relations (C.22), (C.15), (C.8) (see Lemma C.1) and (C.19), it follows that

$$ |v_{\mathcal{U}}^* \mathbf{D}_1 v_{\mathcal{E}}| \le \zeta_{p_1,n} \cdot \max\{\|v_{\mathcal{E}}\|^2, \|v_{\mathcal{U}}\|^2\} \cdot \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} = o_\mathbb{P}(1). \quad (C.30) $$

Hence, the left asymptotic relation in (C.28) is proved.

We now consider the right asymptotic relation in (C.28). Note that, for any $\ell = 1, \dots, p_2$, we can write $\mathbf{U}^*\mathbf{D}_2\mathbf{U}\mathbf{u}_\ell = \eta_{\ell,n}\mathbf{u}_\ell$. Then,

$$ \begin{aligned} |v_{\mathcal{E}}^* \mathbf{U}\mathbf{D}_2\mathbf{U}^* v_{\mathcal{U}}| &= \left| \sum_{\ell=1}^{p_2} \sum_{i=1}^{p_1} \beta_{\ell} \alpha_i (e_i^* \mathbf{U}) \mathbf{D}_2 (\mathbf{U}^* \mathbf{u}_{\ell}) \right| = \left| \sum_{\ell=1}^{p_2} \sum_{i=1}^{p_1} \beta_{\ell} \alpha_i \eta_{\ell,n} e_i^* \mathbf{u}_{\ell} \right| \\ &= \left| \begin{bmatrix} \eta_{1,n} \beta_1 & \dots & \eta_{p_2,n} \beta_{p_2} \end{bmatrix} \begin{bmatrix} e_1^* \mathbf{u}_1 & \dots & e_{p_1}^* \mathbf{u}_1 \\ \vdots & \ddots & \vdots \\ e_1^* \mathbf{u}_{p_2} & \dots & e_{p_1}^* \mathbf{u}_{p_2} \end{bmatrix} \begin{bmatrix} \alpha_1 \\ \vdots \\ \alpha_{p_1} \end{bmatrix} \right| \\ &= |(\mathbf{D}_{2,\text{trunc}} \beta)^* \mathbf{U}_{\text{trunc}} \alpha| \le \|\mathbf{D}_{2,\text{trunc}} \beta\| \cdot \|\mathbf{U}_{\text{trunc}}^* \alpha\| \le \eta_{p_2,n} \|\alpha\| \cdot \|\mathbf{U}_{\text{trunc},p}\|_{\text{op}} \cdot \|\beta\|, \end{aligned} $$

where $\mathbf{D}_{2,\text{trunc}} = \text{diag}(\eta_{1,n}, \dots, \eta_{p_2,n})$. Hence, again by relations (C.22), (C.15), (C.8) (see Lemma C.1) and (C.19), it follows that

$$ v_{\mathcal{E}}^* \mathbf{U}\mathbf{D}_2\mathbf{U}^* v_{\mathcal{U}} = o_{\mathbb{P}}(1). \quad (C.31) $$

This proves the right asymptotic relation in (C.28). Hence, the proposition is established.

□

The following lemma is used in the proof of Proposition B.3.

**Lemma C.3.** Let $\mathbf{U}_1, \mathbf{U}_2 \in O(n)$ be independent random matrices. Further assume that $\mathbf{U}_1 \sim \text{Haar}(O(n))$. Then, $\mathbf{U}_1^* \mathbf{U}_2 \sim \text{Haar}(O(n))$.
**Proof.** Let $\mathcal{B}$ be a set from the Borel $\sigma$-algebra defined on $O(n)$. Then, for any $\tilde{U}_2 \in O(n)$,

$$ P(U_1^* U_2 \in \mathcal{B} | U_2 = \tilde{U}_2) = P(U_1^* \tilde{U}_2 \in \mathcal{B} | U_2 = \tilde{U}_2) = P(U_1^* \in \mathcal{B} | U_2 = \tilde{U}_2) = P(U_1^* \in \mathcal{B}). \qquad (C.32) $$

The second equality in (C.32) follows from the fact that $U_1^* \tilde{U}_2 \sim \text{Haar}(O(n))$ for any fixed $\tilde{U}_2 \in O(n)$. In turn, the third equality follows from the independence of $U_1$ and $U_2$. Therefore,

$$ P(U_1^* U_2 \in \mathcal{B}) = E[P(U_1^* U_2 \in \mathcal{B} | U_2)] = P(U_1^* \in \mathcal{B}). $$

Therefore, $U_1^* U_2 \sim \text{Haar}(O(n))$, as claimed. $\square$

The next lemma is used in the proofs of Propositions B.3 and B.4.

**Lemma C.4.** For $n \in \mathbb{N}$, let $p \in \mathbb{N} \cup \{0\}$ be a random variable such that $p \le n$ a.s. If $p \ge 1$, let $\mathcal{Z} \in \mathcal{M}(n, p, \mathbb{R})$ be a random matrix with independent standard normal random entries, otherwise $\mathcal{Z} = \mathbf{0} \in \mathbb{R}^n$. Let $S_{\ge 0}(n, \mathbb{R}) \ni \mathbf{D} := \text{diag}(\lambda_1(\mathcal{Z}\mathcal{Z}^*), \dots, \lambda_n(\mathcal{Z}\mathcal{Z}^*))$. Then, we can write

$$ \mathcal{Z}\mathcal{Z}^* \stackrel{d}{=} \mathbf{U}\mathbf{D}\mathbf{U}^*, \qquad (C.33) $$

where $\mathbf{U} \sim \text{Haar}(O(n))$.

**Proof.** Note that, by the singular value decomposition of $\mathcal{Z}\mathcal{Z}^*$, there exists a matrix $\mathbf{U} \in O(n)$ such that $\mathcal{Z}\mathcal{Z}^* = \mathbf{U}\mathbf{D}\mathbf{U}^*$ (n.b: when $p=0$, $\mathcal{Z}\mathcal{Z}^* = \mathbf{U}\mathbf{D}\mathbf{U}^*$ holds for any $\mathbf{U} \in O(n)$). So, it remains to show that $\mathbf{U} \sim \text{Haar}(O(n))$.

Observe that, conditionally on $p = \tilde{p} \neq 0$, $\mathcal{Z} \in \mathcal{M}(n, \tilde{p}, \mathbb{R})$ is a random matrix with independent standard normal entries. In addition, $\mathcal{Z}$ is conditionally independent of $p$. Namely, $(\mathcal{Z}|p = \tilde{p}) \stackrel{d}{=} \tilde{\mathcal{Z}}$, where $\tilde{\mathcal{Z}} \in \mathcal{M}(n, \tilde{p}, \mathbb{R})$ is a matrix with independent, standard normal entries. Then, by the rotational invariance of the standard multivariate Gaussian, we can assume that $\mathbf{U} \sim \text{Haar}(O(n))$ conditionally on $p = \tilde{p} \neq 0$ (see Bai and Silverstein [13], pp. 334). In addition, for $\tilde{p} = 0$, we can simply let $\mathbf{U} \sim \text{Haar}(O(n))$. Now let $\mathcal{B}$ be a set from the Borel $\sigma$-algebra defined on $O(n)$. Then, for any $\mathbf{U} \sim \text{Haar}(O(n))$,

$$ P(\mathbf{U} \in \mathcal{B}) = \sum_{\tilde{p}=0}^{n} P(\mathbf{U} \in \mathcal{B} | p = \tilde{p}) \cdot P(p = \tilde{p}) = \sum_{\tilde{p}=0}^{n} P(\mathbf{U} \in \mathcal{B}) \cdot P(p = \tilde{p}) = P(\mathbf{U} \in \mathcal{B}). $$

This shows (C.33). $\square$

The next lemma contains generic facts used in the proof of Lemma B.2.

**Lemma C.5.** Consider three sequences of real-valued random variables $\{X_n\}_{n \in \mathbb{N}}$, $\{L_n\}_{n \in \mathbb{N}}$, and $\{U_n\}_{n \in \mathbb{N}}$ defined on a given probability space. Moreover, assume that for every $n \in \mathbb{N}$, $X_n > 0$ a.s. and

$$ m + L_n \le X_n \le M + U_n \quad \text{a.s.} \qquad (C.34) $$

for constants $M \ge m > 0$. Further assume that

$$ L_n, U_n \stackrel{\mathbb{P}}{\to} 0, \quad n \to \infty. \qquad (C.35) $$

Then,

$$ \frac{\log X_n}{b_n} \stackrel{\mathbb{P}}{\to} 0, \quad n \to \infty, \qquad (C.36) $$

for any sequence $\{b_n\}_{n \in \mathbb{N}}$ such that $b_n \to \infty$ as $n \to \infty$.

**Proof.** For any $n \in \mathbb{N}$, let $\Omega_n = \{\omega : L_n > -m\}$. Also, fix $\varepsilon > 0$. Then,

$$ P\left(\frac{|\log X_n|}{b_n} > \varepsilon\right) \le P\left(\frac{\log X_n}{b_n} > \varepsilon\right) + P\left(\frac{\log X_n}{b_n} < -\varepsilon\right). $$
However,

$$P\left(\frac{\log X_n}{b_n} > \varepsilon\right) \le P\left(\frac{\log(M+U_n)}{b_n} > \varepsilon\right) \to 0, \quad n \to \infty,$$

where the limit is a consequence of condition (C.35). On the other hand,

$$P\left(\frac{\log X_n}{b_n} < -\varepsilon\right) = P\left(\left\{\frac{\log X_n}{b_n} < -\varepsilon\right\} \cap \Omega_n\right) + P\left(\left\{\frac{\log X_n}{b_n} < -\varepsilon\right\} \cap \Omega_n^c\right) \\ \le P\left(\left\{\frac{\log(m+L_n)}{b_n} < -\varepsilon\right\} \cap \Omega_n\right) + o(1) = o(1),$$

where the last inequality and the last equality follow, again, from condition (C.35). This establishes (C.36). □

The next lemma is a generic fact used in the proofs of Proposition B.3 and Lemma B.5.

**Lemma C.6.** Consider three sequences of random variables $\{X_n\}_{n \in \mathbb{N}}$, $\{L_n\}_{n \in \mathbb{N}}$ and $\{U_n\}_{n \in \mathbb{N}}$ defined on a given probability space. Moreover, assume that for every $n \in \mathbb{N}$, $X_n > 0$ a.s. and that there exist constants $M \ge m > 0$ such that condition (C.34) holds. Further suppose that $L_n$, $U_n$ satisfy condition (C.35). Let $\{Y_n\}_{n \in \mathbb{N}}$ be another sequence of random variables, possibly defined on a different probability space, and assume that

$$Y_n \stackrel{d}{=} X_n \quad \text{for each } n \in \mathbb{N}. \qquad (C.37)$$

Then, there exist sequences of random variables $\{L'_n\}_{n \in \mathbb{N}}$, $\{U'_n\}_{n \in \mathbb{N}}$, such that

$$m + L'_n \le Y_n \le M + U'_n \quad \text{a.s.,}$$

where $L'_n = o_P(1)$, $U'_n = o_P(1)$.

Before we prove Lemma C.6, we state and show the following lemma, used in the proof of the former.

**Lemma C.7.** Let $\{X_n\}_{n \in \mathbb{N}}$ be a sequence of random variables defined on a given probability space. Let $m$ and $M$ be two real-valued constants such that $m < M$. Then, the following statements are equivalent.

(i) There exist sequences of random variables $\{L_n\}_{n \in \mathbb{N}}$ and $\{U_n\}_{n \in \mathbb{N}}$ satisfying (C.35) and such that

$$m + L_n \le X_n \le M + U_n \quad \text{a.s.} \qquad (C.38)$$

(ii) For every $\eta_L, \eta_U > 0$,

$$P(m - \eta_L < X_n < M + \eta_U) \to 1, \quad n \to \infty. \qquad (C.39)$$

*Proof.* First, we establish that (i) $\Rightarrow$ (ii). Suppose (i) holds, and fix $\eta_L, \eta_U > 0$. Let $\Omega_{L_n} = \{\omega : |L_n| < \eta_L\}$, $\Omega_{U_n} = \{\omega : |U_n| < \eta_U\}$. Then,

$$P(m - \eta_L < X_n < M + \eta_U) \ge P\left(\left\{m - \eta_L < X_n < M + \eta_U\right\} \cap \Omega_{L_n} \cap \Omega_{U_n}\right) \\ \ge P\left(\left\{m + L_n < X_n < M + U_n\right\} \cap \Omega_{L_n} \cap \Omega_{U_n}\right) = P(\Omega_{L_n} \cap \Omega_{U_n}) \to 1, \quad n \to \infty. \qquad (C.40)$$

In (C.40), the equality follows from condition (C.38). Hence, (ii) holds.

We now establish that (ii) $\Rightarrow$ (i). So, suppose (ii) holds. For $n \in \mathbb{N}$, let $L_n := -\max\{m - X_n, 0\}$, $U_n := \max\{X_n - M, 0\}$. Observe that

$$m - \max\{m - X_n, 0\} \le X_n \le M + \max\{X_n - M, 0\} \quad \text{a.s.,}$$
namely, relation (C.38) holds for this choice of sequences $L_n, U_n$. We now show that these sequences satisfy (C.35). Indeed, fix any $\varepsilon > 0$. Let $\Omega_n^M = \{\omega : X_n - M < \varepsilon/2\}$ and let $\Omega_n^m = \{\omega : m - X_n \le \varepsilon/2\}$. Then,

$$ \begin{aligned} \mathbb{P}(\max\{X_n-M,0\} > \varepsilon) &= \mathbb{P}(\{\max\{X_n-M,0\} > \varepsilon\} \cap \Omega_n^M) + \mathbb{P}(\{\max\{X_n-M,0\} > \varepsilon\} \cap (\Omega_n^M)^c) \\ &= \mathbb{P}(\{\max\{X_n-M,0\} > \varepsilon\} \cap (\Omega_n^M)^c) \le \mathbb{P}((\Omega_n^M)^c) = 1 - \mathbb{P}(X_n < M + \varepsilon/2) = o(1). \end{aligned} \quad (C.41) $$

In (C.41), the second equality follows from the fact that $\mathbb{P}(\{\max\{X_n - M, 0\} > \varepsilon\} \cap \Omega_n^M) = 0$, and the last equality is a consequence of (C.39). By an analogous reasoning, $\mathbb{P}(\max\{m - X_n, 0\} > \varepsilon) = o(1)$. This establishes relation (C.35), as claimed. Hence, (i) holds, as claimed. $\square$

We are now in a position to prove Lemma C.6.

*Proof of Lemma C.6.* By condition (C.37), for every $\eta_L, \eta_U > 0$,

$$ \mathbb{P}(m - \eta_L < X_n < M + \eta_U) = \mathbb{P}(m - \eta_L < Y_n < M + \eta_U). $$

Moreover, by condition (C.34) and Lemma C.7, it follows that $\mathbb{P}(m - \eta_L < Y_n < M + \eta_U) \to 1$ as $n \to \infty$. Consequently, another application of Lemma C.7 provides the existence of the sequences of random variables $\{L'_n\}_{n \in \mathbb{N}}$, $\{U'_n\}_{n \in \mathbb{N}}$. This completes the proof. $\square$

The next lemma is a simple linear algebra fact used in the proof of Lemma C.1. Its proof is omitted.

**Lemma C.8.** Fix $2 \le p \le n$. Consider two fixed integers $p_1, p_2 \ge 1$ such that $p_1 + p_2 = p \le n$. Let $\{v_1, \dots, v_{p_1}\} \subseteq \mathbb{R}^n$ and $\{u_1, \dots, u_{p_2}\} \subseteq \mathbb{R}^n$ each be a collection of linearly independent vectors. Then,

$$ \dim(\operatorname{span}\{v_1, \dots, v_{p_1}, u_1, \dots, u_{p_2}\}) = p \quad (C.42) $$

if and only if

$$ v_j \notin \operatorname{span}\{u_1, \dots, u_{p_2}\}, \quad j = 1, \dots, p_1. \quad (C.43) $$

In the next lemma, for the reader's convenience we state and prove the classical fact that a collection of independent multivariate Gaussian random vectors are linearly independent almost surely. The lemma is used in the proofs of Lemma B.2, Lemma B.5, and Lemma C.1.

**Lemma C.9.** For a fixed $p \le n$, consider the collection of random vectors $\Gamma_1 Z_1, \dots, \Gamma_p Z_p \in \mathbb{R}^n$, where $\mathbb{R}^n \supseteq \{Z_i\}_{i=1,\dots,p} \stackrel{\text{i.i.d}}{=} \mathcal{N}(0, I)$ and

$$ \Gamma_i \in S_{\ge 0}(n, \mathbb{R}) \text{ has full rank for each } i = 1, \dots, p. \quad (C.44) $$

Then,

$$ \dim(\operatorname{span}\{\Gamma_1 Z_1, \dots, \Gamma_p Z_p\}) = p \quad \text{a.s.} \quad (C.45) $$

*Proof.* First, we show that, for $n \ge p$,

$$ \mathbb{P}(\Gamma_p Z_p \in \operatorname{span}\{\Gamma_1 Z_1, \dots, \Gamma_{p-1} Z_{p-1}\}) = 0. \quad (C.46) $$

In fact,

$$ \begin{aligned} \mathbb{P}(\Gamma_p Z_p \in \operatorname{span}\{\Gamma_1 Z_1, \dots, \Gamma_{p-1} Z_{p-1}\}) &= \mathbb{P}(\{Z_p \in \operatorname{span}\{\Gamma_p^{-1} \Gamma_1 Z_1, \dots, \Gamma_p^{-1} \Gamma_{p-1} Z_{p-1}\}\}) \\ &= \int_{\mathbb{R}^n \times \dots \times \mathbb{R}^n} \mathbb{P}(Z_p \in \operatorname{span}\{\Gamma_p^{-1} \Gamma_1 Z_1, \dots, \Gamma_p^{-1} \Gamma_{p-1} Z_{p-1}\} | \{Z_i = z_i\}_{i=1,\dots,p-1}) \prod_{i=1,\dots,p-1} dF_{Z_i}(z_i) \end{aligned} $$
On the e.s.d. of large wavelet random matrices

$$= \int_{\mathbb{R}^n \times \dots \times \mathbb{R}^n} \mathbb{P}(\mathbf{Z}_p \in \operatorname{span}\{\Gamma_p^{-1} \Gamma_1 \mathbf{z}_1, \dots, \Gamma_p^{-1} \Gamma_{p-1} \mathbf{z}_{p-1}\}) \prod_{i=1, \dots, p-1} d\mathbf{Z}_i(\mathbf{z}_i). \quad (C.47)$$

In (C.47), the first equality follows from (C.44), and the second and third equalities follow from the law of total probability and the independence of $\mathbf{Z}_i \in \mathbb{R}^n$, $i = 1, \dots, p$. Since $n \ge p$ and $\dim(\operatorname{span}\{\Gamma_p^{-1} \Gamma_1 \mathbf{z}_1, \dots, \Gamma_p^{-1} \Gamma_{p-1} \mathbf{z}_{p-1}\}) \le p-1$, then the Lebesgue measure of $\operatorname{span}\{\Gamma_p^{-1} \Gamma_1 \mathbf{z}_1, \dots, \Gamma_p^{-1} \Gamma_{p-1} \mathbf{z}_{p-1}\})$ is zero. It follows that, in this case, $\mathbb{P}(\mathbf{Z}_p \in \operatorname{span}\{\Gamma_p^{-1} \Gamma_1 \mathbf{z}_1, \dots, \Gamma_p^{-1} \Gamma_{p-1} \mathbf{z}_{p-1}\}) = 0$. Therefore, (C.46) holds, as claimed.

To show (C.45), we proceed by induction on the collection of vectors. First, by (C.44), observe that

$$\mathbb{P}(\Gamma_1 \mathbf{Z}_1 = \mathbf{0}) = \mathbb{P}(\mathbf{Z}_1 = \mathbf{0}) = 0.$$

Equivalently, $\dim(\operatorname{span}\{\Gamma_1 \mathbf{Z}_1\}) = 1$ a.s., establishing the base case. Now, for $n \ge p$, suppose that

$$\dim(\operatorname{span}\{\Gamma_1 \mathbf{Z}_1, \dots, \Gamma_{p-1} \mathbf{Z}_{p-1}\}) = p-1 \text{ a.s.} \quad (C.48)$$

Then, as a consequence of (C.46), $\dim(\operatorname{span}\{\Gamma_1 \mathbf{Z}_1, \dots, \Gamma_p \mathbf{Z}_p\}) = p$ a.s. This completes the proof. □

In Lemma C.10, stated and shown next, we show the generic fact that a perturbation of a sample does not affect the asymptotic behavior of the associated e.c.d.f. Lemma C.10 is used in the proofs of Theorems 3.1 and 3.3.

**Lemma C.10.** Let $\{X_n\}_{n \in \mathbb{N}}$ be i.i.d random variables with common c.d.f. $F_X$.

(i) Let $\{Y_n\}_{n \in \mathbb{N}}$ be a sequence of random variables such that

$$Y_n \xrightarrow{\text{a.s.}} 0 \quad \text{as} \quad n \to \infty. \quad (C.49)$$

Then, for every $v$ that is a point of continuity of $F_X$,

$$\frac{1}{n} \sum_{\ell=1}^{n} 1_{\{X_{\ell}+Y_n \le v\}} \xrightarrow{\text{a.s.}} F_X(v), \quad n \to \infty. \quad (C.50)$$

(ii) Let $\{E_n\}_{n \in \mathbb{N}}$ be a sequence of random variables such that $E_n \xrightarrow{\mathbb{P}} 0$ as $n \to \infty$. Then, for every $v$ that is a point of continuity of $F_X$,

$$\frac{1}{n} \sum_{\ell=1}^{n} 1_{\{X_{\ell}+E_n \le v\}} \xrightarrow{\mathbb{P}} F_X(v), \quad n \to \infty. \quad (C.51)$$

*Proof.* Towards establishing (i), fix any $v$ that is a point of continuity of $F_X$. Note that, by Kolmogorov's strong law of large numbers, as $n \to \infty$,

$$\frac{1}{n} \sum_{\ell=1}^{n} 1_{\{X_{\ell} \le v\}} \xrightarrow{\text{a.s.}} \mathbb{P}(X_{\ell} \le v). \quad (C.52)$$

However, let $\varepsilon > 0$. Define $\Omega_0 = \{\omega : \text{the limits in (C.49) and (C.52) hold}\}$. In particular,

$$\mathbb{P}(\Omega_0) = 1. \quad (C.53)$$

Then, for $\omega \in \Omega_0$, as a result of relation (C.49) there exists $n_0(\omega)$ such that, for $\ell = 1, \dots, n$ and $n \ge n_0(\omega)$,

$$1_{\{X_{\ell}+\varepsilon \le v\}} \le 1_{\{X_{\ell}+|Y_n| \le v\}} \le 1_{\{X_{\ell}+Y_n \le v\}} \le 1_{\{X_{\ell}-|Y_n| \le v\}} \le 1_{\{X_{\ell}-\varepsilon \le v\}}. \quad (C.54)$$
After taking the sample average $\frac{1}{n} \sum_{\ell=1}^{n}$ across expression (C.54), relations (C.52) and (C.53) imply that, as $n \to \infty$,

$$F_X(v-\epsilon) \le \liminf_{n \to \infty} \frac{1}{n} \sum_{\ell=1}^{n} 1_{\{X_\ell+Y_n \le v\}} \le \limsup_{n \to \infty} \frac{1}{n} \sum_{\ell=1}^{n} 1_{\{X_\ell+Y_n \le v\}} \le F_X(v+\epsilon) \quad \text{a.s. (C.55)}$$

By considering a sequence $\epsilon_i \to 0$ as $i \to \infty$, the continuity of the c.d.f. $F_X$ at $v$ and relation (C.55) imply that (C.50) holds. This establishes (i).

We now turn to (ii). Fix any $v$ that is a point of continuity of $F_X$. Consider any subsequence $\frac{1}{n'} \sum_{\ell=1}^{n'} 1_{\{X_\ell+E_{n'} \le v\}}$. We can pass to any sub-subsequence such that $E_{n''} \xrightarrow{\text{a.s.}} 0$. Thus, by (i),

$$\frac{1}{n''} \sum_{\ell=1}^{n''} 1_{\{X_\ell+E_{n''} \le v\}} \xrightarrow{\text{a.s.}} F_X(v), \quad n'' \to \infty. \qquad (C.56)$$

That is, for any subsequence of $\frac{1}{n} \sum_{\ell=1}^{n} 1_{\{X_\ell+E_n \le v\}}$, there exists a sub-subsequence that almost surely converges to $F_X(v)$. This establishes (ii). $\square$

Consider a nondecreasing random sequence $p = p(n) \in \mathbb{N}$. Let $\mathcal{Z} \in \mathcal{M}(n, p, \mathbb{R})$ be a random matrix with independent standard normal entries, and consider the random matrix $\mathcal{Z}^*\mathcal{Z}/n$. In Proposition C.11, stated and proved next, we establish the almost sure convergence of $\|\mathcal{Z}^*\mathcal{Z}/n\|_{\text{op}}$ in the moderately high-dimensional regime. Proposition C.11 is used in the proofs of Propositions B.3 and B.4.

**Proposition C.11.** Let

$$p = p(n) \in \mathbb{N}, \quad p \le n, \qquad (C.57)$$

be a nondecreasing sequence of random variables such that

$$p/n \to 0 \quad \text{a.s.} \qquad (C.58)$$

For each $p(n)$, let $\mathcal{Z} \in \mathcal{M}(n, p(n), \mathbb{R})$ be a random matrix with i.i.d. standard normal entries. Then, as $n \to \infty$,

$$\left\| \frac{\mathcal{Z}^* \mathcal{Z}}{n} - I_p \right\|_{\text{op}} \xrightarrow{\text{a.s.}} 0. \qquad (C.59)$$

Moreover, as $n \to \infty$,

$$\lambda_1\left(\frac{\mathcal{Z}^*\mathcal{Z}}{n}\right), \lambda_p\left(\frac{\mathcal{Z}^*\mathcal{Z}}{n}\right), \lambda_n\left(\frac{\mathcal{Z}\mathcal{Z}^*}{n}\right) \xrightarrow{\text{a.s.}} 1. \qquad (C.60)$$

**Proof.** For convenience of notation, we let $\mathcal{Z}_p := \mathcal{Z}$. By condition (C.57), the eigenvalues on the left-hand side of (C.60) are well defined. Let $\Omega_0 = \{\{p(n)\}_{n \in \mathbb{N}} : \text{ the limit (C.58) holds}\}$ so that, by (C.58), $P(\Omega_0) = 1$. Write $\mathbf{p} := \{p(n)\}_{n \in \mathbb{N}} \in \mathbb{N}^n$, and let $\tilde{\mathbf{p}} = \{\tilde{p}(n)\}_{n \in \mathbb{N}}$ be a particular realization of the sequence $\mathbf{p}$. Also, let $P(d\tilde{\mathbf{p}})$ be the probability distribution of $\mathbf{p}$. Then, by conditioning on $\mathbf{p}$,

$$\begin{align*}
P\left(\left\|\frac{\mathcal{Z}_p^* \mathcal{Z}_p}{n} - I_p\right\|_{\text{op}} \to 0\right) &= \int_{\Omega_0} P\left(\left\|\frac{\mathcal{Z}_p^* \mathcal{Z}_p}{n} - I_p\right\|_{\text{op}} \to 0\right) \bigg|_{\mathbf{p} = \tilde{\mathbf{p}}} P(d\tilde{\mathbf{p}}) \\
&= \int_{\Omega_0} P\left(\left\|\frac{\mathcal{Z}_{\tilde{p}}^* \mathcal{Z}_{\tilde{p}}}{n} - I_{\tilde{p}(n)}\right\|_{\text{op}} \to 0\right) \bigg|_{\mathbf{p} = \tilde{\mathbf{p}}} P(d\tilde{\mathbf{p}}) \\
&= \int_{\Omega_0} P\left(\left\|\frac{\mathcal{Z}_{\tilde{p}}^* \mathcal{Z}_{\tilde{p}}}{n} - I_{\tilde{p}(n)}\right\|_{\text{op}} \to 0\right) \bigg|_{\mathbf{p} = \tilde{\mathbf{p}}} P(d\tilde{\mathbf{p}}) = 1. \tag{C.61}
\end{align*}$$
In (C.61), the second equality follows from the substitution principle with $Z_{\tilde{p}} \in \mathcal{M}(n, \tilde{p}(n))$, whose entries consist of i.i.d. standard normal random variables. In turn, the third equality results from the independence of $p$ and the entries of $Z_p$. That is, by construction, for each $p$ the sequence of matrices $Z_p$ has independent standard normal entries. Also, the last equality is a consequence of the fact that, for any deterministic sequence $\tilde{p}$ such that $\tilde{p}(n)/n \to 0$,

$$ \frac{1}{2} \sqrt{\frac{n}{\tilde{p}}} \left\| \frac{Z_{\tilde{p}}^* Z_{\tilde{p}}}{n} - I_{\tilde{p}} \right\|_{\text{op}} \xrightarrow{\text{a.s.}} 1 \quad \text{as} \quad n \to \infty $$

(see Chen and Pan [30], Theorem 1). This shows (C.59).

In regard to (C.60), note that, by Weyl's inequality (B.4),

$$ \sup_{1 \le l \le p} \left| \lambda_l \left( \frac{Z^* Z}{n} \right) - 1 \right| \le \left\| \frac{Z^* Z}{n} - I_p \right\|_{\text{op}}. $$

Hence, expression (C.60) is a direct consequence of (C.59) and (B.11). □

The subsequent lemmas are used to establish properties of the rH-fBm $X_H(t)$ (see (2.3)) and also of the univariate, continuous-time wavelet transform

$$ \mathbb{R}^{n_{a,j}} \ni \mathbf{d}_{\mathcal{H}}(2^j) := \{d_{\mathcal{H}}(2^j, k)\}_{k=1,\dots,n_{a,j}} := \left\{2^{-j/2} \int_{\mathbb{R}} \psi(2^{-j}t-k)X_{\mathcal{H}}(t)dt\right\}_{k=1,\dots,n_{a,j}} \quad (C.62) $$

(cf. (2.12)). We start off with the former. The following lemma is used in the proofs of Lemmas C.13 and C.14.

**Lemma C.12.** Let $X_H = \{X_H(t)\}_{t \in \mathbb{R}}$ be as in (2.3). Then,

(i) for any $s, t \in \mathbb{R}$ and $H \in (0, 1)$, the conditional covariance function of $X_H$ is given by

$$ \mathbb{E}[X_{\mathcal{H}}(s)X_{\mathcal{H}}(t)|\mathcal{H}=H] = \frac{|s|^{2H} + |t|^{2H} - |t-s|^{2H}}{2}. \quad (C.63) $$

In addition, the conditional covariance function (C.63) is continuous as a function of $(s, t)$;

(ii) $X_H$ has strictly stationary increments, i.e., for any $\tau \in \mathbb{R}$,

$$ \{X_{\mathcal{H}}(t+\tau) - X_{\mathcal{H}}(\tau)\}_{t \in \mathbb{R}} \stackrel{\text{f.d.d.}}{=} \{X_{\mathcal{H}}(t)\}_{t \in \mathbb{R}}. \quad (C.64) $$

Moreover, $X_H$ is conditionally self-similar, i.e., for any $c > 0$ and $H \in (0, 1)$,

$$ \{X_{\mathcal{H}}(ct)\}_{t \in \mathbb{R}} | \mathcal{H} = H \stackrel{\text{f.d.d.}}{=} \{c^{\mathcal{H}} X_{\mathcal{H}}(t)\}_{t \in \mathbb{R}} | \mathcal{H} = H. \quad (C.65) $$

*Proof.* By definition, conditionally on $\mathcal{H} = H$, $X_H$ is a fBm. Hence, (C.63) holds, as well as the continuity of the conditional covariance function. This establishes (i).

We now prove (ii). In fact, for any $\tau \in \mathbb{R}$, any $m \in \mathbb{N}$, any $t_1, \dots, t_m \in \mathbb{R}$ and any $B_1, \dots, B_m \in \mathcal{B}(\mathbb{R})$, we can write

$$ \begin{aligned} \mathbb{P} \left( \bigcap_{\ell=1}^{m} \{ X_{\mathcal{H}}(t_{\ell}+\tau) - X_{\mathcal{H}}(\tau) \in B_{\ell} \} \right) &= \int_{0}^{1} \mathbb{P} \left( \bigcap_{\ell=1}^{m} \{ X_{\mathcal{H}}(t_{\ell}+\tau) - X_{\mathcal{H}}(\tau) \in B_{\ell} \} \middle| \mathcal{H} = H \right) \pi(dH) \\ &= \int_{0}^{1} \mathbb{P} \left( \bigcap_{\ell=1}^{m} \{ X_{\mathcal{H}}(t_{\ell}) \in B_{\ell} \} \middle| \mathcal{H} = H \right) \pi(dH) = \mathbb{P} \left( \bigcap_{\ell=1}^{m} \{ X_{\mathcal{H}}(t_{\ell}) \in B_{\ell} \} \right). \end{aligned} $$

Equivalently, (C.64) holds. The conditional self-similarity property (C.65) promptly follows from the fact that, conditionally on $\mathcal{H} = H$, $X_H$ is a fBm. This establishes (ii). □
We now turn to the conditional properties of (C.62). The following lemma is used in the proofs of Lemma B.2 and Lemma C.15.

**Lemma C.13.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Fix $H \in (0, 1)$.

(i) Conditionally on $\mathcal{H} = H$, the wavelet transform $d_{\mathcal{H}}(2^j, k)$ is well defined in the mean squared sense for any $j \in \mathbb{N} \cup \{0\}$ and $k \in \mathbb{Z}$. In addition (still conditionally on $\mathcal{H} = H$), for any fixed $j \in \mathbb{N} \cup \{0\}$, the sequence $\{d_{\mathcal{H}}(2^j, k)\}_{k \in \mathbb{Z}}$ is strictly stationary, and for any $\kappa \in \mathbb{Z}$, it satisfies

$$ \mathbb{E}[d_{\mathcal{H}}(2^j, k + \kappa)d_{\mathcal{H}}(2^j, k)|\mathcal{H} = H] = \int_{-\pi}^{\pi} e^{ix\kappa}(2^j)^{1+2H}\alpha(H)^2 \sum_{\ell=-\infty}^{\infty} \frac{|\hat{\psi}(x+2\pi\ell)|^2}{|x+2\pi\ell|^{1+2H}} dx. \quad (C.66) $$

In (C.66), $\hat{\psi}$ denotes the Fourier transform of $\psi$ and

$$ \alpha(H)^2 := \frac{H\Gamma(2H)\sin(H\pi)}{\pi}. \quad (C.67) $$

In particular, conditionally on $\mathcal{H} = H$, $\{d_{\mathcal{H}}(2^j, k)\}_{k \in \mathbb{Z}}$ is wide-sense stationary with spectral density

$$ f_H(x) := (2^j)^{1+2H} \alpha(H)^2 \sum_{\ell=-\infty}^{\infty} \frac{|\hat{\psi}(x+2\pi\ell)|^2}{|x+2\pi\ell|^{1+2H}}. \quad (C.68) $$

(ii) The spectral density $f_H$ as in (C.68) is continuous and strictly positive on $[-\pi, \pi]$. In particular, there exist constants $M_H$ and $m_H$ such that $M_H \ge m_H > 0$ and

$$ m_H := \inf_{x \in [-\pi, \pi]} f_H(x) \quad \text{and} \quad M_H := \sup_{x \in [-\pi, \pi]} f_H(x). \quad (C.69) $$

(iii) For any $H \in (0, 1)$, let

$$ \Sigma_H = \Sigma_{H,n_{a,j}} = \left\{ \mathbb{E} \left[ d_{\mathcal{H}}(2^j, k_1) d_{\mathcal{H}}(2^j, k_2) \middle| \mathcal{H} = H \right] \right\}_{k_1, k_2 = 1, \dots, n_{a,j}} \in S_{\ge 0}(n_{a,j}, \mathbb{R}) \quad (C.70) $$

be the conditional covariance matrix of the wavelet coefficients in (C.62). Then, $\Sigma_H$ has full rank.

**Proof.** First, we show (i). As in Lemma C.12, conditionally on $\mathcal{H} = H$, $X_{\mathcal{H}}$ is a fBm with Hurst exponent $H$. Therefore, $d_{\mathcal{H}}(2^j, k)$ is well defined and the (fixed-scale) conditional process is also strictly stationary (see, for instance, Abry and Didier [4], Proposition 3.1, or Abry et al. [7], expression (39)). In addition, as a consequence of Flandrin [45], expression (34), we can compute the autocovariance of the wavelet coefficients as

$$ \mathbb{E}[d_{\mathcal{H}}(2^j, k + \kappa)d_{\mathcal{H}}(2^j, k)|\mathcal{H} = H] = \int_{\mathbb{R}} e^{ix\kappa}(2^j)^{1+2H}\alpha(H)^2 \frac{|\hat{\psi}(x)|^2}{|x|^{1+2H}} dx. \quad (C.71) $$

After a change of variables, we obtain (C.66). This proves (i).

We now establish (ii). For notational convenience, define the function

$$ g_H(x) := (2^j)^{1+2H} \alpha(H)^2 \frac{|\hat{\psi}(x)|^2}{|x|^{1+2H}} \mathbb{1}_{\{x \neq 0\}}. \quad (C.72) $$

By condition (A.2) (see (W2)),

$$ g_H(x) \text{ is infinitely differentiable at any } x \neq 0. \quad (C.73) $$
Then, the conditional spectral density of {$d_H(2^j, k)$}$_{k \in \mathbb{Z}}$ can be written as

$$f_H(x) = \sum_{\ell=-\infty}^{\infty} g_H(x + 2\pi\ell), \quad x \in [-\pi, \pi]. \qquad (C.74)$$

Next, we aim to establish the continuity of $f_H$. For small $x$, by Mallat [61], Theorem 7.4, we can write

$$\hat{\psi}(x) = \hat{\psi}^{N_{\psi}}(0)x^{N_{\psi}} + O(x^{N_{\psi}}), \quad x \to 0. \qquad (C.75)$$

In addition, observe that, for all $x \in \mathbb{R}$,

$$\frac{(2^j)^{1+2H} \alpha(H)^2}{|x|^{2H}} \le C \cdot \max\{1, |x|^{-2}\}. \qquad (C.76)$$

In (C.76), $C = (2^j)^3 \max_{H \in (0,1)} \alpha(H)^2$, where the maximum exists as a consequence of the continuity of the function $\alpha(H)$ (see (C.67)). Thus, from expressions (C.72), (C.75) and (C.76), we obtain the bound

$$0 \le g_H(x) \le O(x^{2N_{\psi}-3})$$

for small $x$. Since $N_{\psi} \ge 2$ (see (W1)), then the function $g_H(x)$ is continuous at $x = 0$. Hence, in view of (C.73), $g_H(x)$ is continuous for all $x \in \mathbb{R}$. Now, break up the sum in the spectral density as

$$f_H(x) = g_H(x) + \sum_{\ell \neq 0} g_H(x + 2\pi\ell).$$

Observe that, for $\ell \neq 0$ and $x \in [-\pi, \pi]$, by assumption (A.3) (see (W3)) and bound (C.76),

$$g_H(x + 2\pi\ell) \le C \cdot \max\{1, |x + 2\pi\ell|^{-2}\} \cdot \frac{|\hat{\psi}(x + 2\pi\ell)|^2}{|x + 2\pi\ell|} \le \frac{CA}{|x + 2\pi\ell|(1 + |x + 2\pi\ell|)^{2\alpha}}.$$

Thus, by the dominated convergence theorem, the function $x \to \sum_{\ell \neq 0} g_H(x + 2\pi\ell)$ is continuous for all $x \in [-\pi, \pi]$. Therefore, it follows that the spectral density $f_H(x)$ is continuous and bounded from above. In particular, there exists a scalar $M_H$ such that $M_H \ge f_H(x)$ for all $x \in [-\pi, \pi]$ as in (C.69).

We now show that $f_H(x) > 0$ for any $x \in \mathbb{R}$. Assume for the moment that, for any $\varepsilon > 0$, there exists $L = L(\varepsilon) \in \mathbb{N}$ such that

$$\sup_{x \in [-\pi, \pi]} \sum_{|\ell| > L} |\hat{\psi}(x + 2\pi\ell)|^2 < \varepsilon. \qquad (C.77)$$

Now note that, for $L$ as in (C.77), for any $\ell$ such that $|\ell| \le L$ and for $x \in [-\pi, \pi]$,

$$|x + 2\pi\ell|^{1+2H} \le |\pi + 2\pi\ell|||^{1+2H} \le |\pi + 2\pi L|^{1+2H} \Rightarrow \frac{1}{|\pi + 2\pi L|^{1+2H}} \le \frac{1}{|x + 2\pi\ell|^{1+2H}}. \qquad (C.78)$$

So, for $x \neq 0$, relation (C.97) (see Lemma C.17) and the bounds (C.78) and (C.77) imply that

$$\sum_{\ell=-\infty}^{\infty} \frac{|\hat{\psi}(x + 2\pi\ell)|^2}{|x + 2\pi\ell|^{1+2H}} \ge \sum_{|\ell| \le L} \frac{|\hat{\psi}(x + 2\pi\ell)|^2}{|\pi + 2\pi L|^{1+2H}} \ge \frac{1}{|\pi + 2\pi L|^3} \sum_{|\ell| \le L} |\hat{\psi}(x + 2\pi\ell)|^2 \ge \frac{1 - \varepsilon}{|\pi + 2\pi L|^3}.$$

On the other hand, at $x = 0$, the same argument applies, *mutatis mutandis*, starting from the expression

$$\sum_{\ell=-\infty}^{\infty} \frac{|\hat{\psi}(2\pi\ell)|^2}{|2\pi\ell|^{1+2H}} = \sum_{\ell \neq 0} \frac{|\hat{\psi}(2\pi\ell)|^2}{|2\pi\ell|^{1+2H}}.$$
Since $\alpha(H)^2 > 0$ for all $H \in (0, 1)$, $f_H(x)$ as in (C.68) is strictly positive for all $x \in [-\pi, \pi]$, as claimed. In particular, with $f_H(x)$ continuous and strictly positive on a compact set, there exists $m_H > 0$ as in (C.69).

So, we now need to show (C.77). By way of contradiction, suppose there exists $\varepsilon_0 > 0$ and a sequence $\{x_L\}_{L \in \mathbb{N}} \subseteq [-\pi, \pi]$ such that

$$ \sum_{|\ell| > L} |\hat{\psi}(x_L + 2\pi\ell)|^2 \geq \varepsilon_0. \qquad (C.79) $$

Since $[-\pi, \pi]$ is compact, without loss of generality we can assume that $x_L \to x_0 \in [-\pi, \pi]$ as $L \to \infty$. Fix $L' \in \mathbb{N}$. Then, for large enough $L$,

$$ \sum_{|\ell| \le L'} |\hat{\psi}(x_L + 2\pi\ell)|^2 \le \sum_{|\ell| \le L} |\hat{\psi}(x_L + 2\pi\ell)|^2 \le 1, $$

where the second inequality is a consequence of relation (C.97) (see Lemma C.17). Since $\hat{\psi}$ is continuous (see (A.2) in (W2)), we can sequentially take the double limit $\lim_{L' \to \infty} \liminf_{L \to \infty}$ to obtain

$$ 1 = \sum_{\ell=-\infty}^{\infty} |\hat{\psi}(x_0 + 2\pi\ell)|^2 \le \liminf_{L \to \infty} \sum_{|\ell| \le L} |\hat{\psi}(x_L + 2\pi\ell)|^2 \le 1, \qquad (C.80) $$

where the equality again follows from (C.97) (see Lemma C.17). Therefore,

$$ \varepsilon_0 + 1 \le \liminf_{L \to \infty} \sum_{|\ell| > L} |\hat{\psi}(x_L + 2\pi\ell)|^2 + \liminf_{L \to \infty} \sum_{|\ell| \le L} |\hat{\psi}(x_L + 2\pi\ell)|^2 = 1, $$

where the inequality stems from (C.79) and (C.80), and the equality is again a consequence of (C.97) (contradiction). This proves (C.77). Hence, (ii) is established.

We now turn to (iii). Observe that $\Sigma_H$ as in (C.70) is a symmetric Toeplitz matrix (see Definition B.1). Therefore, from part (ii) of this lemma (see also expressions (B.7) and (B.8)), it follows that $\Sigma_H$ is nonsingular for every $n_{a,j}$. This establishes (iii). $\square$

Next, we consider unconditional properties of the wavelet transform vectors. The following lemma proves that $d_H(2^j, k)$ and $D_X(2^j, k)$ are (unconditionally) well defined in the mean squared sense (cf. Lemma C.13). It is used in the proof of Theorem 3.1.

**Lemma C.14.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Then, for any $j \in \mathbb{N} \cup \{0\}$ and $k \in \mathbb{Z}$ we have the following.

(i) The wavelet transform $d_H(2^j, k)$ (see (C.62)) is well defined in the mean squared sense.

(ii) The multivariate wavelet transform $D_X(2^j, k)$ (see (2.12)) is well defined in the mean squared sense.

*Proof.* First, we show (i). Note that, by an application of the dominated convergence theorem, the mapping $\mathbb{R} \ni t \mapsto E|t|^\mathcal{H}$ is continuous. Thus, in view of Lemma C.12, (i), the covariance function $EX_\mathcal{H}(s)X_\mathcal{H}(t)$ is continuous. Therefore, by Cramér and Leadbetter [34], p. 86, it suffices to show that

$$ \int_{\mathbb{R}} \int_{\mathbb{R}} |EX_{\mathcal{H}}(s)X_{\mathcal{H}}(t)||\psi(s)||\psi(t)|dsdt < \infty. \qquad (C.81) $$

In fact, for any $s, t \in \mathbb{R}$, expression (C.63) (see Lemma C.12) implies that

$$ |EX_{\mathcal{H}}(s)X_{\mathcal{H}}(t)| = |E[E(X_{\mathcal{H}}(s)X_{\mathcal{H}}(t)|\mathcal{H})]| $$
$$ \leq \mathbb{E}(\sqrt{\mathbb{E}(X_{\mathcal{H}}(s)^2|\mathcal{H})} \sqrt{\mathbb{E}(X_{\mathcal{H}}(t)^2|\mathcal{H})}) = |\mathbb{E}|s|^\mathcal{H}|t|^\mathcal{H}| \leq \max\{1, |s|\} \cdot \max\{1, |t|\}. $$

Since $\psi$ is integrable with compact support (see assumptions (W1) and (W2)), (C.81) holds. This proves (i).

We now consider the multivariate setting. Recall that, by assumption (A2), $X$ is a $p$-variate process composed of $p$ independent rH-fBms. Thus, conditionally on $\mathbb{H}_n = \check{\mathbb{H}}_n$, $X$ is an ofBm with Hurst exponent (matrix) $\check{\mathbb{H}}_n$. Therefore, a natural adaptation of the proof of part (i) of this lemma establishes (ii). $\square$

The following lemma is used in the proof of Lemma B.2. The lemma establishes properties of the (unconditional) dependence structure of $d_{\mathcal{H}}(2^j)$ as in (C.62).

**Lemma C.15.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Fix $j \in \mathbb{N} \cup \{0\}$ and consider the random vector $d_{\mathcal{H}}(2^j)$ as in (C.62). For $H \in (0, 1)$ and $\Sigma_{\mathcal{H}}$ as in (C.70), let

$$ \Gamma_{\mathcal{H}} := \Sigma_{\mathcal{H}}^{1/2} \in S_{\ge 0}(n_{a,j}, \mathbb{R}). \qquad (C.82) $$

Now define the random matrix

$$ \Gamma_{\mathcal{H}} \in S_{\ge 0}(n_{a,j}, \mathbb{R}) \qquad (C.83) $$

by replacing $H$ with $\mathcal{H}$ in (C.82). The following claims hold.

(i) The random matrix

$$ \Gamma_{\mathcal{H}}^{-1} \text{ is well defined a.s.} \qquad (C.84) $$

and

$$ \mathbb{R}^{n_{a,j}} \ni \Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j) \stackrel{d}{=} \mathcal{N}(0, I). \qquad (C.85) $$

(ii) $\Gamma_{\mathcal{H}}$ and $\Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j)$ are conditionally independent. In other words, fix any (Borel) sets

$$ B_1 \in \mathcal{B}(S_{\ge 0}(n_{a,j}, \mathbb{R})), \quad B_2 \in \mathcal{B}(\mathbb{R}^{n_{a,j}}). \qquad (C.86) $$

Then,

$$ \mathbb{P}(\Gamma_{\mathcal{H}} \in B_1 \cap \Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j) \in B_2 | \mathcal{H}) = \mathbb{P}(\Gamma_{\mathcal{H}} \in B_1 | \mathcal{H}) \cdot \mathbb{P}(\Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j) \in B_2 | \mathcal{H}). \qquad (C.87) $$

(iii) $\Gamma_{\mathcal{H}}$ and $\Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j)$ are (unconditionally) independent.

(iv) Let $D_X(2^j)$ be the random matrix appearing in (2.13). Let

$$ p_i = p_i(n), \quad i = 1, \dots, n, \qquad (C.88) $$

be the random variable which counts the number of random variables $\mathcal{H}_\ell$, $\ell = 1, \dots, p$, equal to $\check{H}_i$. Then, we can write

$$ \frac{1}{n_{a,j}} D_X(2^j)^* D_X(2^j) \stackrel{d}{=} \frac{1}{n_{a,j}} \sum_{i=1}^r \Gamma_{\check{H}_i} Z_i Z_i^* \Gamma_{\check{H}_i}. \qquad (C.89) $$

In (C.89), $Z_i \in \mathcal{M}(n_{a,j}, p_i, \mathbb{R})$, $i = 1, \dots, r$, are independent random matrices containing only independent standard normal entries if $p_i \ge 1$, and $Z_i = 0$ if $p_i = 0$.

**Proof.** To show (i), first note that (C.84) is a consequence of Lemma C.13, (iii). Now let $B \in \mathcal{B}(\mathbb{R}^{n_{a,j}})$. Recall that, by Lemma C.13, the zero-mean random vector $d_{\mathcal{H}}(2^j) \in \mathbb{R}^{n_{a,j}}$ is Gaussian conditionally on $\mathcal{H} = H \in (0, 1)$. Then,

$$ \mathbb{P}(\Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j) \in B | \mathcal{H} = H) = \mathbb{P}(\Gamma_{\mathcal{H}}^{-1} d_{\mathcal{H}}(2^j) \in B | \mathcal{H} = H) = \mathbb{P}(\mathbf{Z} \in B | \mathcal{H} = H) = \mathbb{P}(\mathbf{Z} \in B), \qquad (C.90) $$
where $\mathbb{R}^{n_{a,j}} \ni \mathbf{Z} \sim \mathcal{N}(0, I)$. Then,

$$P(\Gamma_{\mathcal{H}}^{-1} \mathbf{d}_{\mathcal{H}}(2^j) \in B) = \mathrm{E}[P(\Gamma_{\mathcal{H}}^{-1} \mathbf{d}_{\mathcal{H}}(2^j) \in B | \mathcal{H})] = P(\mathbf{Z} \in B).$$

This establishes (C.85), and hence, (i).

We now turn to (ii). Fix $H \in (0, 1)$ and consider $B_1, B_2$ as in (C.86). Then,

$$P(\Gamma_{\mathcal{H}} \in B_1 \cap \Gamma_{\mathcal{H}}^{-1} \mathbf{d}_{\mathcal{H}}(2^j) \in B_2 | \mathcal{H} = H) = P(\Gamma_{\mathcal{H}} \in B_1 \cap \Gamma_{\mathcal{H}}^{-1} \mathbf{d}_{\mathcal{H}}(2^j) \in B_2 | \mathcal{H} = H). \quad (C.91)$$

We arrive at two possible cases. First, suppose $\Gamma_H \in B_1$. Then, the right-hand side of (C.91) can be rewritten as

$$P(\Gamma_{\mathcal{H}}^{-1} \mathbf{d}_{\mathcal{H}}(2^j) \in B_2 | \mathcal{H} = H) = P(\Gamma_H \in B_1 | \mathcal{H} = H) \cdot P(\Gamma_{\mathcal{H}}^{-1} \mathbf{d}_{\mathcal{H}}(2^j) \in B_2 | \mathcal{H} = H).$$

Alternatively, suppose $\Gamma_H \notin B_1$. Then, the right-hand side of (C.91) can be recast in the form

$$\begin{aligned} P(\Gamma_H \in B_1 \cap \Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2 | \mathcal{H} = H) &= 0 \\ &= P(\Gamma_H \in B_1 | \mathcal{H} = H) \cdot P(\Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2 | \mathcal{H} = H). \end{aligned}$$

Thus, in any case, the left-hand side of (C.91) equals

$$P(\Gamma_H \in B_1 | \mathcal{H} = H) \cdot P(\Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2 | \mathcal{H} = H).$$

This establishes (ii).

To show (iii), consider any sets $B_1, B_2$ as in (C.86). Then, by relations (C.87), (C.90) and (C.85),

$$\begin{aligned} P(\Gamma_H \in B_1 \cap \Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2) &= \mathrm{E}[P(\Gamma_H \in B_1 \cap \Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2 | \mathcal{H})] \\ &= \mathrm{E}[P(\Gamma_H \in B_1 | \mathcal{H}) \cdot P(\Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2 | \mathcal{H})] = \mathrm{E}[P(\Gamma_H \in B_1 | \mathcal{H})] \cdot P(\mathbf{Z} \in B_2) \\ &= P(\Gamma_H \in B_1) \cdot P(\mathbf{Z} \in B_2) = P(\Gamma_H \in B_1) \cdot P(\Gamma_H^{-1} \mathbf{d}_H(2^j) \in B_2). \end{aligned}$$

This establishes (iii).

We now turn to (iv). By the independence of the $p$ rH-fBms that make up $\{X(t)\}_{t \in \mathcal{T}}$ (see assumption (A2) and expression (2.6)) and by parts (i) and (iii) in this lemma, we can rewrite

$$D_X(2^j)^* \stackrel{d}{=} (\Gamma_{H_1} \mathbf{Z}_1, \dots, \Gamma_{H_p} \mathbf{Z}_p), \quad (C.92)$$

where $\mathbb{R}^{n_{a,j}} \supseteq \mathbf{Z}_1, \dots, \mathbf{Z}_p \stackrel{\text{i.d.}}{\sim} \mathcal{N}(0, I)$ and $\Gamma_{H_\ell}, \ell = 1, \dots, p$, are as in (C.83). Now, for each $i = 1, \dots, r$, let $p_i = p_i(n)$ be as in (C.88). Note that $\sum_{i=1}^r p_i(n) = p(n)$ a.s. Then,

$$\begin{aligned} S_{\ge 0}(n_{a,j}, \mathbb{R}) &\ni \frac{1}{n_{a,j}} D_X(2^j)^* D_X(2^j) \stackrel{d}{=} \frac{1}{n_{a,j}} \sum_{\ell=1}^p \Gamma_{H_\ell} \mathbf{Z}_\ell \mathbf{Z}_\ell^* \Gamma_{H_\ell} \\ &= \frac{1}{n_{a,j}} \sum_{i=1}^r \sum_{\ell=1}^p \Gamma_{\tilde{H}_i} \mathbf{Z}_{\ell,i} \mathbf{Z}_{\ell,i}^* \Gamma_{\tilde{H}_i} = \frac{1}{n_{a,j}} \sum_{i=1}^r \Gamma_{\tilde{H}_i} \mathbf{Z}_i \mathbf{Z}_i^* \Gamma_{\tilde{H}_i}, \end{aligned} \quad (C.93)$$

where $\mathbf{Z}_1, \dots, \mathbf{Z}_r$ are as in (C.89). In (C.93), for $i = 1, \dots, r$ and $\ell = 1, \dots, p$, $\mathbf{Z}_{\ell,i} \in \mathbb{R}^{n_{a,j}}$ are independent random vectors with independent standard normal entries if $p_i \ge 1$, and $\mathbf{Z}_{\ell,i} = 0$, otherwise. Also, $\mathbf{Z}_i$ is as in (C.89). This establishes (C.89), and hence, (iv). $\square$

The next lemma is used in the proof of Theorem 3.1. Accordingly, the lemma pertains to the continuous-time framework and establishes the scaling properties of $\{D_X(a(n)2^j, k)\}_{k \in \mathbb{Z}}$ and $\mathbf{W}_X(a(n)2^j)$. For the lemma, recall that $H_n$ is given by (2.9).
**Lemma C.16.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. Fix $j \in \mathbb{N} \cup \{0\}$.

(i) The wavelet transform $D_X(a(n)2^j, k)$ as in (2.13) satisfies a random operator self-similarity relation. That is,

$$ \{D_X(a(n)2^j, k)\}_{k \in \mathbb{Z}} \stackrel{\text{f.d.d.}}{=} \{a(n)^{\mathbb{H}_n+1/2 \cdot I} D_X(2^j, k)\}_{k \in \mathbb{Z}}. \quad (C.94) $$

(ii) The wavelet random matrix $W_X(a(n)2^j)$ as in (2.16) satisfies a random operator self-similarity relation. That is,

$$ W_X(a(n)2^j) \stackrel{d}{=} a(n)^{\mathbb{H}_n+1/2 \cdot I} W_X(2^j) a(n)^{\mathbb{H}_n+1/2 \cdot I}. \quad (C.95) $$

*Proof.* By (2.12) and Lemma C.14, (ii), the wavelet transform of the stochastic process $X$ is given by

$$ D_X(a(n)2^j, k) = (a(n)2^j)^{1/2} \int_{\mathbb{R}} \psi(u) X(a(n)2^j(u+k)) du, $$

where we make the change of variables $u = (a(n)2^j)^{-1}t - k$. However, conditionally on $\mathbb{H}_n = \tilde{\mathbb{H}}_n$, $X$ is an ofBm with Hurst (matrix) exponent $\tilde{\mathbb{H}}_n$. Therefore, making use of the operator self-similarity of $X$ (see (1.3)), we obtain

$$ \{D_X(a(n)2^j, k) | \mathbb{H}_n = \tilde{\mathbb{H}}_n\}_{k \in \mathbb{Z}} \stackrel{\text{f.d.d.}}{=} \{a(n)^{\mathbb{H}_n+1/2 \cdot I} 2^{j/2} \int_{\mathbb{R}} \psi(u) X(2^j(u+k)) du\}_{k \in \mathbb{Z}}. $$

Thus, by a second change of variables $t = 2^j(u+k)$, we obtain

$$ \{D_X(a(n)2^j, k) | \mathbb{H}_n = \tilde{\mathbb{H}}_n\}_{k \in \mathbb{Z}} \stackrel{\text{f.d.d.}}{=} \{a(n)^{\mathbb{H}_n+1/2 \cdot I} D_X(2^j, k) | \mathbb{H}_n = \tilde{\mathbb{H}}_n\}_{k \in \mathbb{Z}}. \quad (C.96) $$

Now, for any $m \in \mathbb{N}$, consider a collection of sets $B_1, \dots, B_m \in \mathcal{B}(\mathbb{R}^p)$ and integers $k_1, \dots, k_m$. Then,

$$ \begin{aligned} & P\left(\bigcap_{\ell=1}^{m} \{D_X(a(n)2^j, k_\ell) \in B_\ell\}\right) \\ &= \int_{S_{\ge 0}(p(n), \mathbb{R})} P\left(\bigcap_{\ell=1}^{m} \{D_X(a(n)2^j, k_\ell) \in B_\ell\} \middle| \mathbb{H}_n = \tilde{\mathbb{H}}_n\right) P_{\mathbb{H}_n}(d\tilde{\mathbb{H}}_n) \\ &= \int_{S_{\ge 0}(p(n), \mathbb{R})} P\left(\bigcap_{\ell=1}^{m} \{a(n)^{\mathbb{H}_n+1/2 \cdot I} D_X(2^j, k_\ell) \in B_\ell\} \middle| \mathbb{H}_n = \tilde{\mathbb{H}}_n\right) P_{\mathbb{H}_n}(d\tilde{\mathbb{H}}_n) \\ &= P\left(\bigcap_{\ell=1}^{m} \{a(n)^{\mathbb{H}_n+1/2 \cdot I} D_X(2^j, k_\ell) \in B_\ell\}\right), \end{aligned} $$

where the second equality follows from (C.96). This establishes (i). Moreover, since

$$ W_X(a(n)2^j) = \frac{1}{n_{a,j}} \sum_{k=1}^{n_{a,j}} D_X(a(n)2^j, k) D_X(a(n)2^j, k)^*, $$

then (ii) is a direct consequence of (i).

□

The following result is a technical lemma used in the proof of Lemma C.13.

**Lemma C.17.** Under the assumptions of Lemma C.13,

$$ \sum_{\ell=-\infty}^{\infty} |\hat{\psi}(x + 2\pi\ell)|^2 = 1 \quad \text{for all } x \in [-\pi, \pi]. \quad (C.97) $$
**Proof.** Recall that the orthonormality of {$\psi(t - k)$}$_{k \in \mathbb{Z}}$ is equivalent to the condition that $\sum_{\ell=-\infty}^{\infty} |\hat{\psi}(x + 2\pi\ell)|^2 = 1$ a.e. (see Daubechies [35], p. 134). From condition (A.3) (see (W3)), there exist constants $0 < A < \infty$ and $\alpha > 1$ such that, for each $x \in \mathbb{R}$ and $\ell \in \mathbb{Z}$, $|\hat{\psi}(x + 2\pi\ell)|^2 \le A^2/(1 + |x + 2\pi\ell|)^{2\alpha}$. Since $\hat{\psi}$ is continuous (see conditions (A.1) and (A.2) in (W1) and (W2), then, by the dominated convergence theorem, the function $x \mapsto \sum_{\ell=-\infty}^{\infty} |\hat{\psi}(x + 2\pi\ell)|^2$ is continuous. Hence, (C.97) holds. $\square$

The following lemma is used in the proof of Lemma B.5. In particular, for $X_H$ as in (2.3), it establishes properties of the conditional covariance matrix of the discrete-time wavelet coefficients as given by (2.19), namely,

$$ \mathbb{R}^{n_{a,j}} \ni \tilde{d}_{\mathcal{H}}(2^j) := \{\tilde{d}_{\mathcal{H}}(2^j, k)\}_{k=1,\dots,n_{a,j}} := \left\{ \sum_{\ell \in \mathbb{Z}} X_{\mathcal{H}}(\ell) h_{j,2^j k-\ell} \right\}_{k=1,\dots,n_{a,j}}. \quad (C.98) $$

**Lemma C.18.** Suppose assumptions (A1)-(A5) and (W1)-(W3) hold. For $j \in \mathbb{N} \cup \{0\}$, let $\tilde{d}_{\mathcal{H}}(2^j, k)$ be a univariate discrete-time wavelet coefficient as in (C.98). In addition, for any $H \in (0, 1)$, let

$$ \tilde{\Sigma}_H = \tilde{\Sigma}_{H,n_{a,j}} = \left\{ \mathrm{E} \left[ \tilde{d}_{\mathcal{H}}(2^j, k_1) \tilde{d}_{\mathcal{H}}(2^j, k_2) \mid \mathcal{H} = H \right] \right\}_{k_1,k_2=1,\dots,n_{a,j}} \in \mathcal{S}_{\ge 0}(n_{a,j}, \mathbb{R}) \quad (C.99) $$

be the conditional covariance matrix of the discrete-time wavelet coefficients in (C.62). Then, $\tilde{\Sigma}_H$ has full rank, and conditionally on $\mathcal{H} = H$,

$$ \tilde{d}_{\mathcal{H}}(2^j) |_{\mathcal{H} = H} \stackrel{\text{a.s.}}{=} \tilde{\Gamma}_H \mathbf{Z}, $$

where $\mathbf{Z} \in \mathbb{R}^{n_{a,j}}$ is a standard normal random vector and $\tilde{\Gamma}_H^2 = \tilde{\Sigma}_H$.

**Proof.** By definition, conditionally on $\mathcal{H} = H$, $X_H$ is a fBm with Hurst exponent $H$. Hence, also conditionally, the random vector $\tilde{d}_{\mathcal{H}}(2^j)$ is Gaussian. By expressions (A.10) and (A.11) in Boniece et al. [22], $\mathrm{E}[\tilde{d}_{\mathcal{H}}^2(2^j, 0)|\mathcal{H} = H] > 0$. Equivalently, the main diagonal entries of $\tilde{\Sigma}_H$ are (equal and) strictly positive. Moreover, $\mathrm{E}[\tilde{d}_{\mathcal{H}}(2^j, \kappa)\tilde{d}_{\mathcal{H}}(2^j, 0)|\mathcal{H} = H] \to 0$ as $\kappa \to \infty$ (see, for instance, Boniece and Didier [22], Proposition A.3, (ii)). Therefore, by Proposition 5.1.1 in Brockwell and Davis [25], the matrix $\tilde{\Sigma}_H$ as in (C.99) is nonsingular for every $n_{a,j}$. This completes the proof. $\square$

## References

[1] Abry, P., Boniece, B. C., Didier, G. and Wendt, H. (2024), 'On high-dimensional wavelet eigenanalysis', *Annals of Applied Probability* **34** (6), pp. 5287–5350. MR4840487

[2] Abry, P., Boniece, B. C., Didier, G. and Wendt, H. (2023), 'Wavelet eigenvalue regression in high dimensions', *Statistical Inference for Stochastic Processes* **26**(1), 1–32. MR4562251

[3] Abry, P. and Didier, G. (2018a), 'Wavelet eigenvalue regression for $n$-variate operator fractional Brownian motion', *Journal of Multivariate Analysis* **168**, 75–104. MR3858352

[4] Abry, P. and Didier, G. (2018b), 'Wavelet estimation for operator fractional Brownian motion', *Bernoulli* **24**(2), 895–928. MR3706780

[5] Abry, P. and Didier, G. and Li, H. (2019), 'Two-step wavelet-based estimation for Gaussian mixed fractional processes', *Statistical Inference for Stochastic Processes* **22**, 157–185. MR3959286

[6] Abry, P. and Flandrin, P. (1994), 'On the initialization of the discrete wavelet transform algorithm', *IEEE Signal Processing Letters* **1**(2), 32–34.

[7] Abry, P., Flandrin, P., Taqqu, M. S. and Veitch, D. (2003), Self-similarity and long-range dependence through the wavelet lens, in 'Theory and Applications of Long-Range Dependence', Vol. 1, Birkhäuser, Boston, pp. 527–556. MR1957507
[8] Alshammri, F. and Pan, J. (2021), 'Generalized principal component analysis for moderately non-stationary vector time series', *Journal of Statistical Planning and Inference* **212**, 201–225. MR4180112

[9] Anderson, G., Guionnet, A. and Zeitouni, O. (2010), *An Introduction to Random Matrices*, Cambridge Studies in Advanced Mathematics, volume 118, Cambridge University Press. MR2760897

[10] Anderson, T. (2003), *An Introduction to Multivariate Statistical Analysis*, 3 edn, Wiley. MR1990662

[11] Bai, S. and Taqqu, M. S. (2018), 'How the instability of ranks under long memory affects large-sample inference', *Statistical Science* **33**(1), 96–116. MR3757507

[12] Bai, Z. D. and Yin, Y. Q. (1988), 'Convergence to the semicircle law', *Annals of Probability* **16**(2), 863–875. MR0929083

[13] Bai, Z. and Silverstein, J. (2010), *Spectral Analysis of Large Dimensional Random Matrices*, Vol. 20, Springer. MR2567175

[14] Baik, J., Deift, P. and Johansson, K. (1999), 'On the distribution of the length of the longest increasing subsequence of random permutations', *Journal of the American Mathematical Society* **12**(4), 1119–1178. MR1682248

[15] Baik, J., Deift, P. and Suidan, T. (2016), *Combinatorics and Random Matrix Theory*, Vol. 172 of *Graduate Studies in Mathematics*, American Mathematical Society. MR3468920

[16] Baik, J., Kriecherbauer, T., McLaughlin, K. and Miller, P. D. (2007), *Discrete Orthogonal Polynomials: Asymptotics and Applications*, Vol. 164, Princeton University Press. MR2283089

[17] Balcerek, M., Burnecki, K., Thapa, S., Wylomanska, A. and Chechkin, A. (2022), 'Fractional Brownian motion with random Hurst exponent: accelerating diffusion and persistence transitions', *Chaos* **32**, 093114. MR4480621

[18] Bardet, J.-M. (2002), 'Statistical study of the wavelet analysis of fractional Brownian motion', *IEEE Transactions on Information Theory* **48**(4), 991–999. MR1908463

[19] Becker-Kern, P. and Pap, G. (2008), 'Parameter estimation of selfsimilarity exponents', *Journal of Multivariate Analysis* **99**(1), 117–140. MR2408446

[20] Ben Arous, G. and Guionnet, A. (1997), 'Large deviations for Wigner's law and Voiculescu's non-commutative entropy', *Probability Theory and Related Fields* **108**(4), 517–542. MR1465640

[21] Benson, D. A., Meerschaert, M. M., Baeumer, B. and Scheffler, H.-P. (2006), 'Aquifer operator scaling and the effect on solute mixing and dispersion', *Water Resources Research* **42**(1).

[22] Boniece, B. C., Didier, G. and Sabzikar, F. (2021), 'Tempered fractional Brownian motion: wavelet estimation, modeling and testing', *Applied and Computational Harmonic Analysis* **51**, 461–509. MR4196450

[23] Borodin, A. and Petrov, L. (2014), 'Integrable probability: from representation theory to Macdonald processes', *Probability Surveys* **11**, 1–58. MR3178541

[24] Bourguin, S., Diez, C.-P. and Tudor, C. A. (2021), 'Limiting behavior of large correlated Wishart matrices with chaotic entries', *Bernoulli* **27**(2), 1077–1102. MR4255227

[25] Brockwell, P. and Davis, R. (1991), *Time Series: Theory and Methods*, 2nd edn, Springer. MR1093459

[26] Chakrabarty, A., Hazra, R. S. and Sarkar, D. (2016), 'From random matrices to long range dependence', *Random Matrices: Theory and Applications* **5**(02), 1650008. MR3493554

[27] Chan, N. H., Lu, Y. and Yau, C. Y. (2017), 'Factor modelling for high-dimensional time series: inference and model selection', *Journal of Time Series Analysis* **38**(2), 285–307. MR3611745

[28] Chang, J. and Guo, B. and Yao, Q. (2018), 'Principal component analysis for second-order stationary vector time series', *Annals of Statistics* **46**(5), 2094–2124. MR3845012

[29] Che, Z. (2017), 'Universality of random matrices with correlated entries', *Electronic Journal of Probability* **22**, 1–38. MR3629874

[30] Chen, B. and Pan, G. (2012), 'Convergence of the largest eigenvalue of normalized sample covariance matrices when p and n both tend to infinity with their ratio converging to zero', *Bernoulli* **18**(4), 1405–1420. MR2995802
[31] Cohen, A. (2003), *Numerical Analysis of Wavelet Methods*, Vol. 32, North-Holland, Amsterdam. MR1990555

[32] Comon, P. and Jutten, C. (2010), *Handbook of Blind Source Separation: Independent Component Analysis and Applications*, Academic Press.

[33] Craigmile, P., Guttorp, P. and Percival, D. (2005), 'Wavelet-based parameter estimation for polynomial contaminated fractionally differenced processes', *IEEE Transactions on Signal Processing* **53**(8), 3151-3161. MR2169658

[34] Cramér, H. and Leadbetter, M. R. (1967), *Stationary and Related Stochastic Processes: Sample Function Properties and their Applications*, Courier Dover Publications. MR0217860

[35] Daubechies, I. (1992), *Ten Lectures on Wavelets*, Vol. 61, Society for Industrial and Applied Mathematics, Philadelphia-PA. MR1162107

[36] Deift, P. (2007), *Universality for mathematical and physical systems, in 'Proceedings of the International Congress of Mathematicians Madrid, August 22-30, 2006', pp. 125-152*. MR2334189

[37] Deift, P. (2017), 'Some open problems in random matrix theory and the theory of integrable systems. ii', *SIGMA. Symmetry, Integrability and Geometry: Methods and Applications* **13**, 016. MR3622647

[38] Didier, G. and Pipiras, V. (2012), 'Exponents, symmetry groups and classification of operator fractional Brownian motions', *Journal of Theoretical Probability* **25**, 353-395. MR2914433

[39] Doukhan, P., Oppenheim, G. and Taqqu, M. (2003), *Theory and Applications of Long-Range Dependence*, Birkhäuser, Boston, MA. MR1956041

[40] Dyson, F. J. (1962), 'A Brownian-motion model for the eigenvalues of a random matrix', *Journal of Mathematical Physics* **3**(6), 1191-1198. MR0148397

[41] Embrechts, P. and Maejima, M. (2002), *Selfsimilar Processes*, Princeton Series in Applied Mathematics, Princeton University Press, Princeton, NJ. MR1920153

[42] Engle, R. F. and Granger, C. W. J. (1987), 'Co-integration and error correction: representation, estimation, and testing', *Econometrica* pp. 251-276. MR0882095

[43] Erdős, L., Krüger, T. and Schröder, D. (2019), 'Random matrices with slow correlation decay', *Forum of Mathematics, Sigma* **7**, e8, 89 pages. MR3941370

[44] Erdős, L., Yau, H.-T. and Yin, J. (2012), 'Rigidity of eigenvalues of generalized Wigner matrices', *Advances in Mathematics* **229**(3), 1435-1515. MR2871147

[45] Flandrin, P. (1992), 'Wavelet analysis and synthesis of fractional Brownian motion', *IEEE Transactions on Information Theory* **38**, 910-917. MR1162229

[46] Gotts, S. J., Gilmore, A. W. and Martin, A. (2020), 'Brain networks, dimensionality, and global signal averaging in resting-state fMRI: Hierarchical network structure results in low-dimensional spatiotemporal dynamics', *Neuroimage* **205**, 116289.

[47] Gray, R. M. (2006), *Toeplitz and Circulant Matrices: a Review*, Now Publishers, Inc.

[48] He, J.-H. (2018), 'Fractal calculus and its geometrical explanation', *Results in Physics* **10**, 272-276.

[49] Horn, R. A. and Johnson, C. R. (2012), *Matrix Analysis*, 2 edn, Cambridge University Press. MR2978290

[50] Hualde, J. and Robinson, P. M. (2010), 'Semiparametric inference in multivariate fractionally cointegrated systems', *Journal of Econometrics* **157**(2), 492-511. MR2661617

[51] Hudson, W. and Mason, J. (1982), 'Operator-self-similar processes in a finite-dimensional space', *Transactions of the American Mathematical Society* **273**(1), 281-297. MR0664042

[52] Isotta, F., Frei, C., Weilguni, V., Perčec Tadić, M., Lassegues, P., Rudolf, B., Pavan, V., Cacciamani, C., Antolini, G., Ratto, S. M. and Munari, M. (2014), 'The climate of daily precipitation in the Alps: development and analysis of a high-resolution grid dataset from pan-Alpine rain-gauge data', *International Journal of Climatology* **34**(5), 1657-1675.

[53] Jiang, T. (2008), 'Approximation of Haar distributed matrices and limiting distributions of eigenvalues of Jacobi ensembles', *Probability Theory and Related Fields* **144**(1-2), 221-246. MR2480790
[54] Kolmogorov, A. N. (1940), The Wiener spiral and some other interesting curves in Hilbert space, in 'Dokl. Akad. Nauk SSSR', Vol. 26, pp. 115-118. MR0009098

[55] Kolmogorov, A. N. (1941), The local structure of turbulence in an incompressible fluid at very high Reynolds numbers, in 'Dokl. Akad. Nauk SSSR', Vol. 30, pp. 299-303. MR0004146

[56] Kuijlaars, A. and McLaughlin, K. (2000), 'Generic behavior of the density of states in random matrix theory and equilibrium problems in the presence of real analytic external fields', *Communications on Pure and Applied Mathematics* **53**(6), 736-785. MR1744002

[57] Laha, R. G. and Rohatgi, V. K. (1981), 'Operator self-similar stochastic processes in $R^d$, *Stochastic Processes and their Applications* **12**(1), 73-84. MR0632393

[58] Li, L., Pluta, D., Shahbaba, B., Fortin, N., Ombao, H. and Baldi, P. (2019), 'Modeling dynamic functional connectivity with latent factor Gaussian processes', *Advances in Neural Information Processing Systems* **32**, 8263-8273.

[59] Li, Q., Pan, J. and Yao, Q. (2009), 'On determination of cointegration ranks', *Statistics and Its Interface* **2**(1), 45-56. MR2500767

[60] Maejima, M. and Mason, J. (1994), 'Operator-self-similar stable processes', *Stochastic Processes and their Applications* **54**, 139-163. MR1302699

[61] Mallat, S. (1999), *A Wavelet Tour of Signal Processing*, Academic Press. MR1614527

[62] Mandelbrot, B. (1982), *The Fractal Geometry of Nature*, W. H. Freeman, San Francisco.
MR0665254

[63] Mandelbrot, B. and Van Ness, J. (1968), 'Fractional Brownian Motions, Fractional Noises and Applications', *SIAM Review* **10**(4), 422-437. MR0242239

[64] Mason, J. and Xiao, Y. (2002), 'Sample path properties of operator-self-similiar Gaussian random fields', *Theory of Probability and Its Applications* **46**(1), 58-78. MR1968707

[65] Massopust, P. R. (2014), *Fractal Functions, Fractal Surfaces, and Wavelets*, Academic Press, San Diego, CA. MR3642696

[66] Meckes, E. (2014), *Concentration of Measure and the Compact Classical Matrix Groups*, Lecture Notes (IAS Program for Women and Mathematics).

[67] Meerschaert, M. and Scheffler, H.-P. (1999), 'Moment estimator for random vectors with heavy tails', *Journal of Multivariate Analysis* **71**(1), 145-159. MR1721965

[68] Meerschaert, M. and Scheffler, H.-P. (2003), Portfolio modeling with heavy-tailed random vectors, in 'Handbook of heavy-tailed distributions in Finance (S. T. Rachev (ed.))', Elsevier Science B.V., Amsterdam, pp. 595-640.

[69] Mehta, M. and Gaudin, M. (1960), 'On the density of eigenvalues of a random matrix', *Nuclear Physics* **18**, 420-427. MR0112895

[70] Mehta, M. L. (2004), *Random Matrices*, 3 edn, Elsevier. MR2129906

[71] Menon, G. and Trogdon, T. (2016), 'Smoothed analysis for the conjugate gradient algo-
rithm', SIGMA. Symmetry, Integrability and Geometry: Methods and Applications **12**, 109.
MR3568010

[72] Menon, G. and Trogdon, T. (2023), *Random Matrix Theory and Numerical Algorithms*, https://www.dam.brown.edu/people/menon/publications/rmt-2021.pdf.

[73] Merikoski, J. K. and Kumar, R. (2004), 'Inequalities for spreads of matrix sums and products', *Applied Mathematics E-Notes* **4**, 150-159. MR2077795

[74] Merlevède, F., Najim, J. and Tian, P. (2019), 'Unbounded largest eigenvalue of large sample covariance matrices: Asymptotics, fluctuations and applications', *Linear Algebra and its Applications* **577**, 317-359. MR3948744

[75] Merlevède, F. and Peligrad, M. (2016), 'On the empirical spectral distribution for matrices with long memory and independent rows', *Stochastic Processes and their Applications* **126**(9), 2734-2760. MR3522299

[76] Meyer, Y., Sellan, F. and Taqqu, M. S. (1999), 'Wavelets, generalized white noise and fractional integration: the synthesis of fractional Brownian motion', *Journal of Fourier Analysis and Applications* **5**(5), 465-494. MR1755100
[77] Moulines, E., Roueff, F. and Taqqu, M. (2007a), 'Central limit theorem for the log-regression wavelet estimation of the memory parameter in the gaussian semi-parametric context', *Fractals* **15**(4), 301–313. MR2396718

[78] Moulines, E., Roueff, F. and Taqqu, M. (2007b), 'On the spectral density of the wavelet coefficients of long-memory time series with application to the log-regression estimation of the memory parameter', *Journal of Time Series Analysis* **28**(2), 155–187. MR2345656

[79] Moulines, E., Roueff, F. and Taqqu, M. (2008), 'A wavelet Whittle estimator of the memory parameter of a nonstationary Gaussian time series', *Annals of Statistics*, **36**(4) pp. 1925–1956. MR2435460

[80] NobelPrize.org (2003), 'Clive W.J. Granger: Facts'. https://www.nobelprize.org/prizes/economic-sciences/2003/granger/facts/

[81] Orejola, O., Didier, G., Wendt, H. and Abry, P. (2022), Hurst multimodality detection based on large wavelet random matrices, in '30th European Signal Processing Conference (EUSIPCO), Belgrade, Serbia', pp. 1–5.

[82] Paul, D. and Aue, A. (2014), 'Random matrix theory in statistics: a review', *Journal of Statistical Planning and Inference* **150**, 1–29. MR3206718

[83] Percival, D. B. and Walden, A. (2006), *Wavelet Methods for Time Series Analysis*, Vol. 4, Cambridge University Press. MR2218866

[84] Phillips, P. C. B. and Ouliaris, S. (1988), 'Testing for cointegration using principal components methods', *Journal of Economic Dynamics and Control* **12**(2-3), 205–230. MR0986515

[85] Pipiras, V. and Taqqu, M. S. (2017), *Long-Range Dependence and Self-Similarity*, Cambridge Series on Statistical and Probabilistic Mathematics, Cambridge University Press, Cambridge, United Kingdom. MR3729426

[86] Reif, F. (2009), *Fundamentals of Statistical and Thermal Physics*, Waveland Press.

[87] Roueff, F. and Taqqu, M. S. (2009), 'Asymptotic normality of wavelet estimators of the memory parameter for linear processes', *Journal of Time Series Analysis* **30**(5), 534–558. MR2560417

[88] Samorodnitsky, G. and Taqqu, M. (1994), *Stable Non-Gaussian Processes: Stochastic Models with Infinite Variance*, Chapman & Hall, New York, London. MR1280932

[89] Schmith, T., Johansen, S. and Thejll, P. (2012), 'Statistical analysis of global surface temperature and sea level using cointegration methods', *Journal of Climate* **25**(22), 7822–7833.

[90] Shen, J., Stoev, S. and Hsing, T. (2022), 'Tangent fields, intrinsic stationarity and self similarity', *Electronic Journal of Probability* **27**, 1–56. MR4387842

[91] Shimotsu, K. (2012), 'Exact local Whittle estimation of fractionally cointegrated systems', *Journal of Econometrics* **169**(2), 266–278. MR2947884

[92] Shiryaev, A. N. (2015), *Probability*, Vol. 1, 3rd edn, Springer. MR3467826

[93] Sornette, D. (2006), *Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforganization and Disorder: Concepts and Tools*, Springer Science & Business Media. MR2036307

[94] Soshnikov, A. (1999), 'Universality at the edge of the spectrum in Wigner random matrices', *Communications in Mathematical Physics* **207**(3), 697–733. MR1727234

[95] Steland, A. and Von Sachs, R. (2017), 'Large-sample approximations for variance-covariance matrices of high-dimensional time series', *Bernoulli* **23**(4A), 2299–2329. MR3648032

[96] Stoev, S., Pipiras, V. and Taqqu, M. (2002), 'Estimation of the self-similarity parameter in linear fractional stable motion', *Signal Processing* **82**, 1873–1901.

[97] Tao, T. (2012), 'E pluribus unum: from complexity, universality', *Daedalus* **141**(3), 23–34. MR3287630

[98] Tao, T. and Vu, V. (2011), 'Random matrices: universality of local eigenvalue statistics', *Acta Mathematica* **206**(1), 127–204. MR2784665

[99] Tao, T. and Vu, V. (2012), 'Random covariance matrices: universality of local statistics of eigenvalues', *Annals of Probability* **40**(3), 1285–1315. MR2962092

[100] Taqqu, M. S. (2003), *Fractional Brownian motion and long range dependence*, in 'Theory and Applications of Long-Range Dependence' (P. Doukhan, G. Oppenheim and M. S. Taqqu, eds.)', Birkhäuser, Boston, pp. 5–38. MR1956042
On the e.s.d. of large wavelet random matrices

[101] Tian, P. (2022), 'Joint CLT for top eigenvalues of sample covariance matrices of separable high dimensional long memory processes', *Random Matrices: Theory and Applications* **11**(03), 2250032. MR4440258

[102] Vershynin, R. (2018), *High-Dimensional Probability: an Introduction with Applications in Data Science*, Vol. 47, Cambridge University Press. MR3837109

[103] Wang, L., Aue, A. and Paul, D. (2017), 'Spectral analysis of linear time series in moderately high dimensions', *Bernoulli* **23**(4A), 2181–2209. MR3648029

[104] West, G. B., Brown, J. H. and Enquist, B. J. (1999), 'The fourth dimension of life: fractal geometry and allometric scaling of organisms', *Science* **284**(5420), 1677–1679. MR1693406

[105] Wornell, G. and Oppenheim, A. (1992), 'Estimation of fractal signals from noisy measurements using wavelets', *IEEE Transactions on Signal Processing* **40**(3), 611–623.

[106] Wornell, G. W. (1996), *Signal Processing with Fractals: A Wavelet-Based Approach*, Prentice-Hall, Englewood Cliffs, NJ.

[107] Xia, N., Qin, Y. and Bai, Z. (2013), 'Convergence rates of eigenvector empirical spectral distribution of large dimensional sample covariance matrix', *Annals of Statistics* **41**(5), 2572–2607. MR3161438

[108] Yaskov, P. (2018), 'LLN for quadratic forms of long memory time series and its applications in random matrix theory', *Journal of Theoretical Probability* **31**(5), 2032–2055. MR3866606

[109] Zhang, D. and Wu, W. B. (2017), 'Gaussian approximation for high dimensional time series', *Annals of Statistics* **45**(5), 1895–1919. MR3718156

[110] Zhang, R., Robinson, P. M. and Yao, Q. (2019), 'Identifying cointegration by eigenanalysis', *Journal of the American Statistical Association* **114**(526), 916–927. MR3963191

[111] Zheng, X., Shen, G., Wang, C., Li, Y., Dunphy, D., Hasan, T., Brinker, C. J. and Su, B.-L. (2017), 'Bio-inspired murray materials for mass transfer and activity', *Nature Communications* **8**(1), 1–9.

**Acknowledgments.** The authors would like to thank Jean-Marc Lina, Ken McLaughlin and Govind Menon for their many helpful comments on this work.