*Article* **Ordinal Pattern Dependence in the Context of Long-Range Dependence**

**Ines Nüßgen \* and Alexander Schnurr**

Department of Mathematics, Siegen University, Walter-Flex-Straße 3, 57072 Siegen, Germany; schnurr@mathematik.uni-siegen.de

**\*** Correspondence: nuessgen@mathematik.uni-siegen.de

**Abstract:** Ordinal pattern dependence is a multivariate dependence measure based on the comovement of two time series. In strong connection to ordinal time series analysis, the ordinal information is taken into account to derive robust results on the dependence between the two processes. This article deals with ordinal pattern dependence for a long-range dependent time series including mixed cases of short- and long-range dependence. We investigate the limit distributions for estimators of ordinal pattern dependence. In doing so, we point out the differences that arise for the underlying time series having different dependence structures. Depending on these assumptions, central and non-central limit theorems are proven. The limit distributions for the latter ones can be included in the class of multivariate Rosenblatt processes. Finally, a simulation study is provided to illustrate our theoretical findings.

**Keywords:** ordinal patterns; time series; long-range dependence; multivariate data analysis; limit theorems

**Citation:** Nüßgen, I.; Schnurr, A. Ordinal Pattern Dependence in the Context of Long-Range Dependence. *Entropy* **2021**, *23*, 670. https:// doi.org/10.3390/e23060670

Academic Editor: Christian H. Weiss

Received: 28 April 2021 Accepted: 19 May 2021 Published: 26 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

The origin of the concept of ordinal patterns is in the theory of dynamical systems. The idea is to consider the order of the values within a data vector instead of the full metrical information. The ordinal information is encoded as a permutation (cf. Section 3). Already in the first papers on the subject, the authors considered entropy concepts related to this ordinal structure (cf. [1]). There is an interesting relationship between these concepts and the well-known Komogorov–Sinai entropy (cf. [2,3]). Additionally, an ordinal version of the Feigenbaum diagram has been dealt with e.g., in [4]. In [5], ordinal patterns were used in order to estimate the Hurst parameter in long-range dependent time series. Furthermore, Ref. [6] have proposed a test for independence between time series (cf. also [7]). Hence, the concept made its way into the area of statistics. Instead of long patterns (or even letting the pattern length tend to infinity), rather short patterns have been considered in this new framework. Furthermore, ordinal patterns have been used in the context of ARMA processes [8] and change-point detection within one time series [9]. In [10], ordinal patterns were used for the first time in order to analyze the dependence between two time series. Limit theorems for this new concept were proven in a short-range dependent framework in [11]. Ordinal pattern dependence is a promising tool, which has already been used in financial, biological and hydrological data sets, see in this context, also [12] for an analysis of the co-movement of time series focusing on symbols. In particular, in the context of hydrology, the data sets are known to be long-range dependent. Therefore, it is important to also have limit theorems available in this framework. We close this gap in the present article.

All of the results presented in this article have been established in the Ph.D. thesis of I. Nüßgen written under the supervision of A. Schnurr.

The article is structured as follows: in the subsequent section, we provide the reader with the mathematical framework. The focus is on (multivariate) long-range dependence. In Section 3, we recall the concept of ordinal pattern dependence and prove our main

results. We present a simulation study in Section 4 and close the paper by a short outlook in Section 5.

#### **2. Mathematical Framework**

We consider a stationary *d*-dimensional Gaussian time series - *Yj <sup>j</sup>*∈<sup>Z</sup> (for *<sup>d</sup>* <sup>∈</sup> <sup>N</sup>), with:

$$\mathcal{Y}\_{\vec{j}} := \left( \mathcal{Y}\_{\vec{j}}^{(1)}, \dots, \mathcal{Y}\_{\vec{j}}^{(d)} \right)^{t} \tag{1}$$

such that E *<sup>Y</sup>*(*p*) *j* = 0 and E *<sup>Y</sup>*(*p*) *j* 2 <sup>=</sup> 1 for all *<sup>j</sup>* <sup>∈</sup> <sup>Z</sup> and *<sup>p</sup>* <sup>=</sup> 1, ... , *<sup>d</sup>*. Furthermore, we require the cross-correlation function to fulfil *r*(*p*,*q*)(*k*) <sup>&</sup>lt; 1 for *<sup>p</sup>*, *<sup>q</sup>* <sup>=</sup> 1, ... , *<sup>d</sup>* and *<sup>k</sup>* <sup>≥</sup> 1, where the component-wise cross-correlation functions *r*(*p*,*q*)(*k*) are given by *r*(*p*,*q*)(*k*) = E *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j*+*k* for each *<sup>p</sup>*, *<sup>q</sup>* <sup>=</sup> 1, ... , *<sup>d</sup>* and *<sup>k</sup>* <sup>∈</sup> <sup>Z</sup>. For each random vector *Yj*, we denote the covariance matrix by Σ*d*, since it is independent of *j* due to stationarity. Therefore, we have Σ*<sup>d</sup>* = *r*(*p*,*q*)(0) *p*,*q*=1,...,*d* .

We specify the dependence structure of - *Yj <sup>j</sup>*∈<sup>Z</sup> and turn to long-range dependence: we assume that for the cross-correlation function *r*(*p*,*q*)(*k*) for each *p*, *q* = 1, ... , *d*, it holds that:

$$r^{(p,q)}(k) = L\_{p,q}(k)k^{d\_p + d\_q - 1},\tag{2}$$

with *Lp*,*q*(*k*) → *Lp*,*<sup>q</sup>* (*k* → ∞) for finite constants *Lp*,*<sup>q</sup>* ∈ [0, ∞) with *Lp*,*<sup>p</sup>* = 0, where the matrix *L* = - *Lp*,*q <sup>p</sup>*,*q*=1,...,*<sup>d</sup>* has full rank, is symmetric and positive definite. Furthermore, the parameters *dp*, *dq* ∈ 0, <sup>1</sup> 2 are called long-range dependence parameters. Therefore, - *Yj <sup>j</sup>*∈<sup>Z</sup> is multivariate long-range dependent in the sense of [13], Definition 2.1.

The processes we want to consider have a particular structure, namely for *<sup>h</sup>* <sup>∈</sup> <sup>N</sup>, we obtain for fixed *<sup>j</sup>* <sup>∈</sup> <sup>Z</sup>:

$$\mathcal{Y}\_{j,\mathbb{h}} := \left( \mathcal{Y}\_{j}^{(1)}, \dots, \mathcal{Y}\_{j+\mathfrak{h}-1}^{(1)}, \mathcal{Y}\_{j}^{(2)}, \dots, \mathcal{Y}\_{j+\mathfrak{h}-1}^{(2)}, \dots, \mathcal{Y}\_{j}^{(d)}, \dots, \mathcal{Y}\_{j+\mathfrak{h}-1}^{(d)} \right)^{\mathfrak{t}} \in \mathbb{R}^{d\mathfrak{h}}.\tag{3}$$

The following relation holds between the *extendend process Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> and the primarily regarded process - *Yj <sup>j</sup>*∈Z. For all *<sup>k</sup>* <sup>=</sup> 1, . . . , *dh*, *<sup>j</sup>* <sup>∈</sup> <sup>Z</sup> we have:

$$Y\_{j,h}^{(k)} = Y\_{j+(k \mod h)-1'}^{\left\lfloor \frac{k-1}{h} \right\rfloor + 1} \tag{4}$$

where *x* <sup>=</sup> max{*<sup>k</sup>* <sup>∈</sup> <sup>Z</sup> : *<sup>k</sup>* <sup>≤</sup> *<sup>x</sup>*}. Note that the process *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> is still a centered Gaussian process since all finite-dimensional marginals of - *Yj <sup>j</sup>*∈<sup>Z</sup> follow a normal distribution. Stationarity is also preserved since for all *<sup>p</sup>*, *<sup>q</sup>* <sup>=</sup> 1, ... , *dh*, *<sup>p</sup>* <sup>≤</sup> *<sup>q</sup>* and *<sup>k</sup>* <sup>∈</sup> <sup>Z</sup>, the cross-correlation function *<sup>r</sup>*(*p*,*q*,*h*)(*k*) of the process *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> is given by

$$\begin{aligned} r^{(p,q,h)}(k) &= \mathbb{E}\left(Y\_{j,h}^{(p)}Y\_{j+k,h}^{(q)}\right) \\ &= \mathbb{E}\left(Y\_{j+(p\mod h)-1}^{\left\lfloor\frac{p-1}{h}\right\rfloor+1}Y\_{j+k+(q\mod h)-1}^{\left\lfloor\frac{q-1}{h}\right\rfloor+1}\right) \\ &= r^{(\left\lfloor\frac{p-1}{h}\right\rfloor+1,\left\lfloor\frac{q-1}{h}\right\rfloor+1)}(k+((q-p)\mod h)) \end{aligned} \tag{5}$$

and the last line does not depend on *j*. The covariance matrix Σ*d*,*<sup>h</sup>* of *Yj*,*<sup>h</sup>* has the following structure:

$$\begin{aligned} \left(\Sigma\_{d,\mathrm{h}}\right)\_{\substack{p,q=1,\ldots,d\_{\star} \\ p\leq q}} &= \left(r^{\left(p,q,\mathrm{h}\right)}\left(0\right)\right)\_{\substack{p,q=1,\ldots,d\mathrm{h}\_{\star}\mathrm{-}\\ p\leq q,\end{pmatrix}}\\ \left(\Sigma\_{d,\mathrm{h}}\right)\_{\substack{p,q=1,\ldots,d\_{\star} \\ p>q}} &= \left(r^{\left(q,p,\mathrm{h}\right)}\left(0\right)\right)\_{\substack{p,q=1,\ldots,d\mathrm{h}\_{\star}\mathrm{-}\\ q$$

Hence, we arrive at:

$$
\Sigma\_{\rm dJ} = \left(\Sigma\_{\rm lt}^{(p,q)}\right)\_{1 \le p, q \le d'} \tag{6}
$$

where Σ(*p*,*q*) *<sup>h</sup>* <sup>=</sup> <sup>E</sup> *<sup>Y</sup>*(*p*) <sup>1</sup> ,...,*Y*(*p*) *h <sup>t</sup> <sup>Y</sup>*(*q*) <sup>1</sup> ,...,*Y*(*q*) *h* = *<sup>r</sup>*(*p*,*q*)(*<sup>i</sup>* <sup>−</sup> *<sup>k</sup>*) 1≤*i*,*k*≤*h* , *p*, *q* = 1, ... , *d*. Note that <sup>Σ</sup>(*p*,*q*) *<sup>h</sup>* <sup>∈</sup> <sup>R</sup>*h*×*<sup>h</sup>* and *<sup>r</sup>*(*p*,*q*)(*k*) = *<sup>r</sup>*(*q*,*p*)(−*k*), *<sup>k</sup>* <sup>∈</sup> <sup>Z</sup> since we are studying crosscorrelation functions.

Therefore, we finally have to show that based on the assumptions on - *Yj <sup>j</sup>*∈Z, the extended process is still long-range dependent.

Hence, we have to consider the cross-correlations again:

$$r^{(p,q,h)}(k) = r^{(\left\lfloor \frac{p-1}{h} \right\rfloor + 1, \left\lfloor \frac{q-1}{h} \right\rfloor + 1)}(k + ((q-p) \mod h))$$

$$= r^{(p^\*, q^\*)}(k + m^\*)$$

$$\simeq r^{(p^\*, q^\*)}(k) \ (k \to \infty),\tag{7}$$

since *p*∗, *q*<sup>∗</sup> ∈ {1, ... , *d*} and *m*<sup>∗</sup> ∈ {0, ... , *h* − 1}, with *p*<sup>∗</sup> := \* *<sup>p</sup>*−<sup>1</sup> *h* + + 1, *q*∗ := \* *<sup>q</sup>*−<sup>1</sup> *h* + + 1 and *m*<sup>∗</sup> = (*q* − *p*) mod *h*.

Let us remark that *ak bk* <sup>⇔</sup> lim*k*→<sup>∞</sup> *ak bk* <sup>=</sup> 1.

Therefore, we are still dealing with a multivariate long-range dependent Gaussian process. We see in the proofs of the following limit theorems that the crucial parameters that determine the asymptotic distribution are the long-range dependence parameters *dp*, *p* = 1, ... , *d* of the original process - *Yj <sup>j</sup>*∈<sup>Z</sup> and therefore, we omit the detailed description of the parameters *dp*<sup>∗</sup> herein.

It is important to remark that the extended process *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> is also long-range dependent in the sense of [14], p. 2259, since:

$$\begin{split} \lim\_{k \to \infty} \frac{k^D r^{(p,q,h)}(k)}{L(k)} &= \lim\_{k \to \infty} \frac{k^D r^{(p^\*,q^\*)}(k)}{L(k)} \\ &= \lim\_{k \to \infty} \frac{k^D L\_{p^\*,q^\*} k^{d\_{p^\*} + d\_{q^\*} - 1}}{L(k)} \\ &=: b\_{p^\*,q^\*} \end{split} \tag{8}$$

with:

$$D := \min\_{p^\* \in \{1, \dots, d\}} \{1 - 2d\_{p^\*}\} \in (0, 1) \tag{9}$$

and *L*(*k*) can be chosen as any constant *Lp*,*<sup>q</sup>* that is not equal to zero, so for simplicity, we assume without a loss of generality *L*1,1 = 0, and therefore, *L*(*k*) = *L*1,1, since the condition in [14] only requires convergence to a finite constant *bp*∗,*q*<sup>∗</sup> . Hence, we may apply the results in [14] in the subsequent results.

We define the following set, which is needed in the proofs of the theorems of this section.

$$P^\* := \{ p \in \{ 1, \ldots, d \} \, : \, d\_{\mathcal{P}} \ge d\_{\mathcal{P}}, \text{ for all } q \in \{ 1, \ldots, d \} \} \tag{10}$$

and denote the corresponding long-range dependence parameter to each *p* ∈ *P*<sup>∗</sup> by

*d*<sup>∗</sup> := *dp*, *p* ∈ *P*∗.

We briefly recall the concept of Hermite polynomials as they play a crucial role in determining the limit distribution of functionals of multivariate Gaussian processes.

**Definition 1.** *(Hermite polynomial, [15], Definition 3.1) The j-th Hermite polynomial Hj*(*x*)*, j* = 0, 1, . . .*, is defined as*

$$H\_{\vec{\jmath}}(\mathbf{x}) := (-1)^{\vec{\jmath}} \exp\left(\frac{\mathbf{x}^2}{2}\right) \frac{\mathbf{d}^{\vec{\jmath}}}{\mathbf{d}x^{\vec{\jmath}}} \exp\left(-\frac{\mathbf{x}^2}{2}\right).$$

Their multivariate extension is given by the subsequent definition.

**Definition 2.** *(Multivariate Hermite polynomial, [15], p. 122) Let d* <sup>∈</sup> <sup>N</sup>*. We define as d-dimensional Hermite polynomial:*

$$H\_k(\mathbf{x}) := H\_{k\_1,\dots,k\_d}(\mathbf{x}) := H\_{k\_1,\dots,k\_d}(\mathbf{x}\_1,\dots,\mathbf{x}\_d) = \prod\_{j=1}^d H\_{k\_j}(\mathbf{x}\_j),$$

*with k* <sup>=</sup> (*k*1,..., *kd*) <sup>∈</sup> <sup>N</sup>*<sup>d</sup>* <sup>0</sup> \ {(0, . . . , 0)}*.*

Let us remark that the case *k* = (0, ... , 0) is excluded here due to the assumption E(*f*(*X*)) = 0.

Analogously to the univariate case, the family of multivariate Hermite polynomials , *Hk*1,...,*kd* , *<sup>k</sup>*1,..., *kd* <sup>∈</sup> <sup>N</sup> - forms an orthogonal basis of *L*<sup>2</sup> <sup>R</sup>*d*, *<sup>ϕ</sup>Id* , which is defined as

$$L^2\left(\mathbb{R}^d, q\boldsymbol{\iota}\_d\right) := \left\{ f : \mathbb{R}^d \to \mathbb{R}, \int\_{\mathbb{R}^d} f^2(\mathbf{x}\_1, \dots, \mathbf{x}\_d) \varrho(\mathbf{x}\_1) \dots \varrho(\mathbf{x}\_d) d\mathbf{x}\_d \dots d\mathbf{x}\_1 < \infty \right\}.$$

The parameter *ϕId* denotes the density of the *d*-dimensional standard normal distribution, which is already divided into the product of the univariate densities *ϕ* in the formula above.

We denote the Hermite coefficients by

$$\mathbb{C}(f, X, k) := \mathbb{C}(f, I\_{d'}k) := \langle f, H\_k \rangle = \mathbb{E}(f(X)H\_k(X))\dots$$

The Hermite rank *m*(*f* , *Id*) of *f* with respect to the distribution N (0, *Id*) is defined as the largest integer *m*, such that:

$$\mathbb{E}\left(f(\mathbf{X})\prod\_{j=1}^{d}H\_{k\_j}\Big(\mathbf{X}^{(j)}\Big)\right) = 0 \text{ for all } 0 < k\_1 + \dots \newline k\_d < m.$$

Having these preparatory results in mind, we derive the multivariate Hermite expansion given by

$$f(X) - \mathbb{E}f(X) = \sum\_{k\_1 + \ldots + k\_d \ge m(f, I\_d)} \frac{\mathbb{C}(f, X, k)}{k\_1! \ldots k\_d!} \prod\_{j=1}^d H\_{k\_j} \left( X^{(j)} \right). \tag{11}$$

We focus on the limit theorems for functionals with Hermite rank 2. First, we introduce the matrix-valued Rosenblatt process. This plays a crucial role in the asymptotics of functionals with Hermite rank 2 applied to multivariate long-range dependent Gaussian processes. We begin with the definition of a multivariate Hermitian–Gaussian random measure *B*˜(d*λ*) with independent entries given by

$$\mathcal{B}(\mathbf{d}\lambda) = \left(\mathcal{B}^{(1)}(\mathbf{d}\lambda), \dots, \mathcal{B}^{(d)}(\mathbf{d}\lambda)\right)^{t},\tag{12}$$

where *B*˜(*p*)(d*λ*) is a univariate Hermitian–Gaussian random measure as defined in [16], Definition B.1.3. The multivariate Hermitian–Gaussian random measure *B*˜(d*λ*) satisfies:

$$\begin{aligned} \mathbb{E}\left(\mathcal{B}(\mathbf{d}\lambda)\right) &= 0, \\ \mathbb{E}\left(\vec{B}(\mathbf{d}\lambda)\vec{B}(\mathbf{d}\lambda)^{\*}\right) &= I\_d \,\mathbf{d}\lambda. \end{aligned}$$

and:

$$\mathbb{E}\left(\mathcal{B}^{(p)}(\mathbf{d}\lambda\_1)\overline{\mathcal{B}^{(q)}(\mathbf{d}\lambda\_2)}\right) = 0, \quad |\lambda\_1| \neq |\lambda\_2|, \quad p, q = 1, \dots, d\_{\nu}$$

where *B*˜(d*λ*)<sup>∗</sup> = *B*(1)(d*λ*),..., *B*(*d*)(d*λ*) denotes the Hermitian transpose of *B*˜(d*λ*). Thus, following [14], Theorem 6, we can state the spectral representation of the matrixvalued Rosenblatt process *Z*2,*H*(*t*), *t* ∈ [0, 1] as

$$Z\_{2,H}(t) = \left(Z\_{2,H}^{(p,q)}(t)\right)\_{p,q=1,\ldots,d}$$

where each entry of the matrix is given by

$$Z\_{2,H}^{(p,q)}(t) = \int\_{\mathbb{R}^2}^{\prime\prime} \frac{\exp\left(it(\lambda\_1 + \lambda\_2)\right) - 1}{i(\lambda\_1 + \lambda\_2)} \tilde{B}^{(p)}(\mathbf{d}\lambda\_1)\tilde{B}^{(q)}(\mathbf{d}\lambda\_2) \dots$$

The double prime in <sup>R</sup><sup>2</sup> excludes the diagonals <sup>|</sup>*λi*<sup>|</sup> <sup>=</sup> *λj* , *i* = *j* in the integration. For details on multiple Wiener-Itô integrals, as can be seen in [17].

The following results were taken from [18], Section 3.2. The corresponding proofs were outsourced to the Appendix A.

**Theorem 1.** *Let* - *Yj <sup>j</sup>*∈<sup>Z</sup> *be a stationary Gaussian process as defined in* (1) *that fulfils* (2) *for dp* ∈ 1 4 , 1 2 *, p* <sup>=</sup> 1..., *d. For h* <sup>∈</sup> <sup>N</sup> *we fix:*

$$\mathcal{Y}\_{j,h} := \left( \mathcal{Y}\_j^{(1)}, \dots, \mathcal{Y}\_{j+h-1}^{(1)}, \dots, \mathcal{Y}\_j^{(d)}, \dots, \mathcal{Y}\_{j+h-1}^{(d)} \right)^t \in \mathbb{R}^{dhh}$$

*with Yj*,*<sup>h</sup>* ∼ N (0, <sup>Σ</sup>*d*,*h*) *and* <sup>Σ</sup>*d*,*<sup>h</sup> as described in* (6)*. Let f* : <sup>R</sup>*dh* <sup>→</sup> <sup>R</sup> *be a function with Hermite rank* 2 *such that the set of discontinuity points Df is a Null set with respect to the dh-dimensional Lebesgue measure. Furthermore, we assume f fulfills* E *f* 2 *Yj*,*<sup>h</sup>* <sup>&</sup>lt; <sup>∞</sup>*. Then:*

$$n^{-2d^\*} \left(\mathbb{C}\_2\right)^{-\frac{1}{2}} \sum\_{j=1}^n \left( f\left(Y\_j^{(1)}, \dots, Y\_{j+h-1}^{(d)}\right) - \mathbb{E}\left( f\left(Y\_j^{(1)}, \dots, Y\_{j+h-1}^{(d)}\right) \right) \right)$$

$$\stackrel{\mathcal{D}}{\rightarrow} \sum\_{p,q \in P^\*} \mathbb{A}^{(p,q)} Z\_{2,d^\*+1/2}^{(p,q)}(1),\tag{13}$$

*where:*

$$Z\_{2,d^\*+1/2}^{(p,q)}(1) = K\_{p,q}(d^\*) \int\_{\mathbb{R}^2}^{\prime\prime} \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*} \mathcal{B}\_L^{(p)}(\mathbf{d}\lambda\_1) \mathcal{B}\_L^{(q)}(\mathbf{d}\lambda\_2).$$

*The matrix K*(*d*∗) *is a normalizing constant, as can be seen in [18], Corollary 3.6. Moreover, B*˜*L*(d*λ*) *is a multivariate Hermitian–Gaussian random measure with* E(*BL*(d*λ*)*BL*(d*λ*)∗) = *L* <sup>d</sup>*λ and L as defined in* (2)*. Furthermore, C*<sup>2</sup> := <sup>1</sup> <sup>2</sup>*d*∗(4*d*∗−1) *is a normalizing constant and:*

$$\tilde{\alpha}^{(p\mathcal{A})} := \sum\_{i,k=1}^{h} \alpha\_{i,k}^{(p\mathcal{A})\_i}$$

*where α* (*p*,*q*) *<sup>i</sup>*,*<sup>k</sup>* = *<sup>α</sup>i*+(*p*−1)*h*,*k*+(*q*−1)*<sup>h</sup> for each p*, *<sup>q</sup>* ∈ *<sup>P</sup>*<sup>∗</sup> *and i*, *<sup>k</sup>* = 1, . . . , *h and:*

$$(a\_{i\lambda})\_{1 \le i,k \le dh} = \Sigma\_{d\lambda}^{-1} C \Sigma\_{d\lambda}^{-1}$$

*where C denotes the matrix of second order Hermite coefficients, given by*

$$\mathbb{C} = (\mathfrak{c}\_{i,k})\_{1 \le i,k \le dh} = \mathbb{E}\left(Y\_{1,h}(f(Y\_{1,h}) - \mathbb{E}(f(Y\_{1,h}))) | Y\_{1,h}^t\right).$$

It is possible to soften the assumptions in Theorem 1 to allow for mixed cases of shortand long-range dependence.

**Corollary 1.** *Instead of demanding in the assumptions of Theorem 1 that* (2) *holds for* - *Yj <sup>j</sup>*∈<sup>Z</sup> *with the addition that for all p* = 1, . . . , *d we have dp* ∈ 1 4 , 1 2 *, we may use the following condition. We assume that:*

$$r^{(p,q)}(k) = k^{d\_p + d\_q - 1} L\_{p,q}(k) \quad (k \to \infty),$$

*with Lp*,*q*(*k*) *as given in* (2)*, but we do no longer assume dp* ∈ 1 4 , 1 2 *for all p* = 1, ... , *d but soften the assumption to d*<sup>∗</sup> ∈ 1 4 , 1 2 *and for dp* = *d*∗*, p* = 1, ... , *d we allow for dp* ∈ (−∞, 0) ∪ 0, <sup>1</sup> 4 *. Then, the statement of Theorem 1 remains valid.*

However, with a mild technical assumption on the covariances of the one-dimensional marginal Gaussian processes that is often fulfilled in applications, there is another way of normalizing the partial sum on the right-hand side in Theorem 1, this time explicitly for the case #*P*<sup>∗</sup> <sup>=</sup> 2 and *<sup>h</sup>* <sup>∈</sup> <sup>N</sup>, such that the limit can be expressed in terms of two standard Rosenblatt random variables. This yields the possibility of further studying the dependence structure between these two random variables. In the following theorem, we assume #*P*∗ = *d* = 2 for the reader's convenience.

**Theorem 2.** *Under the same assumptions as in Theorem 1 with* #*P*<sup>∗</sup> = *d* = 2 *and d*<sup>∗</sup> ∈ 1 4 , 1 2 *and the additional condition that <sup>r</sup>*(1,1)(*l*) = *<sup>r</sup>*(2,2)(*l*)*, for <sup>l</sup>* <sup>=</sup> 0, ... , *<sup>h</sup>* <sup>−</sup> <sup>1</sup>*, and <sup>L</sup>*1,1 <sup>+</sup> *<sup>L</sup>*2,2 <sup>=</sup> *L*1,2 + *L*2,1*, it holds that:*

$$\begin{split} n^{-2d^\*} \left( \mathbb{C}\_2 \right)^{-\frac{1}{2}} \sum\_{j=1}^n \left( f \left( Y\_j^{(1)}, \dots, Y\_{j+h-1}^{(d)} \right) - \mathbb{E} f \left( Y\_j^{(1)}, \dots, Y\_{j+h-1}^{(d)} \right) \right) \\ \overset{\mathcal{D}}{\to} \left( \bar{\kappa}^{(1,1)} - \bar{\kappa}^{(1,2)} \right) \frac{L\_{2,2} - L\_{2,1} - L\_{1,2} + L\_{1,1}}{2} Z\_{2,d^\* + 1/2}^\* (1) \\ + \left( \bar{\kappa}^{(1,1)} + \bar{\kappa}^{(1,2)} \right) \frac{L\_{2,2} + L\_{2,1} + L\_{1,2} + L\_{1,1}}{2} Z\_{2,d^\* + 1/2}^{\*\*} (1) \end{split}$$

*with C*<sup>2</sup> := <sup>1</sup> <sup>2</sup>*d*∗(4*d*∗−1) *being the same normalizing factor as in Theorem 1,* (*αi*,*k*)1≤*i*,*k*≤*dh* <sup>=</sup> Σ−<sup>1</sup> *d*,*hC*Σ−<sup>1</sup> *<sup>d</sup>*,*<sup>h</sup> and <sup>C</sup>* <sup>=</sup> (*ci*,*k*)1≤*i*,*k*≤*dh* <sup>=</sup> <sup>E</sup> *<sup>Y</sup>*1,*h*(*f*(*Y*1,*h*) <sup>−</sup> <sup>E</sup>*f*(*Y*1,*h*))*Y<sup>t</sup>* 1,*h . Note that Z*∗ 2,*d*∗+1/2(1) *and Z*∗∗ 2,*d*∗+1/2(1) *are both standard Rosenblatt random variables whose covariance is given by*

$$\mathbb{C}ov\left(Z^\*\_{2,d^\*+1/2}(1), Z^{\*\ast}\_{2,d^\*+1/2}(1)\right) = \frac{\left(L\_{2,2} - L\_{1,1}\right)^2}{\left(L\_{1,1} + L\_{2,2}\right)^2 - \left(L\_{1,2} + L\_{2,1}\right)^2}.\tag{14}$$

#### **3. Ordinal Pattern Dependence**

Ordinal pattern dependence is a multivariate dependence measure that compares the co-movement of two time series based on the ordinal information. First introduced in [10] to analyze financial time series, a mathematical framework including structural breaks and limit theorems for functionals of absolutely regular processes has been built in [11]. In [19], the authors have used the so-called symbolic correlation integral in order to detect the dependence between the components of a multivariate time series. Their considerations focusing on testing independence between two time series are also based on ordinal patterns. They provide limit theorems in the i.i.d.-case and otherwise use bootstrap methods. In contrast, in the mathematical model in the present article, we focus on asymptotic distributions of an estimator of ordinal pattern dependence having a bivariate Gaussian time series in the background but allowing for several dependence structures to arise. As it will turn out in the following, this yields central but also non-central limit theorems.

We start with the definition of an ordinal pattern and the basic mathematical framework that we need to build up the ordinal model.

Let *Sh* denote the set of permutations in {0, ... , *<sup>h</sup>*}, *<sup>h</sup>* <sup>∈</sup> <sup>N</sup><sup>0</sup> that we express as (*<sup>h</sup>* <sup>+</sup> <sup>1</sup>) dimensional tuples, assuring that each tuple contains each of the numbers above exactly once. In mathematical terms, this yields:

$$S\_h = \left\{ \pi \in \mathbb{N}\_0^{h+1} : \quad 0 \le \pi\_i \le h, \text{ and } \pi\_i \ne \pi\_{k'} \text{ whenever } i \ne k, \quad i, k = 0, \dots, h \right\},$$

as can be seen in [11], Section 2.1.

t

The number of permutations in *Sh* is given by #*Sh* = (*h* + 1)!. In order to get a better intuitive understanding of the concept of ordinal patterns, we have a closer look at the following example, before turning to the formal definition.

**Example 1.** *Figure 1 provides an illustrative understanding of the extraction of an ordinal pattern from a data set. The data points of interest are colored in red and we consider a pattern of length h* = 3*, which means that we have to take n* = 4 *data points into consideration. We fix the points in time t*0*, t*1*, t*<sup>2</sup> *and t*<sup>3</sup> *and extract the data points from the time series. Then, we search for the point in time which exhibits the largest value in the resulting data and write down the corresponding time index. In this example, it was given by t* = *t*1*. We order the data points by writing the time position of the largest value as the first entry, the time position of the second largest as the second entry, etc. Hence, the absolute values are ordered from largest to smallest and the ordinal pattern* (1, 0, 3, 2) ∈ *S*<sup>3</sup> *is obtained for the considered data points.*

**Figure 1.** Example of the extraction of an ordinal pattern of a given data set.

Formally, the aforementioned procedure can be defined as follows, as can be seen in [11], Section 2.1.

**Definition 3.** *As the ordinal pattern of a vector <sup>x</sup>* <sup>=</sup> (*x*0,..., *xh*) <sup>∈</sup> <sup>R</sup>*h*+1*, we define the unique permutation π* = (*π*0,..., *πh*) ∈ *Sh:*

$$
\Pi(\mathfrak{x}) = \Pi(\mathfrak{x}\_{0\prime}, \dots, \mathfrak{x}\_{\hbar}) = (\pi\_{0\prime}, \dots, \pi\_{\hbar})\_{\prime\prime}
$$

*such that:*

$$\mathbf{x}\_{\pi\_0} \ge \dots \ge \mathbf{x}\_{\pi\_{h'}}$$

*with <sup>π</sup>i*−<sup>1</sup> < *<sup>π</sup><sup>i</sup> if xπi*−<sup>1</sup> = *<sup>x</sup>π<sup>i</sup> , i* = 1, . . . , *h.*

The last condition assures the uniqueness of *π* if there are ties in the data sets. In particular, this condition is necessary if real-world data are to be considered.

In Figure 2, all ordinal patterns of length *h* = 2 are shown. As already mentioned in the introduction, from the practical point of view, a highly desirable property of ordinal patterns is that they are not affected by monotone transformations, as can be seen in [5], p. 1783.

**Figure 2.** Ordinal patterns for *h* = 2.

Mathematically, this means that if *<sup>f</sup>* : <sup>R</sup> <sup>→</sup> <sup>R</sup> is strictly monotone, then:

$$
\Pi(\mathbf{x}\_0, \dots, \mathbf{x}\_h) = \Pi(f(\mathbf{x}\_0), \dots, f(\mathbf{x}\_h)).\tag{15}
$$

In particular, this includes linear transformations *<sup>f</sup>*(*x*) = *ax* <sup>+</sup> *<sup>b</sup>*, with *<sup>a</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> and *<sup>b</sup>* <sup>∈</sup> <sup>R</sup>.

Following [11], Section 1, the minimal requirement of the data sets we use for ordinal analysis in the time series context, i.e., for ordinal pattern probabilities as well as for ordinal pattern dependence later on, is *ordinal pattern stationarity (of order h)*. This property implies that the probability of observing a certain ordinal pattern of length *h* remains the same when shifting the moving window of length *h* through the entire time series and is not depending on the specific points in time. In the course of this work, the time series, in which the ordinal patterns occur, always have either stationary increments or are even stationary themselves. Note that both properties imply ordinal pattern stationarity. The reason why requiring stationary increments is a sufficient condition is given in the following explanation.

One fundamental property of ordinal patterns is that they are uniquely determined by the increments of the considered time series. As one can imagine in Example 1, the knowledge of the increments between the data points is sufficient to obtain the corresponding ordinal pattern. In mathematical terms, we can define another mapping Π˜ , which assigns the corresponding ordinal pattern to each vector of increments, as can be seen in [5], p. 1783.

**Definition 4.** *We define for y* <sup>=</sup> (*y*1,..., *yh*) <sup>∈</sup> <sup>R</sup>*<sup>h</sup> the mapping* <sup>Π</sup>˜ : <sup>R</sup>*<sup>h</sup>* <sup>→</sup> *Sh:*

$$\tilde{\Pi}(y\_1, \dots, y\_h) := \Pi(0, y\_1, y\_1 + y\_2, \dots, y\_1 + \dots + y\_h),$$

*such that for yi* = *xi* − *xi*−1*, i* = 1, . . . , *h, we obtain:*

$$\begin{aligned} \Pi(y\_1, \dots, y\_h) &= \Pi(0, y\_1, y\_1 + y\_2, \dots, y\_1 + \dots + y\_h) \\ &= \Pi(0, \mathbf{x}\_1 - \mathbf{x}\_0, \mathbf{x}\_2 - \mathbf{x}\_0, \dots, \mathbf{x}\_h - \mathbf{x}\_0) \\ &= \Pi(\mathbf{x}\_0, \mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_h) .\end{aligned}$$

We define the two mappings, following [5], p. 1784:

$$\begin{aligned} \mathcal{S}: \mathcal{S}\_{\mathrm{h}} &\to \mathcal{S}\_{\mathrm{h}\prime} \left( \pi\_{0\prime}, \dots, \pi\_{\mathrm{h}} \right) \to \left( \pi\_{\mathrm{h}\prime}, \dots, \pi\_{\mathrm{0}} \right), \\ \mathcal{T}: \mathcal{S}\_{\mathrm{h}} &\to \mathcal{S}\_{\mathrm{h}\prime} \left( \pi\_{0\prime}, \dots, \pi\_{\mathrm{h}} \right) \to \left( \mathrm{h} - \pi\_{0\prime}, \dots, \mathrm{h} - \pi\_{\mathrm{h}} \right). \end{aligned}$$

An illustrative understanding of these mappings is given as follows. The mapping S(*π*), which is the spatial reversion of the pattern *π*, is the reflection of *π* on a horizontal line, while T (*π*), the time reversal of *π*, is its reflection on a vertical line, as one can observe in Figure 3.

**Figure 3.** Space and time reversion of the pattern *π* = (1, 3, 2, 0).

Based on the spatial reversion, we define a possibility to divide *Sh* into two disjoint sets.

**Definition 5.** *We define S*∗ *<sup>h</sup> as a subset of Sh with the property that for each π* ∈ *Sh, either π or* S(*π*) *are contained in the set, but not both of them.*

Note that this definition does not yield the uniqueness of *S*∗ *h*.

**Example 2.** *We consider the case h* = 2 *again and we want to divide S*<sup>2</sup> *into a possible choice of S*<sup>∗</sup> 2 *and the corresponding spatial reversal. We choose S*∗ <sup>2</sup> = {(2, 1, 0),(2, 0, 1),(1, 2, 0)}*, and therefore, S*<sup>2</sup> \ *S*<sup>∗</sup> <sup>2</sup> = {(0, 1, 2),(1, 0, 2),(0, 2, 1)}*. Remark that S*<sup>∗</sup> <sup>2</sup> = {(0, 1, 2),(2, 0, 1),(1, 2, 0)} *is also a possible choice. The only condition that has to be satisfied is that if one permutation is chosen for S*∗ 2*, then its spatial reverse must not be an element of this set.*

We stick to the formal definition of ordinal pattern dependence, as it is proposed in [11], Section 2.1. The considered moving window consists of *h* + 1 data points, and hence, *h* increments. We define:

$$p := p\_{X^{(1)}, X^{(2)}} := \mathbb{P}\left(\Pi\left(X\_0^{(1)}, \dots, X\_h^{(1)}\right) = \Pi\left(X\_0^{(2)}, \dots, X\_h^{(2)}\right)\right) \tag{16}$$

and:

$$\eta := q\_{X^{(1)}, X^{(2)}} := \sum\_{\pi \in \mathcal{S}\_h} \mathbb{P}\left(\Pi\left(X\_0^{(1)}, \dots, X\_h^{(1)}\right) = \pi\right) \mathbb{P}\left(\Pi\left(X\_0^{(2)}, \dots, X\_h^{(2)}\right) = \pi\right).$$

Then, we define ordinal pattern dependence *OPD* as

$$\text{OPD} := \text{OPD}\_{\mathcal{X}^{(1)}, \mathcal{X}^{(2)}} := \frac{p - q}{1 - q}. \tag{17}$$

The parameter *q* represents the hypothetical case of independence between the two time series. In this case, *p* and *q* would obtain equal values and therefore, *OPD* would equal zero. Regarding the other extreme, the case in which both processes coincide or one is a strictly monotone increasing transform of the other one, we obtain the value 1. However, in the following, we assume *p* ∈ (0, 1) and *q* ∈ (0, 1).

Note that the definition of ordinal pattern dependence in (17) only measures positive dependence. This is no restriction in practice, because negative dependence can be investigated in an analogous way, by considering *OPDX*(1),−*X*(2). If one is interested in both types of dependence simultaneously, in [11], the authors propose to use *OPDX*(1),*X*(2) <sup>+</sup> <sup>−</sup> *OPDX*(1),−*X*(2) . To keep the notation simple, we focus on *OPD* as it is defined in (17).

+ We compare whether the ordinal patterns in *X*(1) *j <sup>j</sup>*∈<sup>Z</sup> coincide with the ones in *X*(2) *j <sup>j</sup>*∈<sup>Z</sup> . Recall that it is an essential property of ordinal patterns that they are uniquely determined by the increment process. Therefore, we have to consider the increment processes - *Yj <sup>j</sup>*∈<sup>Z</sup> <sup>=</sup> *Y*(1) *<sup>j</sup>* ,*Y*(2) *j <sup>j</sup>*∈<sup>Z</sup> as defined in (1) for *<sup>d</sup>* <sup>=</sup> 2, where *<sup>Y</sup>*(*p*) *<sup>j</sup>* <sup>=</sup> *<sup>X</sup>*(*p*) *j* − *<sup>X</sup>*(*p*) *<sup>j</sup>*−1, *<sup>p</sup>* <sup>=</sup> 1, 2. Hence, we can also express *<sup>p</sup>* and *<sup>q</sup>* (and consequently *OPD*) as a probability that only depends on the increments of the considered vectors of the time series. Recall the definition of *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> for *<sup>d</sup>* <sup>=</sup> 2, given by

$$\mathcal{Y}\_{j,\mathfrak{h}} = \left(\mathcal{Y}\_j^{(1)}, \dots, \mathcal{Y}\_{j+\mathfrak{h}-1}^{(1)}, \mathcal{Y}\_j^{(2)}, \dots, \mathcal{Y}\_{j+\mathfrak{h}-1}^{(2)}\right)^t,$$

such that *Yj*,*<sup>h</sup>* ∼ N (0, Σ2,*h*) with Σ2,*<sup>h</sup>* as given in (6).

In the course of this article, we focus on the estimation of *p*. For a detailed investigation of the limit theorems for estimators of *OPD*, we refer to [18]. We define the estimator of *p*, the probability of coincident patterns in both time series in a moving window of fixed length, by

$$\begin{split} \hat{\mathcal{P}}\_{\text{lt}} &= \frac{1}{n-h} \sum\_{j=0}^{n-h-1} \mathbf{1} \left\{ \boldsymbol{\Pi} \left( \boldsymbol{\chi}\_{j}^{(1)}, \dots, \boldsymbol{\chi}\_{j+h}^{(1)} \right) - \boldsymbol{\Pi} \left( \boldsymbol{\chi}\_{j}^{(2)}, \dots, \boldsymbol{\chi}\_{j+h}^{(2)} \right) \right\} \\ &= \frac{1}{n-h} \sum\_{j=1}^{n-h} \mathbf{1} \left\{ \boldsymbol{\Pi} \left( \boldsymbol{\chi}\_{j}^{(1)}, \dots, \boldsymbol{\chi}\_{j+h-1}^{(1)} \right) - \boldsymbol{\Pi} \left( \boldsymbol{\chi}\_{j}^{(2)}, \dots, \boldsymbol{\chi}\_{j+h-1}^{(2)} \right) \right\} \end{split}$$

where:

$$\begin{aligned} \Pi(\mathbf{Y}\_1, \dots, \mathbf{Y}\_h) &:= \Pi(0, \mathbf{Y}\_1, \mathbf{Y}\_1 + \mathbf{Y}\_2, \dots, \mathbf{Y}\_1 + \dots + \mathbf{Y}\_h) \\ &= \Pi(0, \mathbf{X}\_1 - \mathbf{X}\_{0\prime}, \dots, \mathbf{X}\_h - \mathbf{X}\_0) \\ &= \Pi(\mathbf{X}\_{0\prime} \mathbf{X}\_1, \dots, \mathbf{X}\_h). \end{aligned}$$

Figure 4 illustrates the way ordinal pattern dependence is estimated by *p*ˆ*n*. The patterns of interest that are compared in each moving window are colored in red.

**Figure 4.** Illustration of estimation of ordinal pattern dependence.

Having emphasized the crucial importance of the increments, we define the following conditions on the increment process - *Yj <sup>j</sup>*∈Z: let - *Yj <sup>j</sup>*∈<sup>Z</sup> be a bivariate, stationary Gaussian process with *<sup>Y</sup>*(*p*) *<sup>j</sup>* ∼ N (0, 1), *p* = 1, 2:


$$r^{(p,q)}(k) = k^{d\_p + d\_q - 1} L\_{p,q}(k) \quad (k \to \infty),$$

with *Lp*,*q*(*k*) <sup>→</sup> *Lp*,*<sup>q</sup>* and *Lp*,*<sup>q</sup>* <sup>∈</sup> <sup>R</sup> holds.

Furthermore, in both cases, it holds that *r*(*p*,*q*)(*k*) <sup>&</sup>lt; 1 for *<sup>p</sup>*, *<sup>q</sup>* <sup>=</sup> 1, 2 and *<sup>k</sup>* <sup>≥</sup> 1 to exclude ties.

We begin with the investigation of the asymptotics of *p*ˆ*n*. First, we calculate the Hermite rank of *p*ˆ*n*, since the Hermite rank determines for which ranges of *d*∗ the estimator *p*ˆ*<sup>n</sup>* is still long-range dependent. Depending on this range, different limit theorems may hold.

**Lemma 1.** *The Hermite rank of f*(*Yj*,*h*) = **1**. Π˜ *<sup>Y</sup>*(1) *j*+1,...,*Y*(1) *j*+*h* =Π˜ *<sup>Y</sup>*(2) *j*+1,...,*Y*(2) *j*+*h* / *with respect to* Σ2,*<sup>h</sup> is equal to* 2*.*

**Proof.** Following [20], Lemma 5.4 it is sufficient to show the following two properties:

(i) *m*(*f* , Σ2,*h*) ≥ 2, (ii) *m*(*f* , *I*2,*h*) ≤ 2.

Note that the conclusion is not trivial, because *m*(*f* , Σ2,*h*) = *m*(*f* , *I*2,*h*) in general, as can be seen in [15], Lemma 3.7. Lemma 5.4 in [20] can be applied due to the following reasoning. Ordinal patterns are not affected by scaling, therefore, the technical condition that Σ−<sup>1</sup> 2,*<sup>h</sup>* − *I*2,*<sup>h</sup>* is positive semidefinite is fulfilled in our case. We can scale the standard deviation of the random vector *Yj*,*<sup>h</sup>* by any positive real number *<sup>σ</sup>* <sup>&</sup>gt; 0 since for all *<sup>j</sup>* <sup>∈</sup> <sup>Z</sup> we have:

$$\begin{aligned} \left\{ \Pi \left( \chi\_{\stackrel{j}{j}}^{(1)}, \dots, \chi\_{\stackrel{j}{j+h-1}}^{(1)} \right) = \Pi \left( \chi\_{\stackrel{j}{j}}^{(2)}, \dots, \chi\_{\stackrel{j+h-1}{j+h-1}}^{(2)} \right) \right\} \\ = \left\{ \Pi \left( \sigma \chi\_{\stackrel{j}{j}}^{(1)}, \dots, \sigma \chi\_{\stackrel{j+h-1}{j+h-1}}^{(1)} \right) = \Pi \left( \sigma \chi\_{\stackrel{j}{j}}^{(2)}, \dots, \sigma \chi\_{\stackrel{j+h-1}{j+h-1}}^{(2)} \right) \right\}. \end{aligned}$$

To show property (*i*), we need to consider a multivariate random vector:

$$\boldsymbol{\Upsilon}\_{1,h} := \left( \boldsymbol{\Upsilon}\_1^{(1)}, \dots, \boldsymbol{\Upsilon}\_h^{(1)}, \boldsymbol{\Upsilon}\_1^{(2)}, \dots, \boldsymbol{\Upsilon}\_h^{(2)} \right)^t$$

with covariance matrix Σ2,*h*. We fix *i* = 1, ... , 2*h*. We divide the set *Sh* into disjoint sets, namely into *S*∗ *<sup>h</sup>*, as defined in Definition 5 and the complimentary set *Sh* \ *S*<sup>∗</sup> *<sup>h</sup>*. Note that:

$$-\Upsilon\_{j,k} \stackrel{p}{=} \Upsilon\_{j,k}$$

holds. This implies:

E *Y*(*i*) *<sup>j</sup>*,*<sup>h</sup>* **1**. Π˜ *<sup>Y</sup>*(1) <sup>1</sup> ,...,*Y*(1) *h* =Π˜ *<sup>Y</sup>*(2) <sup>1</sup> ,...,*Y*(2) *h* =*π* / <sup>=</sup> <sup>−</sup><sup>E</sup> *Y*(*i*) *<sup>j</sup>*,*<sup>h</sup>* **1**. Π˜ *<sup>Y</sup>*(1) <sup>1</sup> ,...,*Y*(1) *h* =Π˜ *<sup>Y</sup>*(2) <sup>1</sup> ,...,*Y*(2) *h* =S(*π*) / 

for *π* ∈ *Sh*. Hence, we arrive at:

$$\begin{split} \mathbb{E}\left(\boldsymbol{Y}\_{j,h}^{(i)}f(\boldsymbol{Y}\_{j,h})\right) &= \mathbb{E}\left(\boldsymbol{Y}\_{j,h}^{(i)}\mathbf{1}\left\{\boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(1)},\ldots,\boldsymbol{Y}\_{h}^{(1)}\right) - \boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(2)},\ldots,\boldsymbol{Y}\_{h}^{(2)}\right)\right\}\right) \\ &= \sum\_{\pi \in S\_{h}} \mathbb{E}\left(\boldsymbol{Y}\_{j,h}^{(i)}\mathbf{1}\left\{\boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(1)},\ldots,\boldsymbol{Y}\_{h}^{(1)}\right) - \boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(2)},\ldots,\boldsymbol{Y}\_{h}^{(2)}\right) - \pi\right\}\right) \\ &= \sum\_{\pi \in S\_{h}^{\*}} \mathbb{E}\left(\boldsymbol{Y}\_{j,h}^{(i)}\mathbf{1}\left\{\boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(1)},\ldots,\boldsymbol{Y}\_{h}^{(1)}\right) - \boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(2)},\ldots,\boldsymbol{Y}\_{h}^{(2)}\right) - \pi\right\}\right) \\ & \quad - \sum\_{\pi \in S\_{h} \cup S\_{h}^{\*}} \mathbb{E}\left(\boldsymbol{Y}\_{j,h}^{(i)}\mathbf{1}\left\{\boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(1)},\ldots,\boldsymbol{Y}\_{h}^{(1)}\right) - \boldsymbol{\Pi}\left(\boldsymbol{Y}\_{1}^{(2)},\ldots,\boldsymbol{Y}\_{h}^{(2)}\right) - \boldsymbol{S}(\boldsymbol{\pi})\right\}\right) \\ &= 0 \end{split}$$

for *i* = 1, . . . , 2*h*.

Consequently, *m*(*f* , Σ2,*h*) ≥ 2. In order to prove (*ii*), we consider:

$$\mathcal{U}\_{1,h} := \left(\mathcal{U}\_1^{(1)}, \dots, \mathcal{U}\_h^{(1)}, \mathcal{U}\_1^{(2)}, \dots, \mathcal{U}\_h^{(2)}\right)^t$$

to be a random vector with independent N (0, 1) distributed entries. For *i* = 1, ... , *h* and *k* = *h* + 1, . . . , 2*h* such that *k* − *h* = *i*, we obtain:

$$\begin{split} \mathbb{E}\left(\mathsf{U}\_{1,h}^{(1)}\mathcal{U}\_{1,h}^{(1)}f(\mathcal{U}\_{1,h})\right) &= \mathbb{E}\left(\mathsf{U}\_{i}^{(1)}\mathcal{U}\_{k-h}^{(2)}\mathbf{1}\_{\left\{\mathbb{1}\left(\mathcal{U}\_{1}^{(1)},\dots,\mathcal{U}\_{h}^{(1)}\right) - \text{fl}\left(\mathsf{U}\_{1}^{(2)},\dots,\mathsf{U}\_{h}^{(2)}\right)\right\}}\right) \\ &= \sum\_{\pi \in S\_{h}} \mathbb{E}\left(\mathsf{U}\_{i}^{(1)}\mathcal{U}\_{i}^{(2)}\mathbf{1}\_{\left\{\mathbb{1}\left(\mathcal{U}\_{1}^{(1)},\dots,\mathcal{U}\_{h}^{(1)}\right) - \text{fl}\left(\mathsf{U}\_{1}^{(2)},\dots,\mathsf{U}\_{h}^{(2)}\right) - \pi\right\}}\right) \\ &= \sum\_{\pi \in S\_{h}} \left(\mathbb{E}\left(\mathsf{U}\_{i}^{(1)}\mathbf{1}\_{\left\{\mathbb{1}\left(\mathcal{U}\_{1}^{(1)},\dots,\mathcal{U}\_{h}^{(1)}\right) - \pi\right\}}\right)\right)^{2} \\ &\neq 0, \end{split}$$

since E *U*(1) *<sup>i</sup>* **1**. Π˜ *<sup>U</sup>*(1) <sup>1</sup> ,...,*U*(1) *h* =*π* / = 0 for all *π* ∈ *Sh*. This was shown in the proof of Lemma 3.4 in [20].

All in all, we derive *m*(*f* , Σ2,*h*) = 2 and hence, have proven the lemma.

The case *m*(*f* , Σ2,*h*) = 2 exhibits the property that the standard range of the long-range dependence parameter *d*<sup>∗</sup> ∈ 0, <sup>1</sup> 2 has to be divided into two different sets. If *d*<sup>∗</sup> ∈ 1 4 , 1 2 , the transformed process *f Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> is still long-range dependent, as can be seen in [16], Table 5.1. If *d*<sup>∗</sup> ∈ 0, <sup>1</sup> 4 , the transformed process is short-range dependent, which means by definition that the autocorrelations of the transformed process are summable, as can be seen in [13], Remark 2.3. Therefore, we have two different asymptotic distributions that have to be considered for the estimator *p*ˆ*<sup>n</sup>* of coincident patterns.

#### *3.1. Limit Theorem for the Estimator of p in Case of Long-Range Dependence*

First, we restrict ourselves to the case that at least one of the two parameters *d*<sup>1</sup> and *d*<sup>2</sup> is in <sup>1</sup> 4 , 1 2 . This assures *d*<sup>∗</sup> ∈ 1 4 , 1 2 . We explicitly include mixing cases where the process corresponding to min{*d*1, *d*2} is allowed to be long-range as well as short-range dependent.

Note that this setting includes the pure long-range dependence case, which means that for *p* = 1, 2, we have *dp* ∈ 1 4 , 1 2 , or even *d*<sup>1</sup> = *d*<sup>2</sup> = *d*∗. However, in general, the assumptions are lower, such that we only require *dp* ∈ 1 4 , 1 2 for either *p* = 1 or *p* = 2 and the other parameter is also allowed to be in (−∞, 0) or 0, <sup>1</sup> 4 .

We can, therefore, apply the results of Corollary 1 and obtain the following asymptotic distribution for *p*ˆ*n*:

**Theorem 3.** *Under the assumption in (L), we obtain:*

$$n^{1-2d^\*} (C\_2)^{-\frac{1}{2}} (\not p\_n - p) \stackrel{\mathcal{D}}{\longrightarrow} \sum\_{p,q \in P^\*} \overline{a}^{(p,q)} Z\_{2,d^\*+1/2}^{(p,q)}(1) \tag{18}$$

*with <sup>Z</sup>*(*p*,*q*) 2,*d*∗+1/2(1) *as given in Theorem <sup>1</sup> for <sup>p</sup>*, *<sup>q</sup>* <sup>∈</sup> *<sup>P</sup>*<sup>∗</sup> *and <sup>C</sup>*<sup>2</sup> :<sup>=</sup> <sup>1</sup> <sup>2</sup>*d*∗(4*d*∗−1) *being a normalizing constant. We have:*

$$\tilde{\alpha}^{(p,q)} := \sum\_{i,k=1}^{h} \alpha\_{i,k}^{(p,q)}, \text{ where } \alpha\_{i,k}^{(p,q)} = \alpha\_{i+(p-1)h,k+(q-1)h}.$$

*for each p*, *<sup>q</sup>* <sup>∈</sup> *<sup>P</sup>*<sup>∗</sup> *and i*, *<sup>k</sup>* <sup>=</sup> 1, . . . , *h and* (*αi*,*k*)1≤*i*,*k*≤*dh* <sup>=</sup> <sup>Σ</sup>−<sup>1</sup> 2,*hC*Σ−<sup>1</sup> 2,*<sup>h</sup> , where the variable:*

$$\mathcal{C} = (\mathfrak{c}\_{i,k})\_{1 \le i,k \le 2h} = \mathbb{E}\left(\mathbf{Y}\_{1,h}\left(\mathbf{1}\_{\left\{\mathrm{\hat{\mathbf{I}}\left(\mathbf{Y}\_1^{(1)},\dots,\mathbf{Y}\_h^{(1)}\right) - \mathrm{\hat{\mathbf{I}}\left(\mathbf{Y}\_1^{(2)},\dots,\mathbf{Y}\_h^{(2)}\right)}\right\}} - p\right)\mathbf{Y}\_{1,h}^t\right)$$

*denotes the matrix of second order Hermite coefficients.*

**Proof.** The proof of this theorem is an immediate application of the Corollary 1 and Lemma 1. Note that for *p*ˆ*<sup>n</sup>* it holds that it is square integrable with respect to *Yj*,*<sup>h</sup>* and that the set of discontinuity points is a Null set with respect to the 2*h*-dimensional Lebesgue measure. This is shown in [18], Equation (4.5).

Following Theorem 2, we are also able to express the limit distribution above in terms of two standard Rosenblatt random variables by modifying the weighting factors in the limit distribution. Note that this requires slightly stronger assumptions as in Theorem 1.

**Theorem 4.** *Let (L) hold with d*<sup>1</sup> = *d*2*. Additionally, we assume that r*(1,1)(*l*) = *r*(2,2)(*l*)*, for l* = 0, . . . , *h* − 1*, and L*1,1 + *L*2,2 = *L*1,2 + *L*2,1*. Then, we obtain:*

$$\begin{split} n^{1-2d^\*} (\mathsf{C}\_2)^{-\frac{1}{2}} (\mathfrak{P}\_n - p) \xrightarrow{\mathcal{D}} \left( \mathring{\mathfrak{a}}^{(1,1)} - \mathring{\mathfrak{a}}^{(1,2)} \right) \frac{L\_{2,2} - L\_{2,1} - L\_{1,2} + L\_{1,1}}{2} \mathcal{Z}^\*\_{2,d^\* + 1/2} (1) \\ + \left( \mathring{\mathfrak{a}}^{(1,1)} + \mathring{\mathfrak{a}}^{(1,2)} \right) \frac{L\_{2,2} + L\_{2,1} + L\_{1,2} + L\_{1,1}}{2} \mathcal{Z}^{\*\*}\_{2,d^\* + 1/2} (1), \end{split}$$

*with C*<sup>2</sup> *and α*˜(*p*,*q*) *as given in Theorem 3. Note that Z*<sup>∗</sup> 2,*d*∗+1/2(1) *and Z*∗∗ 2,*d*∗+1/2(1) *are both standard Rosenblatt random variables, whose covariance is given by*

$$\mathbb{C}ov\left(Z^\*\_{2,d^\*+1/2}(1), Z^{\*\ast}\_{2,d^\*+1/2}(1)\right) = \frac{\left(L\_{2,2} - L\_{1,1}\right)^2}{\left(L\_{1,1} + L\_{2,2}\right)^2 - \left(L\_{1,2} + L\_{2,1}\right)^2}.\tag{19}$$

**Remark 1.** *Following [18], Corollary 3.14, if additionally r*(1,1)(*k*) = *r*(2,2)(*k*) *and r*(1,2)(*k*) = *<sup>r</sup>*(2,1)(*k*) *is fulfilled for all <sup>k</sup>* <sup>∈</sup> <sup>Z</sup>*, then the two limit random variables following a standard Rosenblatt distribution in Theorem 4 are independent. Note that due to the considerations in [21], Equation (10), we know that the distribution of the sum of two independent standard Rosenblatt random variables is not standard Rosenblatt. However, this yields a computational benefit, as it is possible to efficiently simulate the standard Rosenblatt distribution, for details, as can be seen in [21].*

We turn to an example that deals with the asymptotic variance of the estimator of *p* in Theorem 3 in the case *h* = 1.

**Example 3.** *We focus on the case h* = 1 *and consider the underlying process* - *Yj*,1 *<sup>j</sup>*∈<sup>Z</sup> <sup>=</sup> *Y*(1) *<sup>j</sup>* ,*Y*(2) *j <sup>j</sup>*∈<sup>Z</sup> *. It is possible to determine the asymptotic variance depending on the correlation r*(1,2)(0) *between these two increment variables.*

*We start with the calculation of the second order Hermite coefficients in the case π* = (1, 0)*. This corresponds to the event* . *Y*(1) *<sup>j</sup>* <sup>≥</sup> 0,*Y*(2) *<sup>j</sup>* ≥ 0 / *, which yields:*

$$\mathcal{L}\_{1,1}^{\pi,2} = \mathbb{E}\left(\left(\left(\boldsymbol{Y}\_{\boldsymbol{\jmath}}^{(1)}\right)^2 - 1\right)\mathbf{1}\_{\left\{\boldsymbol{Y}\_{\boldsymbol{\jmath}}^{(1)} \ge 0, \boldsymbol{Y}\_{\boldsymbol{\jmath}}^{(2)} \ge 0\right\}}\right)^2$$

*and:*

$$c\_{1,2}^{\pi,2} = \mathbb{E}\left( \left( \mathbf{Y}\_j^{(1)} \mathbf{Y}\_j^{(2)} \right) \mathbf{1}\_{\left\{ \mathbf{Y}\_j^{(1)} \succeq 0, \mathbf{Y}\_j^{(2)} \succeq 0 \right\}} \right).$$

*Due to <sup>r</sup>*(1,2)(0) = *<sup>r</sup>*(2,1)(0)*, we have Y*(1) *<sup>j</sup>* ,*Y*(2) *j* <sup>D</sup> = *Y*(2) *<sup>j</sup>* ,*Y*(1) *j and therefore, c π*,2 1,1 = *c π*,2 2,2 *. We identify the second order Hermite coefficients as the ones already calculated in [20], Example 3.13, although we are considering two consecutive increments of a univariate Gaussian process there. However, since the corresponding values are only determined by the correlation between the Gaussian variables, we can simply replace the autocorrelation at lag* 1 *by the cross-correlation at lag* 0*. Hence, we obtain:*

$$\begin{aligned} c\_{1,1}^{\pi,2} &= \varphi^2(0)r^{(1,2)}(0)\sqrt{1 - \left(r^{(1,2)}(0)\right)^2}, \\ c\_{1,2}^{\pi,2} &= \varphi^2(0)\sqrt{1 - \left(r^{(1,2)}(0)\right)^2}. \end{aligned}$$

*Recall that the inverse* Σ−<sup>1</sup> 2,1 <sup>=</sup> - *gi*,*<sup>j</sup> <sup>i</sup>*,*j*=1,2 *of the correlation matrix of Y*(1) *<sup>j</sup>* ,*Y*(2) *j is given by*

,

.

$$
\Sigma\_{2,1}^{-1} = \frac{1}{1 - \left(r^{(1,2)}(0)\right)^2} \begin{pmatrix} 1 & -r^{(1,2)}(0) \\ -r^{(1,2)}(0) & 1 \end{pmatrix},
$$

*By using the formula for α*˜(*p*,*q*) *obtained in [18], Equation (4.23), we derive:*

$$\begin{aligned} \mathfrak{a}\_{\pi,2}^{(1,1)} &= \mathfrak{a}\_{1,1}^{\pi,2} = \left(\mathfrak{g}\_{1,1}^{2} + \mathfrak{g}\_{1,2}^{2}\right) \mathfrak{c}\_{1,1}^{\pi,2} + 2 \mathfrak{g}\_{1,1} \mathfrak{g}\_{1,2} \mathfrak{c}\_{1,2}^{\pi,2}, \\ \mathfrak{a}\_{\pi,2}^{(1,2)} &= \mathfrak{a}\_{1,2}^{\pi,2} = \left(\mathfrak{g}\_{1,1}^{2} + \mathfrak{g}\_{1,2}^{2}\right) \mathfrak{c}\_{1,2}^{\pi,2} + 2 \mathfrak{g}\_{1,1} \mathfrak{g}\_{1,2} \mathfrak{c}\_{1,1}^{\pi,2}. \end{aligned}$$

*Plugging the second order Hermite coefficients and the entries of the inverse of the covariance matrix depending on r*(1,2)(0) *into the formulas, we arrive at:*

$$
\bar{\mathfrak{a}}\_{
\pi,2}^{(1,1)} = \frac{-q^2(0)r^{(1,2)}(0)}{\left(1 - \left(r^{(1,2)}(0)\right)^2\right)^{1/2}}
$$

*and:*

$$
\bar{\alpha}\_{\pi,2}^{(1,2)} = \frac{q^2(0)}{\left(1 - \left(r^{(1,2)}(0)\right)^2\right)^{1/2}}.
$$

*Therefore, in the case h* = 1*, we obtain the following factors in the limit variance in Theorem 3:*

$$
\begin{aligned}
\tilde{\mathfrak{a}}^{(1,1)} &= \tilde{\mathfrak{a}}^{(2,2)} = \frac{-2\varrho^2(0)r^{(1,2)}(0)}{\left(1 - \left(r^{(1,2)}(0)\right)^2\right)^{1/2}}, \\
\tilde{\mathfrak{a}}^{(1,2)} &= \tilde{\mathfrak{a}}^{(2,1)} = \frac{2\varrho^2(0)}{\left(1 - \left(r^{(1,2)}(0)\right)^2\right)^{1/2}}.
\end{aligned}
$$

**Remark 2.** *It is not possible to analytically determine the limit variance for h* = 2*, as this includes orthant probabilities of a four-dimensional Gaussian distribution. Following [22], no closed formulas are available for these probabilities. However, there are fast algorithms at hand that calculate the limit variance efficiently. It is possible to take advantage of the symmetry properties of the multivariate Gaussian distribution to keep the computational cost of these algorithms low. For detail, as can be seen in [18], Section 4.3.1.*

#### *3.2. Limit Theorem for the Estimator of p in Case of Short-Range Dependence*

In this section, we focus on the case of *d*<sup>∗</sup> ∈ (−∞, 0) ∪ 0, <sup>1</sup> 4 . If *d*<sup>∗</sup> ∈ 0, <sup>1</sup> 4 , we are still dealing with a long-range dependent multivariate Gaussian process *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> . However, the transformed process *p*ˆ*<sup>n</sup>* − *p* is no longer long-range dependent, since we are considering a function with Hermite rank 2, see also [16], Table 5.1. Otherwise, if *<sup>d</sup>*<sup>∗</sup> <sup>∈</sup> (−∞, 0), the process *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> itself is already short-range dependent, since the crosscorrelations are summable. Therefore, we obtain the following central limit theorem by applying Theorem 4 in [14].

**Theorem 5.** *Under the assumptions in (S), we obtain:*

$$n^{\frac{1}{2}}(\not p\_n - p) \xrightarrow{\mathcal{D}} \mathcal{N}\left(0, \sigma^2\right),$$

*with:*

$$\begin{split} \sigma^{2} = \sum\_{k=-\infty}^{\infty} \mathbb{E}\left[ \left( \mathbf{1}\_{\left\{ \Pi\left( Y\_{1}^{(1)}, \dots, Y\_{h}^{(1)} \right) - \Pi\left( Y\_{1}^{(2)}, \dots, Y\_{h}^{(2)} \right) \right) \right\} - p \right) \\ \times \left( \mathbf{1}\_{\left\{ \Pi\left( Y\_{1+k}^{(1)}, \dots, Y\_{h+k}^{(1)} \right) - \Pi\left( Y\_{1+k}^{(2)}, \dots, Y\_{h+k}^{(2)} \right) \right\} - p} \right) \right]. \end{split}$$

We close this section with a brief retrospect of the results obtained. We established limit theorems for the estimator of *p* as probability of coincident pattern in both time series and hence, on the most important parameter in the context of ordinal pattern dependence. The long-range dependent case as well as the mixed case of short- and long-range dependence was considered. Finally, we provided a central limit theorem for a multivariate Gaussian

time series that is short-range dependent if transformed by *p*ˆ*n*. In the subsequent section, we provide a simulation study that illustrates our theoretical findings. In doing so, we shed light on the Rosenblatt distribution and the distribution of the sum of Rosenblatt distributed random variables.

#### **4. Simulation Study**

We begin with the generation of a bivariate long-range dependent fractional Gaussian noise series *Y*(1) *<sup>j</sup>* ,*Y*(2) *j j*=1,...,*n* .

First, we simulate two independent fractional Gaussian noise processes *U*(1) *j* 

*j*=1,...,*n* and *U*(2) *j j*=1,...,*n* derived by the R-package "longmemo", for a fixed parameter *H* ∈ 1 <sup>2</sup> , 1 in both time series. For the reader's convenience, we denote the long-range dependence parameter *d* by *H* = *d* + <sup>1</sup> <sup>2</sup> as it is common, when dealing with fractional Gaussian noise and fractional Brownian motion. We refer to *H* as *Hurst parameter*, tracing back to the work of [23]. For *H* = 0.7 and *H* = 0.8 we generate *n* = 106 samples, for *H* = 0.9, we choose *<sup>n</sup>* <sup>=</sup> <sup>2</sup> · <sup>10</sup>6. We denote the correlation function of univariate fractional Gaussian noise by *r* (1,1) *<sup>H</sup>* (*k*), *<sup>k</sup>* <sup>≥</sup> 0. Then, we obtain *Y*(1) *<sup>j</sup>* ,*Y*(2) *j j* for *j* = 1, . . . , *n*:

$$\begin{aligned} \mathcal{Y}\_{j}^{(1)} &= \mathcal{U}\_{j}^{(1)},\\ \mathcal{Y}\_{j}^{(2)} &= \psi \mathcal{U}\_{j}^{(1)} + \phi \mathcal{U}\_{j}^{(2)}, \end{aligned} \tag{20}$$

for *<sup>ψ</sup>*, *<sup>φ</sup>* <sup>∈</sup> <sup>R</sup>.

Note that this yields the following properties for the cross-correlations of the two processes for *k* ≥ 0:

$$\begin{aligned} r\_H^{(1,2)}(k) &= \mathbb{E} \left( Y\_j^{(1)} Y\_{j+k}^{(2)} \right) = \psi r\_H^{(1,1)}(k) \\ r\_H^{(2,1)}(k) &= r^{(1,2)}(-k) = \psi r\_H^{(1,1)}(k) \\ r\_H^{(2,2)}(k) &= \mathbb{E} \left( Y\_j^{(2)} Y\_{j+k}^{(2)} \right) = \left( \psi^2 + \phi^2 \right) r\_H^{(1,1)}(k) .\end{aligned}$$

We use *ψ* = 0.6 and *φ* = 0.8 to obtain unit variance in the second process.

Note that we chose the same Hurst parameter in both processes to get a better simulation result. The simulations of the processes *Y*(1) *j <sup>j</sup>*∈<sup>Z</sup> and *Y*(2) *j <sup>j</sup>*∈<sup>Z</sup> are visualized in Figure 5. On the left-hand side, the different fractional Gaussian noises depending on the Hurst parameter *H* are displayed. They represent the stationary long-range dependent Gaussian *increment* processes we need in the view of the limit theorems we derived in Section 3. The processes in which we are comparing the coincident ordinal patterns, namely *X*(1) *j <sup>j</sup>*∈<sup>Z</sup> and *X*(2) *j <sup>j</sup>*∈<sup>Z</sup> , are shown on the right-hand side in Figure 5. The long-range dependent behavior of the increment processes is very illustrative in these processes: roughly speaking, they become smoother the larger the Hurst parameter gets.

**Figure 5.** Plots of 500 data points of one path of two dependent fractional Gaussian noise processes (**left**) and the paths of the corresponding fractional Brownian motions (**right**) for different Hurst parameters: *H* = 0.7 (**top**), *H* = 0.8 (**middle**), *H* = 0.9 (**bottom**).

We turn to the simulation results for the asymptotic distribution of the estimator *p*ˆ*n*. The first limit theorem is given in Theorem 3 for *H* = 0.8 and *H* = 0.9. In the case of *H* = 0.7, a different limit theorem holds, see Theorem 5. Therefore, we turn to the simulation results of the asymptotic distribution of the estimator *p*ˆ*<sup>n</sup>* of *p*, as shown in Figure 6 for pattern length *h* = 2. The asymptotic normality in case *H* = 0.7 can be clearly observed. We turn to the interpretation of the simulation results of the distribution of *p*ˆ*<sup>n</sup>* − *p* for *H* = 0.8 and *H* = 0.9 as the weighted sum of the sample (cross-)correlations: we observe in the Q–Q plot for *H* = 0.8 that the samples in the upper and lower tail deviate from the reference line. For *H* = 0.9, a similar behavior in the Q–Q plot is observed.

**Figure 6.** Histogram, kernel density estimation and Q–Q plot with respect to the normal distribution (*H* = 0.7) or to the Rosenblatt distribution of *p*ˆ*<sup>n</sup>* − *p* with *h* = 2 for different Hurst parameters: *H* = 0.7 (**top**); *H* = 0.8 (**middle**); *H* = 0.9 (**bottom**).

We want to verify the result in Theorem 4 that it is possible, by a different weighting, to express the limit distribution of *p*ˆ*<sup>n</sup>* − *p* as the distribution of the sum of two independent standard Rosenblatt random variables. The simulated convergence result is provided in Figure 7. We observed the standard Rosenblatt distribution.

**Figure 7.** Histogram, kernel density estimation and Q–Q plot with respect to the Rosenblatt distribution of <sup>1</sup> *<sup>n</sup>* <sup>∑</sup>*<sup>n</sup> <sup>j</sup>*=<sup>1</sup> *H*<sup>2</sup> *Y*∗ *j* for different Hurst parameters: *H* = 0.8 (**top**); *H* = 0.9 (**bottom**).

#### **5. Conclusions and Outlook**

We considered limit theorems in the context of the estimation of ordinal pattern dependence in the long-range dependence setting. Pure long-range dependence, as well as mixed cases of short- and long-range dependence, were considered alongside the transformed short-range dependent case. Therefore, we complemented the asymptotic results in [11]. Hence, we made ordinal pattern dependence applicable for long-range dependent data sets as they arise in the context of neurology, as can be seen in [24] or artificial intelligence, as can be seen in [25]. As these kinds of data were already investigated using ordinal patterns, as can be seen, for example, in [26], this emphasizes the large practical impact of the ordinal approach in analyzing the dependence structure multivariate time series. This yields various research opportunities in these fields in the future.

Our results rely on the assumption of Gaussianity of the considered multivariate time series. If we focus on comparing the coincident ordinal patterns in a stationary longrange dependent bivariate time series, we highly benefit from the property of ordinal patterns not being affected by monotone transformations. It is possible to transform the data set to the Gaussian framework without losing the necessary ordinal information. In applications, this property is highly desirable. If we consider the more general setting, that is, stationary increments, the mathematical theory in the background gets a lot more complex leading to the limitations of our results. A crucial argument used in the proofs of the results in Section 2 is given in the *Reduction Theorem*, originally proven in Theorem 4.1 in [27] in the univariate case and extended to the multivariate setting in Theorem 6 in [14]. For further details, we refer the reader to the Appendix A. However, this result only holds in the Gaussian case. Limit theorems for the sample cross-correlation process of multivariate linear long-range dependent processes with Hermite rank 2 have recently been proven in Theorem 4 in [28]. This is possibly an interesting starting point to adapt the proofs in the Appendix A to this larger class of processes without requiring Gaussianity. Considering the property of having a discrete bivariate time series in the background, an interesting extension is given in time continuous processes and the associated techniques of discretization to still regard the ordinal perspective. To think even further beyond our scope, a generalization to categorical data is conceivable and yields an interesting open research opportunity.

**Author Contributions:** Conceptualization, I.N. and A.S.; methodology and mathematical theory, I.N.; simulations, I.N.; validation, I.N. and A.S.; writing—original draft preparation, I.N.; writing—review and editing, A.S.; funding acquisition, A.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the German Research Foundation (DFG) grant number SCHN 1231/3-2.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author, Ines Nüßgen, upon reasonable request.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Technical Appendix**

All proofs in this appendix were taken from [18], Chapter 3.

#### *Appendix A.1. Preliminary Results*

Before turning to limit theorems, we introduce a possibility to decompose the *d*dimensional Gaussian process - *Yj <sup>j</sup>*∈<sup>Z</sup> using the Cholesky decomposition, as can be seen in [29]. Based on the definition of the multivariate normal distribution, as can be seen in [30], Definition 1.6.1, we find an upper triangular matrix *A*˜, such that *A*˜ *A*˜*<sup>t</sup>* = Σ*d*. Then, it holds that:

$$
\gamma\_{\dot{j}} \stackrel{\mathcal{D}}{=} \bar{A} \mathcal{U}\_{\dot{j}}^\*,\tag{A1}
$$

where *U*∗ *<sup>j</sup>* is a *d*-dimensional Gaussian process where each *U*<sup>∗</sup> *<sup>j</sup>* has independent and identically <sup>N</sup> (0, 1) distributed entries. We want to assure that *U*∗ *j <sup>j</sup>*∈<sup>Z</sup> preserves the long-range dependent structure of - *Yj <sup>j</sup>*∈Z. Since we know from (2) that:

$$\mathbb{E}\left(\mathbf{Y}\_{\vec{\boldsymbol{\gamma}}}\mathbf{Y}\_{\vec{\boldsymbol{\gamma}}+\boldsymbol{k}}\right) = \Gamma\_{\mathbf{Y}}(\boldsymbol{k}) \simeq k^{D-\frac{1}{2}I\_{d}}Lk^{D-\frac{1}{2}I\_{d}} \quad (\boldsymbol{k}\to\infty),$$

the process *U*∗ *j* has to fulfill:

$$\mathbb{E}\left(\mathcal{U}\_{j}^{\*}\mathcal{U}\_{j+k}^{\*}\right) = \Gamma\_{\mathcal{U}^{\*}}(k) \simeq k^{D-\frac{1}{2}I\_{d}}L\_{\mathcal{U}}k^{D-\frac{1}{2}I\_{d}} \quad (k \to \infty),\tag{A2}$$

with *L* = *AL*˜ *<sup>U</sup>*<sup>∗</sup> *A*˜*<sup>t</sup>*

. Then, it holds for all *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> that:

$$\mathbb{P}\left(Y\_{\mathbf{j}}, \mathbf{j} = 1, \ldots, n\right) \stackrel{\mathcal{D}}{=} \left(\tilde{A}! I\_{\mathbf{j}}^{\*}, \, \mathbf{j} = 1, \ldots, n\right). \tag{A.3}$$

Note that the assumption in (A2) is only well-defined because we assumed *r*(*p*,*q*)(*k*) < 1 for *k* ≥ 1 and *p*, *q* = 1, ... , *d* in (1). This becomes clear in the following considerations. In the proofs of the theorems in this chapter, we do not only need a decomposition of *Yj*, but also of *Yj*,*h*. As *Yj*,*<sup>h</sup>* is still a multivariate Gaussian process, the covariance matrix of *Yj*,*<sup>h</sup>*

given by Σ*d*,*<sup>h</sup>* is positive definite. Hence, it is possible to find an upper triangular matrix *A*, such that *AA<sup>t</sup>* = Σ*d*,*h*. It holds thatL

$$\chi\_{j,h} \stackrel{\mathcal{D}}{=} A \mathcal{U}\_{j,k}$$

for:

$$\mathcal{U}\_{j,\mathfrak{h}} = \left( \mathcal{U}\_{(j-1)\mathfrak{h}+1'}^{(1)}, \dots, \mathcal{U}\_{j\mathfrak{h}}^{(1)}, \dots, \mathcal{U}\_{(j-1)\mathfrak{h}+1'}^{(d)}, \dots, \mathcal{U}\_{j\mathfrak{h}}^{(d)} \right)^t.$$

The random vector *Uj*,*<sup>h</sup>* consists of (*d* · *h*) independent and standard normally distributed random variables. We notice the different structure of *Uj*,*<sup>h</sup>* compared to *Yj*,*h*. We assure that for consecutive *j*, the entries in *Uj*,*<sup>h</sup>* are all different while there are identical entries, for example in *Y*1,*<sup>h</sup>* = *Y*(1) <sup>1</sup> ,*Y*(1) <sup>2</sup> ,...,*Y*(*d*) *h t* and *Y*2,*<sup>h</sup>* = *Y*(1) <sup>2</sup> ,...,*Y*(*d*) *<sup>h</sup>* ,*Y*(*d*) *h*+1 *t* . This complicates our aim that:

$$\left(\boldsymbol{Y}\_{\boldsymbol{j},\boldsymbol{h}\prime} \boldsymbol{j} = 1, \ldots, \boldsymbol{n}\right)^{t} \stackrel{\mathcal{D}}{=} \left(\boldsymbol{A}\boldsymbol{U}\_{\boldsymbol{j},\boldsymbol{h}\prime} \boldsymbol{j} = 1, \ldots, \boldsymbol{n}\right)^{t} \tag{A4}$$

holds.

The special structure of *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> , namely that consisting of *h* consecutive entries of each marginal process *<sup>Y</sup>*(*p*) *j* , *p* = 1, ... , *d*, alongside the dependence between two random vectors in the process *Yj*,*<sup>h</sup>* , has to be reflected in the covariance matrix of *Uj*,*h*, *j* = 1, . . . , *n* . Hence, we need to check whether such a vector *Uj*,*h*, *j* = 1, . . . , *n* exists, i.e., if there is a positive semi-definite matrix that fulfills these conditions. We define **A** as a block diagonal matrix with *A* as main-diagonal blocks and all off-diagonal blocks as *dh* × *dh*-zero matrix.

We denote the covariance matrix of *Yj*,*h*, *j* = 1, . . . , *n t* by Σ*Y*,*<sup>n</sup>* and define the following matrix:

$$
\Sigma\_{lI,\mathfrak{u}} := \text{inv}(\mathbf{A}) \Sigma\_{\mathbf{Y},\mathfrak{u}} \text{inv}(\mathbf{A}^{\mathfrak{f}}).\tag{A5}
$$

We know that <sup>Σ</sup>*Y*,*<sup>n</sup>* is a positive semi-definite for all *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> because - *Yj* is a Gaussian process. Mathematically described, this means that:

$$\mathbf{x}^{\dagger} \Sigma\_{Y,\mathcal{U}} \mathbf{x} \ge \mathbf{0},\tag{A6}$$

for all *x* = (*x*1,..., *xnhd*) *<sup>t</sup>* <sup>∈</sup> <sup>R</sup>*nhd*. We conclude:

$$\begin{aligned} \mathbf{x}^t \boldsymbol{\Sigma}\_{lI, \mathbf{u}} \mathbf{x} &= \mathbf{x}^t \text{inv}(\mathbf{A}) \boldsymbol{\Sigma}\_{Y, \mathbf{u}} \text{inv} \begin{pmatrix} \mathbf{A}^t \end{pmatrix} \mathbf{x} \\ &= \begin{pmatrix} \text{inv} \left( \mathbf{A}^t \right) \mathbf{x} \end{pmatrix}^t \boldsymbol{\Sigma}\_{Y, \mathbf{u}} \begin{pmatrix} \mathbf{x}^t \text{inv}(\mathbf{A}) \end{pmatrix}, \\ \mathbf{x}^{\text{(A6)}} &\geq 0. \end{aligned}$$

Therefore, <sup>Σ</sup>*U*,*<sup>n</sup>* is a positive semi-definite matrix for all *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> and the random vector:

$$\left(\mathcal{U}\_{j,h\prime} \circ = 1, \ldots, n\right)^t \mathcal{N} \sim \left(0, \Sigma\_{\mathcal{U},h}\right)^t$$

exists and (A4) holds. Note that we do not have any further information on the dependence structure within the process - *Uj* , in general, this process neither exhibits long-range dependence, nor independence, nor stationarity.

We continue with two preparatory results that are also necessary for proving Theorem 2.1.

**Lemma A1.** *Let* - *Yj <sup>j</sup>*∈<sup>Z</sup> *be a d-dimensional Gaussian process as defined in (1) that fulfills (2) with d*<sup>1</sup> = ... = *dd* = *d*∗*, such that:*

$$
\Gamma\_Y(k) = \mathbb{E}\left(\mathcal{Y}\_{\vec{\jmath}}\mathcal{Y}\_{\vec{\jmath}+k}^t\right) \simeq Lk^{2d^\*-1}, \quad (k \to \infty).
$$

*Let C*<sup>2</sup> *be a normalization constant:*

$$\mathcal{C}\_2 = \frac{1}{2d^\*(4d^\*-1)}$$

*and let BY be an upper triangular matrix, such that:*

$$B\_Y B\_Y^t = L.$$

*Furthermore, for l* <sup>∈</sup> <sup>N</sup> *we have:*

$$\Upsilon\_{Y\mu}(l) = \frac{1}{n-l} \sum\_{j=1}^{n-l} \Upsilon\_j \Upsilon\_{j+l}^t$$

.

*Then, for h* <sup>∈</sup> <sup>N</sup> *it holds that:*

$$\begin{aligned} \left(n^{1-2d^\*} \left(\mathbb{C}\_2\right)^{-1/2} \left(B\_Y \otimes B\_Y\right)^{-1} \text{vec} \left(\mathbb{\hat{\Gamma}\_n(l) - \Gamma(l)\right), \ l = 0, \dots, h - 1}\right) \\ \stackrel{D}{\to} \left(\text{vec} \left(Z\_{2, d^\* + 1/2}^{(p, q)}(1)\right)\_{p, q = 1, \dots, d'}, l = 0, \dots, h - 1\right), \end{aligned}$$

*where Z*(*p*,*q*) 2,*d*∗+1/2(1) *has the spectral domain representation:*

$$Z\_{2,d^\*+1/2}^{(p,q)}(1) = K\_{p,q}(d^\*) \int\_{\mathbb{R}^2}^{\prime\prime} \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*} \mathcal{B}^{(p)}(\mathbf{d}\lambda\_1) \mathcal{B}^{(q)}(\mathbf{d}\lambda\_2)$$

*where:*

$$K\_{p,q}^2(d^\*) = \begin{cases} \begin{array}{c} \frac{1}{2\mathbb{C}\_2\left(2\Gamma\left(1-2d^\*\right)\sin\left(\pi d^\*\right)\right)^2}, & p=q\\ \frac{1}{\mathbb{C}\_2\left(2\Gamma\left(1-2d^\*\right)\sin\left(\pi d^\*\right)\right)^2}, & p\neq q. \end{array} \end{cases}$$

*and <sup>B</sup>*˜(d*λ*) = *B*˜(1)(d*λ*),..., *B*˜(*d*)(d*λ*) *is a multivariate Hermitian–Gaussian random measure as defined in (12).*

**Proof.** First, we can use (A1):

$$\chi\_j \stackrel{\mathcal{D}}{=} \tilde{A} \mathcal{U}\_j^\* \mathcal{V}$$

such that *U*∗ *j* is a multivariate Gaussian process with *U*∗ *<sup>j</sup>* ∼ N (0, *Id*) and *U*∗ *j* is still long-range dependent, as can be seen in (A2). It is possible to decompose the sample cross-covariance matrix <sup>Γ</sup>ˆ*Y*,*n*(*l*) <sup>−</sup> <sup>Γ</sup>*Y*(*l*) with respect to - *Yj* at lag *l* given by

$$
\Gamma\_{Y,n}(l) - \Gamma\_Y(l) = \frac{1}{n-l} \sum\_{j=1}^{n-l} \Upsilon\_{\vec{j}} \Upsilon\_{\vec{j}+l}^t - \mathbb{E}\left(\mathcal{Y}\_{\vec{j}} \mathcal{Y}\_{\vec{j}+l}^t\right),
$$

to:

$$
\Gamma\_{Y,n}(l) - \Gamma\_Y(l) \stackrel{\mathcal{D}}{=} \check{A} \left( \Gamma\_{U^\*,n}(l) - \Gamma\_{U^\*}(l) \right) \check{A}^t,
$$

where we define the sample cross-covariance matrix <sup>Γ</sup><sup>ˆ</sup> *<sup>U</sup>*∗,*n*(*l*) <sup>−</sup> <sup>Γ</sup>*U*<sup>∗</sup> (*l*) with respect to *U*∗ *j* at lag *l* by

$$
\hat{\Gamma}\_{\mathcal{U}^\*,\mathcal{U}}(l) - \Gamma\_{\mathcal{U}^\*}(l) = \frac{1}{n-l} \sum\_{j=1}^{n-l} \mathcal{U}\_j^\* \mathcal{U}\_{j+l}^\* - \mathbb{E}\left(\mathcal{U}\_j^\* \mathcal{U}\_{j+l}^\*\right).
$$

Each entry of:

$$
\Gamma\_{ll^\*,n}(l) - \Gamma\_{ll^\*}(l) = \left( \mathfrak{f}\_{n,\tilde{U}}^{(p,q)}(l) - r\_{\tilde{U}^\*}^{(p,q)}(l) \right)\_{p,q=1,\dots,d}
$$

is given by

$$\left(\mathcal{I}\_{n,II^\*}^{(p,q)}(l) - r\_{II^\*}^{(p,q)}(l)\right) := \sum\_{j=1}^n \mathcal{U}\_j^{\*\ (p)} \mathcal{U}\_{j+l}^{\*\ (q)} - \mathbb{E}\left(\mathcal{U}\_j^{\*\ (p)} \mathcal{U}\_{j+l}^{\*\ (q)}\right).$$

Following [31], proof of Lemma 7.4, the limit distribution of:

$$\left(\Gamma\_{ll^\*,n}(l) - \Gamma\_{ll^\*}(l), \ l = 0, \dots, h - 1\right)$$

is equal to the limit distribution of:

$$(\Upsilon\_{L^{\mathbf{L}^\*,\mathbf{n}}}(0) - \Gamma\_{L^{\mathbf{L}^\*}}(0), \; l = 0, \dots, h - 1).$$

We recall the assumption that *d*<sup>∗</sup> = *dp* for all *p* = 1, ... , *d*. We follow [14], Theorem 6 and use the Cramer–Wold device: Let *<sup>a</sup>*1,1, *<sup>a</sup>*1,2, ... , *ad*,*<sup>d</sup>* <sup>∈</sup> <sup>R</sup>. We are interested in the asymptotic behavior of:

$$\begin{aligned} &n^{1-2d^\*}\sum\_{p,q=1}^d a\_{p,q} \Big(\hat{r}\_{n,\mathcal{U}}^{(p,q)}(0) - r\_{\mathcal{U}}^{(p,q)}(0)\Big) \\ &= n^{-2d^\*}\sum\_{j=1}^n \sum\_{p,q=1}^d a\_{p,q} \Big(\mathcal{U}\_j^{\* \ (p)} \mathcal{U}\_j^{\* \ (q)} - \mathbb{E}\Big(\mathcal{U}\_j^{\* \ (p)} \mathcal{U}\_j^{\* \ (q)}\Big)\Big). \end{aligned}$$

We consider the function:

$$f\left(\mathcal{U}\_{\vec{j}}^{\*}\right) = \sum\_{p,q=1}^{d} a\_{p,q} \left(\mathcal{U}\_{\vec{j}}^{\*\ \left(p\right)}\mathcal{U}\_{\vec{j}}^{\*\ \left(q\right)} - \mathbb{E}\left(\mathcal{U}\_{\vec{j}}^{\*\ \left(p\right)}\mathcal{U}\_{\vec{j}}^{\*\ \left(q\right)}\right)\right) \tag{A7}$$

and may apply Theorem 6 in [14]. Using the Hermite decomposition of *f* as given in (11), we observe that *f* and therefore, *ap*,*q*, *p*, *q* = 1, ... , *d*, only affects the Hermite coefficients. Indeed, using [15], Lemma 3.5, the Hermite coefficients reduce to *ap*,*<sup>q</sup>* for each summand on the right-hand side in (A7). Hence, we can state:

$$m^{-2d^\*} \sum\_{j=1}^n \sum\_{p,q=1}^d a\_{p,q} \left( \mathcal{U}\_j^{\*\ (p)} \mathcal{U}\_j^{\*\ (q)} - \mathbb{E} \left( \mathcal{U}\_j^{\*\ (p)} \mathcal{U}\_j^{\*\ (q)} \right) \right) \tag{A8}$$

$$\stackrel{\mathcal{D}}{\rightarrow} \sum\_{p,q=1}^{d} a\_{p,q} Z\_{2,d^\* + 1/2}^{(p,q)}(1),\tag{A9}$$

where *<sup>Z</sup>*(*p*,*q*) 2,*d*∗+1/2(1) has the spectral domain representation:

$$Z\_{2,d^\*+1/2}^{(p,q)}(1) = K\_{p,q}(d^\*) \int\_{\mathbb{R}^2} \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*} \tilde{B}^{(p)}(\mathrm{d}\lambda\_1) \tilde{B}^{(q)}(\mathrm{d}\lambda\_2) \tag{A10}$$

where:

$$K\_{p,q}^2(d^\*) = \begin{cases} \frac{1}{2C\_2(2\Gamma(1-2d^\*)\sin(\pi d^\*))^2} & p=q\\\frac{1}{C\_2(2\Gamma(1-2d^\*)\sin(\pi d^\*))^2} & p \neq q. \end{cases}$$

and *<sup>B</sup>*˜(d*λ*) = *B*˜(1)(d*λ*),..., *B*˜(*d*)(d*λ*) is an appropriate multivariate Hermitian–Gaussian random measure. Thus, we proved convergence in the distribution of the sample-cross correlation matrix:

$$m^{1-2d^\*} \left( \mathbb{\hat{\Gamma}}\_{\operatorname{LI^\*},n}(0) - \Gamma\_{\operatorname{LI^\*}}(0) \right) \xrightarrow{D} \left( Z^{(p,q)}\_{2,d^\* + 1/2}(1) \right)\_{p,q=1,\dots,d}.$$

We take a closer look at the covariance matrix of vec- <sup>Γ</sup><sup>ˆ</sup> *<sup>U</sup>*∗,*n*(0) <sup>−</sup> <sup>Γ</sup>*U*<sup>∗</sup> (0) . Following [28], Lemma 5.7, we observe:

$$\begin{aligned} &n^{1-2d^\*} \left(4d^\*(4d^\*-1)\right)^{1/2} \text{Cov}\left(\text{vec}\left(\mathbb{T}\_{\varPi^\*,n}(0)-\varGamma\_{\varPi^\*}(0)\right), \text{vec}\left(\mathbb{T}\_{\varPi^\*,n}(0)-\varGamma\_{\varPi^\*}(0)\right)\right) \\ &= (I\_{d^2}+K\_{d^2})(L\_{\varPi^\*}\otimes L\_{\varPi^\*}), \end{aligned}$$

with *LU*<sup>∗</sup> as defined in (A3) and ⊗ denotes the Kronecker product. Furthermore, *Kd* denotes the commutation matrix that that transforms vec(*A*) into vec- *At* for *<sup>A</sup>* <sup>∈</sup> <sup>R</sup>*d*×*d*. Further details can be found in [32].

Hence, the covariance matrix of the vector of the sample cross-covariances is fully specified by the knowledge of *LU*<sup>∗</sup> as it arises in the context of long-range dependence in (A3).

We obtain a relation between *L* and *LU*<sup>∗</sup> , since:

$$
\Gamma\_Y(\cdot) = \vec{A} \Gamma\_U(\cdot) \vec{A}^t.
$$

Both:

$$
\Gamma\_Y(k) \simeq L k^{2d^\*-1} \left( k \to \infty \right),
$$

and:

$$
\Gamma\_{ll^\*} (k) \simeq L\_{ll^\*} k^{2d^\*-1} \ (k \to \infty),
$$

hold and we obtain:

$$L = \vec{A}L\_{U^\*}\vec{A}^t.$$

We study the covariance matrix of: vec- <sup>Γ</sup>ˆ*Y*,*n*(0) <sup>−</sup> <sup>Γ</sup>*Y*(0) :

$$\begin{split} &n^{1-2d^\*} \left(4d^\*(4d^\*-1)\right)^{1/2} \text{Cov}\left(\text{vec}(\text{Pr}\_{\mathcal{I},\mathbb{R}}(0)-\Gamma\_{\mathcal{I}}(0)\right), \text{vec}(\text{Pr}\_{\mathcal{I},\mathbb{R}}(0)-\Gamma\_{\mathcal{I}}(0)\right)^t \\ &\to \quad \left(I\_{d^2}+K\_{d^2}\right) \left(L\otimes L\right) \\ &= \quad \left(I\_{d^2}+K\_{d^2}\right) \left(\bar{A}L\_{l\mathcal{I}^\*}\bar{A}^t\right)\otimes \left(\bar{A}L\_{l\mathcal{I}^\*}\bar{A}^t\right) \\ &= \quad \left(I\_{d^2}+K\_{d^2}\right) \left(\bar{A}\otimes\bar{A}\right)\cdot \left(L\_{l\mathcal{I}^\*}\otimes L\_{l\mathcal{I}^\*}\right)\cdot \left(\bar{A}^t\otimes\bar{A}^t\right). \end{split} \tag{A11}$$

Let *BU*<sup>∗</sup> be an upper triangular matrix, such that:

$$B\_{LU^\*} B\_{LU^\*}^t := L\_{LU^\*}.$$

We know that such a matrix exists because *LU*<sup>∗</sup> is positive definite. Analogously, we define *BY*:

$$B\_Y := AB\_{U^\*} \cdot$$

Then, it holds that:

$$B\_Y B\_Y^t = L.$$

We arrive at:

$$\begin{split} &n^{1-2d^\*} \left(\mathrm{C}\_2\right)^{-1/2} \left(\mathrm{B}\_Y \otimes \mathrm{B}\_Y\right)^{-1} \mathrm{vec} \left(\widehat{\Gamma}\_{Y,\mathfrak{n}}(0) - \Gamma\_Y(0)\right) \\ &\stackrel{\scriptstyle D}{=} n^{1-2d^\*} \left(\mathrm{C}\_2\right)^{-1/2} \left(\mathrm{B}\_{\mathrm{I}^\*} \otimes \mathrm{B}\_{\mathrm{II}^\*}\right)^{-1} \left(A \otimes A\right)^{-1} \mathrm{vec} \left(\tilde{A} \left(\Gamma\_{\mathrm{II}^\*,\mathfrak{n}}(0) - \Gamma\_{\mathrm{II}^\*}(0)\right) \tilde{A}^{\mathfrak{t}}\right) \\ &= n^{1-2d^\*} \left(\mathrm{C}\_2\right)^{-1/2} \left(\mathrm{B}\_{\mathrm{II}^\*} \otimes \mathrm{B}\_{\mathrm{II}^\*}\right)^{-1} \mathrm{vec} \left(\Gamma\_{\mathrm{II}^\*,\mathfrak{n}}(0) - \Gamma\_{\mathrm{II}^\*}(0)\right) \\ &\stackrel{\scriptstyle D}{\to} \mathrm{vec} \left(Z\_{2,d^\*}^{(p,q)} + 1/2\right) \Biggr)\_{p,q=1,\ldots,d'} \end{split}$$

where *<sup>Z</sup>*(*p*,*q*) 2,*d*∗+1/2(1) has the spectral domain representation:

$$Z\_{2,d^\*+1/2}^{(p,q)}(1) = K\_{p,q}(d^\*) \int\_{\mathbb{R}^2}^{\prime\prime} \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*} \mathcal{B}^{(p)}(\mathbf{d}\lambda\_1) \mathcal{B}^{(q)}(\mathbf{d}\lambda\_2)$$

where:

$$K\_{p,q}^2(d^\*) = \begin{cases} \frac{1}{2\mathbb{C}\_2(2\Gamma(1-2d^\*)\sin(\pi d^\*))^2}, & p=q\\\frac{1}{\mathbb{C}\_2(2\Gamma(1-2d^\*)\sin(\pi d^\*))^2}, & p \neq q. \end{cases}$$

and *<sup>B</sup>*˜(d*λ*) = *B*˜(1)(d*λ*),..., *B*˜(*d*)(d*λ*) is a multivariate Hermitian–Gaussian random measure as defined in (12). Note that the standardization on the left-hand side is appropriate since the covariance matrix of vec(*Z*2,*d*∗+1/2(1)) is given by

$$\begin{split} \mathbb{E}\left(\boldsymbol{K}^{2}(\boldsymbol{d}^{\*})\int\_{\mathbb{R}^{2}}^{\prime\prime} \int\_{\mathbb{R}^{2}} \boldsymbol{E}\_{\lambda\_{1},\lambda\_{2}} \overline{\boldsymbol{E}\_{\lambda\_{3},\lambda\_{4}}} \text{vec}\left(\boldsymbol{\tilde{\boldsymbol{B}}}(\mathrm{d}\lambda\_{1})\left(\boldsymbol{\tilde{\boldsymbol{B}}}(\mathrm{d}\lambda\_{2})\right)^{\prime}\right) \\\\ \left(\mathrm{vec}\left(\overline{\boldsymbol{\tilde{\boldsymbol{B}}}(\mathrm{d}\lambda\_{3})\left(\overline{\boldsymbol{\tilde{\boldsymbol{B}}}(\mathrm{d}\lambda\_{4})}\right)^{\prime}}\right)\right)^{\prime} \end{split} \tag{A12}$$

by denoting:

$$E\_{\lambda\_1, \lambda\_2} := \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*}.$$

We observe:

$$\begin{aligned} &\mathbb{E}\left(\text{vec}\left(\bar{B}(\text{d}\lambda\_1)\bar{B}(\text{d}\lambda\_2)^t\right)\left(\text{vec}\left(\overline{\bar{B}(\text{d}\lambda\_3)\left(\bar{B}(\text{d}\lambda\_4)\right)^t}\right)\right)^t\right) \\ &= \left\{ \begin{array}{ll} I\_d \text{2d}\lambda\_1 \text{d}\lambda\_2, & |\lambda\_1| = |\lambda\_3| \neq |\lambda\_2| = |\lambda\_4|, \\\ &\mathbb{K}\_{d^2} \text{d}\lambda\_1 \text{d}\lambda\_2, & |\lambda\_1| = |\lambda\_4| \neq |\lambda\_2| = |\lambda\_3|, \end{array} \right. \end{aligned} \tag{A13}$$

following [28], (27). Neither the case |*λ*1| = |*λ*2| nor |*λ*3| = |*λ*4| has to be incorporated as the diagonals are excluded in the integration in (A12).

**Corollary A1.** *Under the assumptions of Lemma A1, there is a different representation of the limit random vector. For h* <sup>∈</sup> <sup>N</sup>*, we obtain:*

$$\left( \left( n^{1 - 2d^\*} \left( \mathbb{C}\_2 \right)^{-1/2} \text{vec} \left( \mathbb{\Gamma}\_n(l) - \Gamma(l) \right) , \ l = 0, \dots, h - 1 \right) \xrightarrow{D} \left( \text{vec} \left( \mathbb{Z}\_{2, d^\* + 1/2}(1) \right) \, l = 0, \dots, h - 1 \right) ,$$

*where* vec(*Z*2,*d*∗+1/2(1)) *has the spectral domain representation:*

$$\text{vec}(Z\_{2,d^\*+1/2}(1)) = D\_{\mathcal{K}(d^\*)} \int\_{\mathbb{R}^2}^{\prime\prime} \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*} \text{vec}\left(\mathcal{B}\_L(\text{d}\lambda\_1)\mathcal{B}\_L(\text{d}\lambda\_2)^t\right).$$

*The matrix DK*(*d*∗) *is a diagonal matrix:*

$$D\_{K(d^\*)} = \text{diag}(\text{vec}(K(d^\*))),$$

*and K*(*d*∗) = - *Kp*,*q*(*d*∗) *<sup>p</sup>*,*q*=1,...,*<sup>d</sup> is such that:*

$$K^2(d^\*)\_{p,q} = \begin{cases} \frac{1}{2C\_2\left(2\Gamma(1-2d^\*)\sin\left(\pi d^\*\right)\right)^2}, & p=q\\\frac{1}{C\_2\left(2\Gamma(1-2d^\*)\sin\left(\pi d^\*\right)\right)^2}, & p\neq q. \end{cases}$$

*Furthermore, B*˜*L*(d*λ*) *is a multivariate Hermitian–Gaussian random measure that fulfills:*

$$\mathbb{E}\left(\mathcal{B}\_L(\mathbf{d}\lambda)\mathcal{B}\_L(\mathbf{d}\lambda)^\*\right) = L \,\mathbf{d}\lambda.$$

**Proof.** The proof is an immediate consequence of Lemma A1 using *B*˜*L*(d*λ*) = *BYB*˜(d*λ*) with *BYB<sup>t</sup> <sup>Y</sup>* <sup>=</sup> *<sup>L</sup>* and *<sup>B</sup>*˜(d*λ*) as defined in (12).

*Appendix A.2. Proof of Theorem 2.1*

**Proof.** Without loss of generality, we assume E *f Yj*,*<sup>h</sup>* <sup>=</sup> 0. Following the argumentation in [20], Theorem 5.9, we first remark that *Yj*,*<sup>h</sup>* D = *AUj*,*<sup>h</sup>* with *Uj*,*<sup>h</sup>* and *A* as described in (A4) and (A5). We now want to study the asymptotic behavior of the partial sum *n* ∑ *j*=1 *f* <sup>∗</sup>- *Uj* where *f* ∗ *Uj*,*<sup>h</sup>* := *f AUj*,*<sup>h</sup>* <sup>D</sup> = *f Yj*,*<sup>h</sup>* . Since *m*(*f* <sup>∗</sup>, *Idh*) = *m*(*f* ◦ *A*, *Idh*) = *m*(*f* , Σ*d*,*h*) = 2, as can be seen in [15], Lemma 3.7, hence, we know by [14], Theorem 6, that these partial sums are dominated by the second order terms in the Hermite expansion of *f* ∗:

$$\sum\_{j=1}^{n} f^\* \left( \mathcal{U}\_{j,h} \right) \sum\_{j=1}^{n} \sum\_{l\_1 + \ldots + l\_{dh} = 2} \mathbb{E} \left( f^\* \left( \mathcal{U}\_{j,h} \right) H\_{l\_1, \ldots, l\_{dh}} \left( \mathcal{U}\_{j,h} \right) \right) H\_{l\_1, \ldots, l\_{dh}} \left( \mathcal{U}\_{j,h} \right) \\ + o\_{\mathbb{P}} \left( n^{2d^\*} \right).$$

This follows from the multivariate extension of the Reduction Theorem as proven in [14]. We obtain:

$$\begin{split} &\sum\_{l\_{1}+\ldots+l\_{h}=2} \mathbb{E}\left(f^{\*}\left(\mathcal{U}\_{j,h}\right)H\_{l\_{1},\ldots,l\_{h}}\left(\mathcal{U}\_{j,h}\right)\right)H\_{l\_{1},\ldots,l\_{h}}\left(\mathcal{U}\_{j,h}\right) \\ = &\sum\_{l=1}^{dh} \mathbb{E}\left(f^{\*}\left(\mathcal{U}\_{j,h}\right)\left(\left(\mathcal{U}\_{j,h}^{(l)}\right)^{2}-1\right)\right)\left(\left(\mathcal{U}\_{j,h}^{(l)}\right)^{2}-1\right) + \sum\_{1\leq j,k\leq dh, i\neq k} \mathbb{E}\left(f^{\*}\left(\mathcal{U}\_{j,h}\right)\mathcal{U}\_{j,h}^{(l)}\mathcal{U}\_{j,h}^{(k)}\right)\mathcal{U}\_{j,h}^{(l)}\mathcal{U}\_{j,h}^{(k)} \\ = &\sum\_{l=1}^{dh} \mathbb{E}\left(f^{\*}\left(\mathcal{U}\_{j,h}\right)\left(\mathcal{U}\_{j,h}^{(l)}\right)^{2}\right)\left(\left(\mathcal{U}\_{j,h}^{(l)}\right)^{2}-1\right) + \sum\_{1\leq j,k\leq dh, i\neq k} \mathbb{E}\left(f^{\*}\left(\mathcal{U}\_{j,h}\right)\mathcal{U}\_{j,h}^{(l)}\mathcal{U}\_{j,h}^{(k)}\right)\mathcal{U}\_{j,h}^{(l)}\mathcal{U}\_{j,h}^{(k)}.\end{split}$$

since E *f* ∗ *Uj*,*<sup>h</sup>* <sup>=</sup> <sup>E</sup> *f Yj*,*<sup>h</sup>* <sup>=</sup> 0. This results in:

$$\sum\_{l=1}^{dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{l,h}\right)\left(\mathcal{U}\_{l,h}^{(l)}\right)^2\right) \left(\left(\mathcal{U}\_{l,h}^{(l)}\right)^2 - 1\right) + \sum\_{1 \le j,k \le dh, l \ne k} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{l,h}\right)\mathcal{U}\_{j,h}^{(l)}\mathcal{U}\_{j,h}^{(k)}\right) \mathcal{U}\_{j,h}^{(l)} \mathcal{U}\_{j,h}^{(k)}$$

$$= \sum\_{1 \le j,k \le dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,h}\right)\mathcal{U}\_{j,h}^{(l)}\mathcal{U}\_{j,h}^{(k)}\right) \mathcal{U}\_{j,h}^{(l)} \mathcal{U}\_{j,h}^{(k)} - \sum\_{l=1}^{dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{l,h}\right)\left(\mathcal{U}\_{l,h}^{(l)}\right)^2\right). \tag{A14}$$

Note that:

$$\sum\_{1 \le i,k \le dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,h}\right)\mathcal{U}\_{j,h}^{(i)}\mathcal{U}\_{j,h}^{(k)}\right)\mathbb{E}\left(\mathcal{U}\_{j,h}^{(i)}\mathcal{U}\_{j,h}^{(k)}\right) = \sum\_{i=1}^{dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,h}\right)\left(\mathcal{U}\_{j,h}^{(i)}\right)^2\right) \tag{A15}$$

since the entries of *Uj*,*<sup>h</sup>* are independent for fixed *j* and identically N (0, 1) distributed. Thus, the subtrahend in (A14) equals the expected value of the minuend.

Define *<sup>B</sup>* :<sup>=</sup> (*bi*,*k*)1≤*i*,*k*≤*dh* <sup>∈</sup> <sup>R</sup>(*dh*)×(*dh*) with:

$$b\_{i,k} := \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,\mathcal{h}}\right)\mathcal{U}\_{j,\mathcal{h}}^{(i)}\mathcal{U}\_{j,\mathcal{h}}^{(k)}\right) = \mathbb{E}\left(f^\*\left(\mathcal{U}\_1\right)\mathcal{U}\_1^{(i)}\mathcal{U}\_1^{(k)}\right),$$

since we are considering a stationary process. We obtain:

$$B = \mathbb{E}\left(\mathcal{U}\_{\bar{j},h}f^\*\left(\mathcal{U}\_{\bar{j},h}\right)\mathcal{U}\_{\bar{j},h}^t\right) = \mathbb{E}\left(A^{-1}\mathcal{Y}\_{\bar{j},h}f\left(\mathcal{Y}\_{\bar{j},h}\right)\mathcal{Y}\_{\bar{j},h}^t\left(A^{-1}\right)^t\right).$$

Hence, we can state the following:

$$\begin{split} &\sum\_{1\le j,k\le dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{jh}\right)\mathcal{U}\_{jh}^{(i)}\mathcal{U}\_{jh}^{(k)}\right)\mathcal{U}\_{jh}^{(i)}\mathcal{U}\_{jh}^{(k)} \\ &=\mathcal{U}\_{jh}^{\dagger}B\mathcal{U}\_{jh} \\ &\overset{\mathcal{D}}{=}\mathcal{Y}\_{jh}^{\dagger}\left(A^{-1}\right)^{\dagger}BA^{-1}\mathcal{Y}\_{j,h} \\ &=\mathcal{Y}\_{jh}^{\dagger}\left(A^{-1}\right)^{\dagger}A^{-1}\mathbb{E}\left(\mathcal{Y}\_{jh}\mathcal{f}\left(\mathcal{Y}\_{jh}\right)\mathcal{Y}\_{jh}^{\dagger}\right)\left(A^{-1}\right)^{\dagger}A^{-1}\mathcal{Y}\_{j,h} \\ &=\mathcal{Y}\_{jh}^{\dagger}\Sigma\_{d,h}^{-1}\mathbb{E}\left(\mathcal{Y}\_{jh}\mathcal{f}\left(\mathcal{Y}\_{jh}\right)\mathcal{Y}\_{j,h}^{\dagger}\right)\Sigma\_{d,h}^{-1}\mathcal{Y}\_{j,h} \\ &=\mathcal{Y}\_{jh}^{\dagger}\mathbf{A}\mathcal{Y}\_{j,h} \\ &=\sum\_{1\le j,k\le dh} \mathcal{Y}\_{j}^{(i)}\mathcal{Y}\_{j}^{(k)}\pi\_{ik\nu} \end{split} \tag{A16}$$

,

where we define <sup>A</sup> :<sup>=</sup> (*αik*)1≤*ik*≤*dh* :<sup>=</sup> <sup>Σ</sup>−<sup>1</sup> *d*,*hC*Σ−<sup>1</sup> *<sup>d</sup>*,*h*, with *<sup>C</sup>* :<sup>=</sup> <sup>E</sup> *Yj*,*<sup>h</sup> f Yj*,*<sup>h</sup> Yt j*,*h* as the matrix of second order Hermite coefficients, in contrast to *B* now with respect to the original considered process *Yj*,*<sup>h</sup> <sup>j</sup>*∈<sup>Z</sup> .

Remembering the special structure of *Yj*,*<sup>h</sup>* = *Y*(1) *<sup>j</sup>* ,...,*Y*(1) *<sup>j</sup>*+*h*−1,...,*Y*(*d*) *<sup>j</sup>* ,*Y*(*d*) *j*+*h*−1 *t* namely that *Y*(*k*) *<sup>j</sup>*,*<sup>h</sup>* <sup>=</sup> *<sup>Y</sup>*( *<sup>k</sup>*−<sup>1</sup> *<sup>h</sup>* +1) *<sup>j</sup>*+(*<sup>k</sup>* mod *<sup>h</sup>*)−1, *<sup>k</sup>* <sup>=</sup> 1, . . . , *dh*, we can see that:

$$\sum\_{j=1}^{n} \sum\_{1 \le \vec{k} \le dh} Y\_{j,h}^{(i)} Y\_{j,h}^{(k)} \mathfrak{a}\_{ik} = \sum\_{j=1}^{n} \sum\_{1 \le \vec{k} \le dh} Y\_{j+(i \mod h)-1}^{\left(\left\lfloor \frac{k-1}{h} \right\rfloor + 1\right)}^{\left(\left\lfloor \frac{k-1}{h} \right\rfloor + 1\right)} Y\_{j+(k \mod h)-1}^{\left(\left\lfloor \frac{k-1}{h} \right\rfloor + 1\right)} \mathfrak{a}\_{ik}$$

$$= \sum\_{j=1}^{n} \sum\_{p,q=1}^{d} \sum\_{i,k=1}^{h} Y\_{j+i-1}^{(p)} Y\_{j+k-1}^{(q)} \mathfrak{a}\_{ik}^{(p,q)},\tag{A17}$$

,

where we divide:

$$\mathbb{A} = \begin{pmatrix} \mathbb{A}^{(1,1)} & \mathbb{A}^{(1,2)} & \dots & \mathbb{A}^{(1,d)} \\ \mathbb{A}^{(2,1)} & \mathbb{A}^{(2,2)} & \dots & \mathbb{A}^{(2,d)} \\ \vdots & \vdots & \vdots \\ \mathbb{A}^{(d,1)} & \mathbb{A}^{(d,2)} & \dots & \mathbb{A}^{(d,d)} \end{pmatrix}.$$

with A(*p*,*q*) = *α* (*p*,*q*) *i*,*k* 1≤*i*,*k*≤*h* <sup>∈</sup> <sup>R</sup>*h*×*<sup>h</sup>* such that *<sup>α</sup>* (*p*,*q*) *<sup>i</sup>*,*<sup>k</sup>* = *<sup>α</sup>i*+(*p*−1)*h*,*k*+(*q*−1)*<sup>h</sup>* for each *<sup>p</sup>*, *<sup>q</sup>* = 1, . . . , *d* and *i*, *k* = 1, . . . , *h*.

We can now split the considered sum in (A17) in a way such that we are afterwards able to express it in terms of sample cross-covariances. In order to do so, we define the sample cross-covariance at lag *l* by

$$\mathfrak{r}\_n^{(p,q)}(l) := \frac{1}{n} \sum\_{j=1}^{n-l} X\_j^{(p)} X\_{j+l}^{(q)}$$

for *p*, *q* = 1, . . . , *d*.

Note that in the case *h* = 1, it follows directly that:

$$\sum\_{p=1}^{n} \sum\_{p,q=1}^{d} \sum\_{i,k=1}^{h} \chi\_{j+i-1}^{(p)} \chi\_{j+k-1}^{(q)} a\_{ik}^{(p,q)} = \sum\_{p,q=1}^{d} a\_{1,1}^{(p,q)} \sum\_{j=1}^{n} \chi\_{j}^{(p)} \chi\_{j}^{(q)} = n \sum\_{p,q=1}^{d} \hat{r}\_{n}^{(p,q)}(0).$$

The case *h* = 2 has to be regarded separately, too, and we obtain:

$$\begin{split} &\sum\_{j=1}^{n}\sum\_{p,q=1}^{d}\sum\_{j,k=1}^{2}Y\_{j+j-1}^{(p)}Y\_{j+k-1}^{(q)}a\_{ik}^{(p,q)}\\ &=\sum\_{p,q=1}^{d}\left(a\_{1,1}^{(p,q)}\sum\_{j=1}^{n}Y\_{j}^{(p)}Y\_{j}^{(q)} + a\_{1,2}^{(p,q)}\sum\_{j=1}^{n}Y\_{j}^{(p)}Y\_{j+1}^{(q)} + a\_{2,1}^{(p,q)}\sum\_{j=1}^{n}Y\_{j+1}^{(p)}Y\_{j}^{(q)} + a\_{2,2}^{(p,q)}\sum\_{j=1}^{n}Y\_{j+1}^{(p)}Y\_{j+1}^{(q)}\right)\\ &=\sum\_{p,q=1}^{d}\left(a\_{1,1}^{(p,q)}n\_{n}^{(p,q)}(0) + a\_{1,2}^{(p,q)}\left(n!\mathfrak{r}\_{n}^{(p,q)}(1) + \underbrace{Y\_{n}^{(p)}Y\_{n+1}^{(q)}}\_{\star\star}\right) + a\_{2,1}^{(p,q)}\left(n!\mathfrak{r}\_{n}^{(p,q)}(1) + \underbrace{Y\_{n+1}^{(p)}Y\_{n}^{(q)}}\_{\star\star}\right)\right)\\ &\qquad+ a\_{2,2}^{(p,q)}\left(n!\mathfrak{r}\_{n}^{(p,q)}(0) + \underbrace{Y\_{n+1}^{(p)}Y\_{n+1}^{(q)}}\_{\star\star} - \underbrace{Y\_{1}^{(p)}Y\_{1}^{(q)}}\_{\star\star}\right)\right),\end{split}$$

Note that for each of the terms labeled by -, the following holds for *d*<sup>∗</sup> ∈ 1 4 , 1 2 :

$$n^{-2d^\*} \star \stackrel{\mathbb{P}}{\rightarrow} 0, \ (n \rightarrow \infty).$$

We use this property later on when dealing with the asymptotics of the term in (A17). Finally, we consider the term in (A17) for *h* ≥ 3 and arrive at:

$$\begin{split} &\sum\_{j=1}^{n}\sum\_{p,q=1}^{d}\sum\_{i,k=1}^{h}Y\_{j+i-1}^{(p)}Y\_{j+k-1}^{(q)}a\_{ik}^{(p,q)}\\ &=\sum\_{p,q=1}^{d}\sum\_{i,k=1}^{h}a\_{ik}^{(p,q)}\sum\_{j=i}^{n+i-1}Y\_{j}^{(p)}Y\_{j+k-i}^{(q)}\\ &=\sum\_{p,q=1}^{d}\sum\_{i=0}^{h-1}\sum\_{i=1}^{h-1}a\_{i,i+l}^{(p,q)}\sum\_{j=i}^{n+i-1}Y\_{j}^{(p)}Y\_{j+l}^{(q)} + \sum\_{p,q=1}^{d}\sum\_{l=-(h-1)}^{-1}\sum\_{i=1}^{h}a\_{i,j+l}^{(p,q)}\sum\_{j=i}^{n+i-1}Y\_{j}^{(p)}Y\_{j+l}^{(q)}\\ &=\sum\_{p,q=1}^{d}\sum\_{l=0}^{h-1}\sum\_{i=0}^{h-1}a\_{i,i+l}^{(p,q)}\sum\_{j=i}^{n+i-1}Y\_{j}^{(p)}Y\_{j+l}^{(q)}\\ &+\sum\_{p,q=1}^{d}\sum\_{l=1}^{h-1}\sum\_{i=1}^{h-1}a\_{i,i+l}^{(p,q)}\sum\_{j=i}^{n+i-1}Y\_{j+l}^{(p)}Y\_{j}^{(q)} \end{split} \tag{A18}$$

= *d* ∑ *p*,*q*=1 *h* ∑ *i*=1 *α* (*p*,*q*) *i*,*i n*+*i*−1 ∑ *j*=*i <sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j* + *d* ∑ *p*,*q*=1 *h*−1 ∑ *l*=1 *h*−*l* ∑ *i*=1 \$ *α* (*p*,*q*) *i*,*i*+*l n*+*i*−1 ∑ *j*=*i <sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *<sup>j</sup>*+*<sup>l</sup>* + *α* (*p*,*q*) *i*+*l*,*i n*+*i*−1 ∑ *j*=*i <sup>Y</sup>*(*p*) *j*+*l <sup>Y</sup>*(*q*) *j* % (A19) = *d* ∑ *p*,*q*=1 \$ *α* (*p*,*q*) 1,1 *n* ∑ *j*=1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *<sup>j</sup>* + *h* ∑ *i*=2 *α* (*p*,*q*) *i*,*i n*+*i*−1 ∑ *j*=*i <sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j* % + *d* ∑ *p*,*q*=1 *h*−2 ∑ *l*=1 \$\$*α* (*p*,*q*) 1,1+*l n* ∑ *j*=1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *<sup>j</sup>*+*<sup>l</sup>* + *α* (*p*,*q*) 1+*l*,1 *n* ∑ *j*=1 *<sup>Y</sup>*(*p*) *j*+*l <sup>Y</sup>*(*q*) *j* % + *h*−*l* ∑ *i*=2 \$ *α* (*p*,*q*) *i*,*i*+*l n*+*i*−1 ∑ *j*=*i <sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *<sup>j</sup>*+*<sup>l</sup>* + *α* (*p*,*q*) *i*+*l*,*i n*+*i*−1 ∑ *j*=*i <sup>Y</sup>*(*p*) *j*+*l <sup>Y</sup>*(*q*) *j* %% + *d* ∑ *p*,*q*=1 \$ *α* (*p*,*q*) 1,*h n* ∑ *j*=1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *<sup>j</sup>*+*h*−<sup>1</sup> <sup>+</sup> *<sup>α</sup>* (*p*,*q*) *<sup>h</sup>*,1 *<sup>Y</sup>*(*p*) *<sup>j</sup>*+*h*−1*Y*(*q*) *j* % (A20) = *d* ∑ *p*,*q*=1 \$ *α* (*p*,*q*) 1,1 *nr*ˆ (*p*,*q*) *<sup>n</sup>* (0) + *h* ∑ *i*=2 *α* (*p*,*q*) *i*,*i* \$ *<sup>n</sup>*+*i*−<sup>1</sup> ∑ *j*=*n*+1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j* 1 23 4 - +*nr*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *i*−1 ∑ *j*=1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j* 1 23 4 - %% + *d* ∑ *p*,*q*=1 *h*−2 ∑ *l*=1 \$\$*α* (*p*,*q*) 1,1+*l nr*ˆ (*p*,*q*) *<sup>n</sup>* (*l*) + *<sup>α</sup>* (*p*,*q*) <sup>1</sup>+*l*,1*nr*ˆ (*q*,*p*) *<sup>n</sup>* (*l*) % + *h*−*l* ∑ *i*=2 \$ *α* (*p*,*q*) *i*,*i*+*l* \$ *<sup>n</sup>*+*i*−<sup>1</sup> ∑ *j*=*n*−*l*+1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j*+*l* 1 23 4 - +*nr*ˆ (*p*,*q*) *<sup>n</sup>* (*l*) <sup>−</sup> *i*−1 ∑ *j*=1 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j*+*l* 1 23 4 - % + *α* (*p*,*q*) *i*+*l*,*i* \$ *<sup>n</sup>*+*i*−<sup>1</sup> ∑ *j*=*n*−*l*+1 *<sup>Y</sup>*(*p*) *j*+*l <sup>Y</sup>*(*q*) *j* 1 23 4 - +*nr*ˆ (*q*,*p*) *<sup>n</sup>* (*l*) <sup>−</sup> *i*−1 ∑ *j*=1 *<sup>Y</sup>*(*p*) *j*+*l <sup>Y</sup>*(*q*) *j* 1 23 4 - %%% + *d* ∑ *p*,*q*=1 \$ *α* (*p*,*q*) 1,*h* \$ *n* ∑ *j*=*n*−*h*+2 *<sup>Y</sup>*(*p*) *<sup>j</sup> <sup>Y</sup>*(*q*) *j*+*h*−1 1 23 4 - +*nr*ˆ (*p*,*q*) *<sup>n</sup>* (*<sup>h</sup>* <sup>−</sup> <sup>1</sup>) % (A21) + *α* (*p*,*q*) *<sup>h</sup>*,1 \$ *<sup>n</sup>* ∑ *j*=*n*−*h*+2 *<sup>Y</sup>*(*p*) *<sup>j</sup>*+*h*−1*Y*(*q*) *j* 1 23 4 - +*nr*ˆ (*q*,*p*) *<sup>n</sup>* (*<sup>h</sup>* <sup>−</sup> <sup>1</sup>) %%. (A22)

Again for each of the terms labeled by it holds for *d*<sup>∗</sup> ∈ 1 4 , 1 2 :

$$n^{-2d^\*} \star \stackrel{\mathbb{P}}{\rightarrow} 0, \ (n \rightarrow \infty),$$

since each describes a sum with a finite number (independent of *n*) of summands. Therefore, we continue to express the terms denoted by by *o*P *n*2*d*<sup>∗</sup> .

With these calculations, we are able to re-express the partial sum, whose asymptotics we are interested in, in terms of the sample cross-correlations of the original long-range dependent process - *Yj j*∈Z.

Finally, the previous calculations lead to:

$$\begin{aligned} &\sum\_{j=1}^{n} f\left(Y\_{j,h}\right) \\ \stackrel{\scriptstyle D}{=} &\sum\_{j=1}^{n} f^\*\left(\mathcal{U}\_{j,h}\right) \\ \stackrel{\scriptstyle D}{=} &\sum\_{j=1}^{n} \left(\sum\_{1\le j,k\le \rm dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,h}\right)\mathcal{U}\_{j,h}^{(j)}\mathcal{U}\_{j,h}^{(k)}\right) \mathcal{U}\_{j,h}^{(j)} \mathcal{U}\_{j,h}^{(k)} - \sum\_{i=1}^{\rm dh} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,h}\right)\left(\mathcal{U}\_{j,h}^{(i)}\right)^2\right)\right) + o\_{\mathbb{P}}\left(n^{2d^\*}\right) \\ \stackrel{\scriptstyle D}{=} &\sum\_{j=1}^{n} \sum\_{p,q=1}^{n} \sum\_{i,k=1}^{h} a\_{ik}^{(p,q)} \left(Y\_{j+i-1}^{(p)} Y\_{j+k-1}^{(q)} - \mathbb{E}\left(Y\_{j+i-1}^{(p)} Y\_{j+k-1}^{(q)}\right)\right) + o\_{\mathbb{P}}\left(n^{2d^\*}\right), \end{aligned} \tag{A23}$$

where (A23) follows, since (A15) yields:

$$\begin{split} \sum\_{i=1}^{\mathrm{d}h} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,\mathrm{h}}\right)\left(\mathcal{U}\_{j,\mathrm{h}}^{(i)}\right)^2\right) &= \sum\_{1 \le i,j,k \le \mathrm{d}h} \mathbb{E}\left(f^\*\left(\mathcal{U}\_{j,\mathrm{h}}\right)\mathcal{U}\_{j,\mathrm{h}}^{(i)}\mathcal{U}\_{j,\mathrm{h}}^{(k)}\right) \mathbb{E}\left(\mathcal{U}\_{j,\mathrm{h}}^{(i)}\mathcal{U}\_{j,\mathrm{h}}^{(k)}\right),\\ &\stackrel{(\mathrm{A17})}{=} \sum\_{p,q=1}^{d} \sum\_{i,k=1}^{h} a\_{ik}^{(p,q)} \mathbb{E}\left(Y\_{j+i-1}^{(p)} Y\_{j+k-1}^{(q)}\right). \end{split}$$

Taking the parts containing the sample cross-correlations into account, we derive:

*n* ∑ *j*=1 *d* ∑ *p*,*q*=1 *h* ∑ *i*,*k*=1 *α* (*p*,*q*) *ik <sup>Y</sup>*(*p*) *<sup>j</sup>*+*i*−1*Y*(*q*) *<sup>j</sup>*+*k*−<sup>1</sup> <sup>−</sup> <sup>E</sup> *<sup>Y</sup>*(*p*) *<sup>j</sup>*+*i*−1*Y*(*q*) *j*+*k*−1 <sup>+</sup> *<sup>o</sup>*P(*n*2*d*<sup>∗</sup> ) (A22) = *d* ∑ *p*,*q*=1 \$ *α* (*p*,*q*) 1,1 *n r*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (0) + *h* ∑ *i*=2 *α* (*p*,*q*) *<sup>i</sup>*,*<sup>i</sup> n r*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (0) % + *d* ∑ *p*,*q*=1 *h*−2 ∑ *l*=1 \$ *α* (*p*,*q*) 1,1+*l n r*ˆ (*p*,*q*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (*l*) + *α* (*p*,*q*) <sup>1</sup>+*l*,1*n r*ˆ (*q*,*p*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*q*,*p*) (*l*) + *h*−*l* ∑ *i*=2 \$ *α* (*p*,*q*) *<sup>i</sup>*,*i*+*<sup>l</sup> n r*ˆ (*p*,*q*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (*l*) + *α* (*p*,*q*) *<sup>i</sup>*+*l*,*<sup>i</sup> n r*ˆ (*q*,*p*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*q*,*p*) (*l*) %% + *d* ∑ *p*,*q*=1 *α* (*p*,*q*) 1,*<sup>h</sup> n r*ˆ (*p*,*q*) *<sup>n</sup>* (*<sup>h</sup>* <sup>−</sup> <sup>1</sup>) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (*h* − 1) + *α* (*p*,*q*) *<sup>h</sup>*,1 *n r*ˆ (*q*,*p*) *<sup>n</sup>* (*<sup>h</sup>* <sup>−</sup> <sup>1</sup>) <sup>−</sup> *<sup>r</sup>*(*q*,*p*) (*h* − 1) (A24) + *o*P(*n*2*d*<sup>∗</sup> )

$$\begin{split} \mathcal{I} &= \; n \sum\_{p,q=1}^{d} \left( \sum\_{l=0}^{h-1} \sum\_{i=1}^{h-l} a\_{i,i+l}^{(p,q)} \left( \hat{\mathfrak{r}}\_{n}^{(p,q)}(l) - r^{(p,q)}(l) \right) + \sum\_{l=1}^{h-1} \sum\_{i=1}^{h-l} a\_{i+l,i}^{(p,q)} \left( \hat{\mathfrak{r}}\_{n}^{(q,p)}(l) - r^{(q,p)}(l) \right) \right) \\ &+ o\_{\mathbb{P}}(n^{2d^\*}). \end{split} \tag{A25}$$

We take a closer look at the impact of each long-range dependence parameter *dp*, *p* = 1, ... , *d* to the convergence of this sum. The setting we are considering does not allow for a normalization depending on *<sup>p</sup>* and *<sup>q</sup>* for each cross-correlation *r*ˆ (*p*,*q*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*p*,*q*)(*l*) , *l* = 0, ... , *h* − 1, but we need to find a normalization for all *p*, *q* = 1, ... , *d*. Hence, we need to remember the set *P*<sup>∗</sup> := {*p* ∈ {1, ... , *d*} : *dp* ≥ *dq*∀*q* ∈ {1, ... , *d*}} and the parameter

*d*∗ = max *p*=1,...,*d dp*, such that for each *p* ∈ *P*∗, we have *dp* = *d*∗. For each *p*, *q* ∈ {1, ... , *d*} with (*p*, *q*) ∈/ (*P*<sup>∗</sup> × *P*∗) and *l* = 0, . . . , *h* − 1, we conclude that:

$$\begin{split} \mathbb{E}\left(\left(n^{1-2d^\*}\left(\mathfrak{f}\_n^{\langle p,q\rangle}(l)-r^{\langle p,q\rangle}(l)\right)\right)^2\right) &= n^{2\langle d\_p+d\_q-2d^\*\rangle}\mathbb{E}\left(\left(n^{1-d\_p-d\_q}\left(\mathfrak{f}\_n^{\langle p,q\rangle}(l)-r^{\langle p,q\rangle}(l)\right)\right)^2\right) \\ &= n^{2d\_p+2d\_q-4d^\*}\mathbb{C}\_2\left(L\_{p,p}L\_{q,q}+L\_{p,q}L\_{q,p}\right) \end{split} \tag{A26}$$
 
$$\xrightarrow{(n\to\infty)} 0,\tag{A26}$$

since *dp* + *dq* − 2*d*<sup>∗</sup> < 0.

This implies that:

$$n^{1-2d^\*} \left( \mathfrak{f}\_n^{(p,q)}(0) - r^{(p,q)}(0) \right) \xrightarrow{\mathbb{P}} 0$$

Hence, using Slutsky's theorem, the crucial parameters that determine the normalization and therefore, the limit distribution of (A27) are given in *P*∗. We have an equal long-range dependence parameter *d*<sup>∗</sup> to regard for all *p* ∈ *P*∗. Applying Lemma A1, we obtain the following, by using the symmetry in *l* = 0 of the cross-correlation function *<sup>r</sup>*(*p*,*q*)(0) = *<sup>r</sup>*(*q*,*p*)(0) for *<sup>p</sup>*, *<sup>q</sup>* <sup>∈</sup> *<sup>P</sup>*∗:

*d* ∑ *p*,*q*=1 \$ *<sup>h</sup>*−<sup>1</sup> ∑ *l*=0 *h*−*l* ∑ *i*=1 *α* (*p*,*q*) *i*,*i*+*l r*ˆ (*p*,*q*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (*l*) + *h*−1 ∑ *l*=1 *h*−*l* ∑ *i*=1 *α* (*p*,*q*) *i*+*l*,*i r*ˆ (*q*,*p*) *<sup>n</sup>* (*l*) <sup>−</sup> *<sup>r</sup>*(*q*,*p*) (*l*) %% = ∑ *p*,*q*∈*P*<sup>∗</sup> \$ *<sup>h</sup>*−<sup>1</sup> ∑ *l*=0 *h*−*l* ∑ *i*=1 *α* (*p*,*q*) *i*,*i*+*l r*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (0) + *h*−1 ∑ *l*=1 *h*−*l* ∑ *i*=1 *α* (*p*,*q*) *i*+*l*,*i r*ˆ (*q*,*p*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*q*,*p*) (0) %% + *o*P(*n*2*d*∗−1) = ∑ *p*,*q*∈*P*<sup>∗</sup> *r*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (0) \$ *<sup>h</sup>*−<sup>1</sup> ∑ *l*=0 *h*−*l* ∑ *i*=1 *α* (*p*,*q*) *<sup>i</sup>*,*i*+*<sup>l</sup>* + *h*−1 ∑ *l*=1 *h*−*l* ∑ *i*=1 *α* (*p*,*q*) *i*+*l*,*i* % + *o*P(*n*2*d*<sup>∗</sup> ) = ∑ *p*,*q*∈*P*<sup>∗</sup> *r*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (0) \$ *<sup>h</sup>* ∑ *i*,*k*=1 *α* (*p*,*q*) *i*,*k* % + *o*P(*n*2*d*<sup>∗</sup> ) = ∑ *p*,*q*∈*P*<sup>∗</sup> *α*˜(*p*,*q*) *r*ˆ (*p*,*q*) *<sup>n</sup>* (0) <sup>−</sup> *<sup>r</sup>*(*p*,*q*) (0) + *o*P(*n*2*d*<sup>∗</sup> ), (A27)

by defining *<sup>α</sup>*˜(*p*,*q*) :<sup>=</sup> *<sup>h</sup>* ∑ *i*,*k*=1 *α* (*p*,*q*) *<sup>i</sup>*,*<sup>k</sup>* . Applying the continuous mapping theorem given in [33], Theorem 2.3, to the result in Corollary A1, we arrive at:

$$\begin{split} n^{-2d^\*} \left( \mathbb{C}\_2 \right)^{-1/2} \sum\_{j=1}^n f \left( Y\_{j,h} \right) &= n^{-2d^\*} \left( n \sum\_{p,q=1}^d \mathbb{A}^{(p,q)} \left( \hat{r}\_n^{(p,q)}(0) - r^{(p,q)}(0) \right) + o\_\mathbb{P} \left( n^{2d^\*} \right) \right) \\ &= n^{1-2d^\*} \left( \mathbb{C}\_2 \right)^{-1/2} \sum\_{p,q=1}^d \mathbb{A}^{(p,q)} \left( \hat{r}\_n^{(p,q)}(0) - r^{(p,q)}(0) \right) + o\_\mathbb{P} (1) \\ &\stackrel{\mathcal{D}}{\rightarrow} \sum\_{p,q \in P^\*} \mathbb{A}^{(p,q)} Z\_{2,d^\*+1/2}^{(p,q)}(1), \end{split}$$

where:

$$Z\_{2,d^\*+1/2}^{(p,q)}(1) = K\_{p,q}(d^\*) \int\_{\mathbb{R}^2}^{\prime\prime} \frac{\exp(i(\lambda\_1 + \lambda\_2)) - 1}{i(\lambda\_1 + \lambda\_2)} |\lambda\_1 \lambda\_2|^{-d^\*} \tilde{B}\_L^{(p)}(\mathbf{d}\lambda\_1) \tilde{B}\_L^{(q)}(\mathbf{d}\lambda\_2).$$

The matrix *K*(*d*∗) is given in Corollary A1. Moreover, *B*˜*L*(d*λ*) is a multivariate Hermitian–Gaussian random measure with E(*BL*(d*λ*)*BL*(d*λ*)∗) = *L* d*λ* and *L* as defined in (2).

#### *Appendix A.3. Proof of Corollary 2.2*

**Proof.** We assumed *d*<sup>∗</sup> ∈ 1 4 , 1 2 , because otherwise we leave the long-range dependent setting, since we are studying functionals with Hermite rank 2 and the transformed process would no longer be long-range dependent and limit theorems for functionals of short-range dependent processes would hold, as can be seen in Theorem 4 in [14]. This choice of *d*∗ assures that the multivariate generalization of the Reduction Theorem as it is used in the proof of Theorem 2.1 still holds for these softened assumptions, as explained in (8).

We turn to the asymptotics of *g*(*p*,*q*) -*Yj* . We obtain for all *p*, *q* ∈ {1, ... , *d*} \ *P*∗, i.e., excluding *dp* = *dq* = *d*<sup>∗</sup> and for all *l* = 0, . . . , *h* − 1 as in (A26), that:

$$\begin{split} \mathbb{E}\left(\left(n^{1-2d^\*}\left(\mathfrak{f}\_{\boldsymbol{n}}^{\left(p,q\right)}(\boldsymbol{l})-r^{\left(p,q\right)}(\boldsymbol{l})\right)\right)^2\right) &= n^{2\left(d\_p+d\_q-2d^\*\right)}\mathbb{E}\left(\left(n^{1-d\_p-d\_q}\left(\mathfrak{f}\_{\boldsymbol{n}}^{\left(p,q\right)}(\boldsymbol{l})-r^{\left(p,q\right)}(\boldsymbol{l})\right)\right)^2\right) \\ &= n^{2d\_p+2d\_q-4d^\*}\mathbb{C}\_2\left(L\_{p,p}L\_{q,q}+L\_{p,q}L\_{q,p}\right) \\ &\xrightarrow{(\mathfrak{n}\to\infty)} 0,\end{split} \tag{A28}$$

since *dp* + *dq* − 2*d*<sup>∗</sup> < 0.

This implies that:

$$n^{1-2d^\*} \left( \mathfrak{h}\_n^{(p,q)}(0) - r^{(p,q)}(0) \right) \xrightarrow{\mathbb{P}} 0.$$

Applying Slutsky's theorem, we observe that only *p*, *q* ∈ *P*<sup>∗</sup> have an impact on the convergence behavior as it is given in (A27) and hence, the result in Theorem 2.1 holds.

#### *Appendix A.4. Proof of Theorem 2.3*

**Proof.** We follow the proof of Theorem 2.1 until (A27), in order to obtain a limit distribution that can be expressed by the sum of two standard Rosenblatt random variables:

$$\begin{split} &\sum\_{p,q=1}^{2}\tilde{\alpha}^{\left(p,q\right)}\left(\hat{\boldsymbol{\boldsymbol{\gamma}}}\_{n}^{\left(p,q\right)}\left(\boldsymbol{0}\right)-r^{\left(p,q\right)}\left(\boldsymbol{0}\right)\right) \\ &=\frac{1}{n}\sum\_{j=1}^{n}\sum\_{p,q=1}^{2}\tilde{\alpha}^{\left(p,q\right)}\left(\boldsymbol{Y}\_{j}^{\left(p\right)}\boldsymbol{Y}\_{j}^{\left(q\right)}-r^{\left(p,q\right)}\left(\boldsymbol{0}\right)\right) \\ &=\frac{1}{n}\sum\_{j=1}^{n}\left(\boldsymbol{Y}\_{j}^{\left(1\right)},\boldsymbol{Y}\_{j}^{\left(2\right)}\right)\left(\begin{matrix}\tilde{\alpha}^{\left(1,1\right)} & \tilde{\alpha}^{\left(1,2\right)}\\ \tilde{\alpha}^{\left(2,1\right)} & \tilde{\alpha}^{\left(2,2\right)}\end{matrix}\right)\left(\boldsymbol{Y}\_{j}^{\left(1\right)},\boldsymbol{Y}\_{j}^{\left(2\right)}\right)^{t} \\ & \qquad \quad -\mathbb{E}\left(\left(\boldsymbol{Y}\_{j}^{\left(1\right)},\boldsymbol{Y}\_{j}^{\left(2\right)}\right)\begin{pmatrix}\tilde{\alpha}^{\left(1,1\right)} & \tilde{\alpha}^{\left(1,2\right)}\\ \tilde{\alpha}^{\left(2,1\right)} & \tilde{\alpha}^{\left(2,2\right)}\end{matrix}\right)\left(\boldsymbol{Y}\_{j}^{\left(1\right)},\boldsymbol{Y}\_{j}^{\left(2\right)}\right)^{t}\right). \tag{A29} \end{split} \tag{A20}$$

We remember that *α*˜(*p*,*q*) = ∑*<sup>h</sup> <sup>i</sup>*,*k*=<sup>1</sup> *α* (*p*,*q*) *<sup>i</sup>*,*<sup>k</sup>* <sup>=</sup> <sup>∑</sup>*<sup>h</sup> <sup>i</sup>*,*k*=<sup>1</sup> *<sup>α</sup>i*+(*p*−1)*h*,*k*+(*q*−1)*<sup>h</sup>* for *<sup>p</sup>*, *<sup>q</sup>* = 1, 2 and <sup>A</sup> <sup>=</sup> (*αi*,*k*)1≤*i*,*k*≤2*<sup>h</sup>* <sup>=</sup> <sup>Σ</sup>−<sup>1</sup> 2,*hC*Σ−<sup>1</sup> 2,*<sup>h</sup>* . Since <sup>Σ</sup>−<sup>1</sup> 2,*<sup>h</sup>* is the inverse of the covariance matrix Σ2,*<sup>h</sup>* of *Y*1,*<sup>h</sup>* it is a symmetric matrix. The matrix of second order Hermite coefficients *C* has the representation *C* = E *Yj*,*<sup>h</sup> f Yj*,*<sup>h</sup> Yt j*,*h* and therefore, *ci*,*<sup>k</sup>* = E *Y*(*i*) *<sup>j</sup>*,*<sup>h</sup> <sup>Y</sup>*(*k*) *<sup>j</sup>*,*<sup>h</sup> f Yj*,*<sup>h</sup>* <sup>=</sup> *ck*,*<sup>i</sup>* for each *i*, *k* = 1, ... , 2*h*. Then, A is also a symmetric matrix, since A*<sup>t</sup>* = Σ−<sup>1</sup> 2,*hC*Σ−<sup>1</sup> 2,*h t* = Σ−<sup>1</sup> 2,*h t Ct* Σ−<sup>1</sup> 2,*h t* <sup>=</sup> <sup>A</sup>. We can now show that \$ *α*˜(1,1) *α*˜(1,2) *α*˜(2,1) *α*˜(2,2) % is a symmetric matrix, i.e., *α*˜(1,2) = *α*˜(2,1). To this end, we define I*<sup>p</sup>* = (0, 0, . . . , 0, 1, . . . , 1, 0, . . . , 0) *<sup>t</sup>* <sup>∈</sup> <sup>R</sup>2*<sup>h</sup>* such that I (*i*) *<sup>p</sup>* = 1 only if *i* = (*p* − 1)*h* + 1, . . . , *ph*, *p* = 1, 2. Then, we obtain:

$$\mathbb{R}^{(1,2)} = \sum\_{i,k=1}^{h} \mathbb{a}^{(1,2)}\_{i,k} = \left(\mathbb{\bar{a}}^{(1,2)}\right)^{t} = \left(\mathbb{I}^{t}\_{1}\mathbb{A}\mathbb{I}\_{2}\right)^{t} = \mathbb{I}^{t}\_{2}\mathbb{A}\mathbb{I}\_{1} = \mathbb{\bar{a}}^{(2,1)}.$$

We now apply the new assumption that *<sup>r</sup>*(1,1)(*l*) = *<sup>r</sup>*(2,2)(*l*), for *<sup>l</sup>* <sup>=</sup> 0, ... , *<sup>h</sup>* <sup>−</sup> 1 and show *α*˜(1,1) = *α*˜(2,2) with the symmetry features of the multivariate normal distribution discussed in (2.2) and in (2.3) in [18], since *ci*,*<sup>j</sup>* = *<sup>c</sup>*2*h*−*i*+1,2*h*−*j*<sup>+</sup>1, *<sup>i</sup>*, *<sup>j</sup>* = 1, . . . , 2*h*. We have to study:

$$\mathfrak{a}^{(2,2)} = \left(\mathbb{I}\_2^t \mathbf{A} \mathbb{I}\_2\right)^\mathbf{f} = \mathbb{I}\_2^t \Sigma\_{2,h}^{-1} \mathbb{C} \Sigma\_{2,h}^{-1} \mathbb{I}\_2.$$

Since Σ−<sup>1</sup> 2,*<sup>h</sup>* <sup>=</sup> (*gi*,*k*)1≤*i*,*k*≤2*<sup>h</sup>* is a symmetric and persymmetric matrix, we have *gi*,*<sup>k</sup>* <sup>=</sup> *gk*,*<sup>i</sup>* and *gi*,*<sup>k</sup>* = *<sup>g</sup>*2*h*−*i*+1,2*h*−*k*+<sup>1</sup> for *<sup>i</sup>*, *<sup>k</sup>* = 1, . . . , 2*h*. Then, we obtain:

$$\begin{split} \Pi\_{2}^{t}\Sigma\_{2,\mathbb{H}}^{-1} &= \left(\sum\_{i=k+1}^{2h} \mathcal{G}\_{i,1}, \ldots, \sum\_{i=h+1}^{2h} \mathcal{G}\_{i,2h}\right) \\ &= \left(\sum\_{i=1}^{h} \mathcal{G}\_{i+h,1}, \ldots, \sum\_{i=1}^{h} \mathcal{G}\_{i+h,2h}\right) \\ &= \left(\sum\_{i=1}^{h} \mathcal{G}\_{h-i+1,2h}, \ldots, \sum\_{i=1}^{h} \mathcal{G}\_{h-i+1,1}\right) \\ &= \left(\sum\_{i=1}^{h} \mathcal{G}\_{i,2h}, \ldots, \sum\_{i=1}^{h} \mathcal{G}\_{i,1}\right) \\ &= \left(\sum\_{i=1}^{h} \mathcal{G}\_{2h,i}, \ldots, \sum\_{i=1}^{h} \mathcal{G}\_{1,i}\right) \\ &=: \left(\tilde{\mathcal{G}}\_{2h}, \ldots, \tilde{\mathcal{G}}\_{1}\right). \end{split}$$

Note that:

$$\Sigma\_{2,h}^{-1}\mathbb{I}\_1 = \left(\sum\_{i=1}^h \mathbb{g}\_{1,i}, \dots, \sum\_{i=1}^h \mathbb{g}\_{2h,i}\right)^t = \left(\mathbb{g}\_1, \dots, \mathbb{g}\_{2h}\right)^t.$$

Then, we arrive at:

$$\begin{split} \bar{\mathbb{A}}^{(2,2)} &= \left(\mathbb{I}\_{2}^{t}\mathbb{A}\mathbb{I}\_{2}\right)^{t} = \mathbb{I}\_{2}^{t}\Sigma\_{2,l}^{-1}\mathbb{C}\Sigma\_{2,h}^{-1}\mathbb{I}\_{2} \\ &= \sum\_{i,k=1}^{2h} \tilde{\xi}2^{h-i+1}\tilde{\xi}2^{h-k+1}\mathbb{C}i\_{k} \\ &= \sum\_{i,k=1}^{2h} \tilde{\xi}2^{h-i+1}\tilde{\xi}2^{h-k+1}\mathbb{C}2^{h-i+1}2^{h-k+1} \\ &= \sum\_{i,k=1}^{2h} \tilde{\xi}i\tilde{\xi}k^{c}\mathbb{C}i\_{k} \\ &= \mathbb{I}\_{1}^{t}\Sigma\_{2,h}^{-1}\mathbb{C}\Sigma\_{2,h}^{-1}\mathbb{I}\_{1} \\ &= \bar{\mathbb{A}}^{(1,1)}. \end{split}$$

Therefore, we have to deal with a special type of 2 × 2 matrix, since the original matrix in the formula (A27), namely \$ *α*˜(1,1) *α*˜(1,2) *α*˜(2,1) *α*˜(2,2) % , has now reduced to \$ *α*˜(1,1) *α*˜(1,2) *α*˜(1,2) *α*˜(1,1) % .

Finally, we know that any real-valued symmetric matrix *A* can be decomposed via diagonalization into an orthogonal matrix *V* and a diagonal matrix *D*, where the entries of the latter one are determined via the eigenvalues of *A*, for details, as can be seen in [34], p. 327.

We can explicitly give formulas for the entries of these matrices here:

$$V = \begin{pmatrix} -2^{-1/2} & 2^{-1/2} \\ 2^{-1/2} & 2^{-1/2} \end{pmatrix}, \quad D = \begin{pmatrix} \lambda\_1 = \mathbb{\bar{n}}^{(1,1)} - \mathbb{\bar{n}}^{(1,2)} & 0 \\ 0 & \lambda\_2 = \mathbb{\bar{n}}^{(1,1)} + \mathbb{\bar{n}}^{(1,2)} \end{pmatrix},$$

such that:

$$VDV = \begin{pmatrix} \frac{\lambda\_1 + \lambda\_2}{2} & \frac{\lambda\_2 - \lambda\_1}{2} \\ \frac{\lambda\_2 - \lambda\_1}{2} & \frac{\lambda\_1 + \lambda\_2}{2} \end{pmatrix} = \begin{pmatrix} \tilde{a}^{(1,1)} & \tilde{a}^{(1,2)} \\ \tilde{a}^{(1,2)} & \tilde{a}^{(1,1)} \end{pmatrix}.$$

Therefore, continuing with (A29), we now have the representation:

1 *n n* ∑ *j*=1 *Y*(1) *<sup>j</sup>* ,*Y*(2) *j* \$ *α*˜(1,1) *α*˜(1,2) *α*˜(2,1) *α*˜(2,2) % *Y*(1) *<sup>j</sup>* ,*Y*(2) *j t* <sup>−</sup> <sup>E</sup> \$ *Y*(1) *<sup>j</sup>* ,*Y*(2) *j* \$ *α*˜(1,1) *α*˜(1,2) *α*˜(2,1) *α*˜(2,2) % *Y*(1) *<sup>j</sup>* ,*Y*(2) *j t* % = 1 *n n* ∑ *j*=1 *Y*(1) *<sup>j</sup>* ,*Y*(2) *j VDV Y*(1) *<sup>j</sup>* ,*Y*(2) *j t* <sup>−</sup> <sup>E</sup> *Y*(1) *<sup>j</sup>* ,*Y*(2) *j VDV Y*(1) *<sup>j</sup>* ,*Y*(2) *j t* = 1 *n n* ∑ *j*=1 *<sup>α</sup>*˜(1,1) <sup>−</sup> *<sup>α</sup>*˜(1,2) 2 *Y*(2) *<sup>j</sup>* <sup>−</sup> *<sup>Y</sup>*(1) *j* 2 <sup>−</sup> <sup>E</sup> *Y*(2) *<sup>j</sup>* <sup>−</sup> *<sup>Y</sup>*(1) *j* 2 + 1 *n n* ∑ *j*=1 *α*˜(1,1) + *α*˜(1,2) 2 *Y*(1) *<sup>j</sup>* <sup>+</sup> *<sup>Y</sup>*(2) *j* 2 <sup>−</sup> <sup>E</sup> *Y*(1) *<sup>j</sup>* <sup>+</sup> *<sup>Y</sup>*(2) *j* 2 = 1 *n n* ∑ *j*=1 *<sup>α</sup>*˜(1,1) <sup>−</sup> *<sup>α</sup>*˜(1,2) <sup>1</sup> <sup>−</sup> *<sup>r</sup>*(1,2) (0) ⎛ ⎜⎝ ⎛ ⎝ *Y*(2) *<sup>j</sup>* <sup>−</sup> *<sup>Y</sup>*(1) *j* 0 <sup>2</sup> <sup>−</sup> <sup>2</sup>*r*(1,2)(0) ⎞ ⎠ 2 − 1 ⎞ ⎟⎠ + 1 *n n* ∑ *j*=1 *α*˜(1,1) + *α*˜(1,2) <sup>1</sup> <sup>+</sup> *<sup>r</sup>*(1,2) (0) ⎛ ⎜⎝ ⎛ ⎝ *Y*(1) *<sup>j</sup>* <sup>+</sup> *<sup>Y</sup>*(2) *j* 0 2 + 2*r*(1,2)(0) ⎞ ⎠ 2 − 1 ⎞ ⎟⎠ <sup>=</sup> <sup>1</sup> *n <sup>α</sup>*˜(1,1) <sup>−</sup> *<sup>α</sup>*˜(1,2) <sup>1</sup> <sup>−</sup> *<sup>r</sup>*(1,2) (0) *<sup>n</sup>* ∑ *j*=1 *H*<sup>2</sup> *Y*∗ *j* + 1 *n α*˜(1,1) + *α*˜(1,2) <sup>1</sup> <sup>+</sup> *<sup>r</sup>*(1,2) (0) *<sup>n</sup>* ∑ *j*=1 *H*<sup>2</sup> *Y*∗∗ *j* , (A30)

with *Y*∗ *<sup>j</sup>* :<sup>=</sup> *<sup>Y</sup>*(2) *<sup>j</sup>* <sup>−</sup>*Y*(1) √ *j* <sup>2</sup>−2*r*(1,2)(0) and *Y*∗∗ *<sup>j</sup>* :<sup>=</sup> *<sup>Y</sup>*(1) *<sup>j</sup>* <sup>+</sup>*Y*(2) √ *j* 2+2*r*(1,2)(0)

Now, note that:

$$\mathbb{E}\left(Y\_{\vec{j}}^{\*}Y\_{\vec{j}}^{\*\*}\right) = \mathbb{E}\left(\frac{Y\_{\vec{j}}^{(2)} - Y\_{\vec{j}}^{(1)}}{\sqrt{2 - 2r^{(1,2)}(0)}} \frac{Y\_{\vec{j}}^{(1)} + Y\_{\vec{j}}^{(2)}}{\sqrt{2 + 2r^{(1,2)}(0)}}\right) = 0.1$$

.

Therefore, we created a bivariate long-range dependent Gaussian process, since:

$$
\begin{pmatrix} -1 & 1 \\ 1 & -1 \end{pmatrix} \begin{pmatrix} \chi^{(1)}\_{\dot{\jmath}}, \chi^{(2)}\_{\dot{\jmath}} \end{pmatrix}^{t} = \begin{pmatrix} \chi^\*\_{\dot{\jmath}}, \chi^{\*\*}\_{\dot{\jmath}} \end{pmatrix}^{t} \sim \mathcal{N}(0, I\_2),
$$

with cross-covariance function:

$$\begin{split} r\_{\*}^{(1,2)}(k) := \mathbb{E} \left( Y\_{j}^{\*} Y\_{j+k}^{\*\ast} \right) &= \mathbb{E} \left( \frac{Y\_{j}^{(2)} - Y\_{j}^{(1)}}{\sqrt{2 - 2r^{(1,2)}(0)}} \frac{Y\_{j+k}^{(1)} + Y\_{j+k}^{(2)}}{\sqrt{2 + 2r^{(1,2)}(0)}} \right) \\ &= \frac{r^{(2,1)}(k) + r^{(2,2)}(k) - r^{(1,1)}(k) - r^{(1,2)}(k)}{2\sqrt{\left(1 - r^{(1,2)}(0)\right) \left(1 + r^{(1,2)}(0)\right)}} \\ &\simeq \frac{L\_{2,2} + L\_{2,1} - L\_{1,2} - L\_{1,1}}{2\sqrt{\left(1 - r^{(1,2)}(0)\right) \left(1 + r^{(1,2)}(0)\right)}} k^{2d^{\ast} - 1} . \end{split} \tag{A31}$$

Note that the covariance functions have the following asymptotic behavior:

$$\begin{split} r^{(1,1)}\_{\*}(k) := \mathbb{E}\left(Y^{\*}\_{j}Y^{\*}\_{j+k}\right) &= \mathbb{E}\left(\frac{Y^{(2)}\_{j} - Y^{(1)}\_{j}}{\sqrt{2 - 2r^{(1,2)}(0)}} \frac{Y^{(2)}\_{j+k} - Y^{(1)}\_{j+k}}{\sqrt{2 - 2r^{(1,2)}(0)}}\right) \\ &= \frac{r^{(2,2)}(k) - r^{(2,1)}(k) - r^{(1,2)}(k) + r^{(1,1)}(k)}{2 - 2r^{(1,2)}(0)} \\ &\simeq \underbrace{\frac{L\_{2,2} - L\_{2,1} - L\_{1,2} + L\_{1,1}}{2 - 2r^{(1,2)}(0)}}\_{-:L^{\*}\_{1,1}} k^{2d^{\*} - 1} \end{split}$$

and analogously:

$$r\_\*^{(2,2)}(k) := \mathbb{E}\left(\mathbf{Y}\_j^{\*\*}\mathbf{Y}\_{j+k}^{\*\*}\right) \simeq \underbrace{\frac{L\_{2,2} + L\_{2,1} + L\_{1,2} + L\_{1,1}}{2 + 2r^{(1,2)}(0)}}\_{=: L\_{2,2}^s} k^{2d^\*-1}.$$

We can now apply the result of [14], Theorem 6, since we created a bivariate Gaussian process with independent entries for fixed *j*. Note that for the function we apply here, namely ˜ *f Y*∗ *<sup>j</sup>* ,*Y*∗∗ *j* = *H*<sup>2</sup> *Y*∗ *j* + *H*<sup>2</sup> *Y*∗∗ *j* the weighting factors in [14], Theorem 6, reduce to *e*1,1 = *e*2,2 = 1 and *e*1,2 = *e*2,1 = 0. These weighting factors fit into the result in [14], (3.6) and (3.7), which even yields the joint convergence of the vector of both univariate summands, *H*<sup>2</sup> *Y*∗ *j* , *H*<sup>2</sup> *Y*∗∗ *j* , suitably normalized to a vector of two (dependent) Rosenblatt random variables. Since the long-range dependence property in [13], Definition 2.1 is more specific than in [14], p. 2259, (3.1) (see the considerations in (8)), we are able to scale the variances of each Rosenblatt random variable to 1 and give the covariance between them, by using the normalization given in [15], Theorem 4.3. We obtain:

$$\begin{split} &n^{-2d\*}(2\mathsf{C}\_{2})^{-1/2}\Big(\tilde{\mathfrak{a}}^{(1,1)}-\tilde{\mathfrak{a}}^{(1,2)}\Big)\Big(1-r^{(1,2)}(0)\Big)\sum\_{j=1}^{n}H\_{2}\Big(Y\_{j}^{+}\Big) \\ &+n^{-2d\*}(2\mathsf{C}\_{2})^{-1/2}\Big(\tilde{\mathfrak{a}}^{(1,1)}+\tilde{\mathfrak{a}}^{(1,2)}\Big)\Big(1+r^{(1,2)}(0)\Big)\sum\_{j=1}^{n}H\_{2}\Big(Y\_{j}^{+\*}\Big) \\ \overset{\mathcal{D}}{\rightarrow}&\Big(\tilde{\mathfrak{a}}^{(1,1)}-\tilde{\mathfrak{a}}^{(1,2)}\Big)\Big(1-r^{(1,2)}(0)\Big)L\_{1,1}^{\*}Z\_{2,d^{\*}+1/2}^{\*}(1) \\ &+\Big(\tilde{\mathfrak{a}}^{(1,1)}+\tilde{\mathfrak{a}}^{(1,2)}\Big)\Big(1+r^{(1,2)}(0)\Big)L\_{2,2}^{\*}Z\_{2,d^{\*}+1/2}^{\*\*} \\ = &\Big(\tilde{\mathfrak{a}}^{(1,1)}-\tilde{\mathfrak{a}}^{(1,2)}\Big)\frac{L\_{2,2}-L\_{2,1}-L\_{1,2}+L\_{1,1}}{2}L\_{2,d^{\*}+1/2}^{\*}(1) \\ &+\Big(\tilde{\mathfrak{a}}^{(1,1)}+\tilde{\mathfrak{a}}^{(1,2)}\Big)\frac{L\_{2,2}+L\_{2,1}+L\_{1,2}+L\_{1,1}}{2}L\_{2,d^{\*}+1/2}^{\*}(1) \end{split}$$

with *C*<sup>2</sup> := <sup>1</sup> <sup>2</sup>*d*∗(4*d*∗−1) being the same normalizing factor as in Theorem 2.1.

We observe that *Z*∗ 2,*d*∗+1/2(1) and *Z*∗∗ 2,*d*∗+1/2(1) are both standard Rosenblatt random variables. Following Corollary A1, their covariance is given by

$$\begin{split} \text{Cov}\left(Z^{\*}\_{2,d^{\*}+1/2}(1), Z^{\* \* \*}\_{2,d^{\*}+1/2}(1)\right) &= \frac{\left(L^{\*}\_{1,2} + L^{\*}\_{2,1}\right)^{2}}{L^{\*}\_{1,1}L^{\*}\_{2,2}} \\ &= \frac{2\left(\left(L\_{2,2} - L\_{1,1}\right)^{2}\right)}{4\left(1 - r^{(1,2)}(0)\right)\left(1 + r^{(1,2)}(0)\right)} \left(L^{\*}\_{1,1}L^{\*}\_{2,2}\right)^{-1} \\ &= \frac{\left(L\_{2,2} - L\_{1,1}\right)^{2}}{\left(L\_{1,1} + L\_{2,2}\right)^{2} - \left(L\_{1,2} + L\_{2,1}\right)^{2}}. \end{split}$$

Note that (*L*1,1 + *L*2,2) <sup>2</sup> <sup>−</sup> (*L*1,2 <sup>+</sup> *<sup>L</sup>*2,1) <sup>2</sup> <sup>=</sup> 0 is fulfilled since *<sup>L</sup>*1,1 <sup>+</sup> *<sup>L</sup>*2,2 <sup>=</sup> *<sup>L</sup>*1,2 <sup>+</sup> *<sup>L</sup>*2,1.

#### **References**

