[Testing equality of spectral density operators for functional linear processes]{}

Anne Leucht\
Technische Universität Braunschweig\
Institut für Mathematische Stochastik\
Universitätsplatz 2, D-38106 Braunschweig, Germany\
E-mail: a.leucht@tu-braunschweig.de\
Efstathios Paparoditis\
University of Cyprus\
Department of Mathematics and Statistics\
P.O.Box 20537, CY-1678, Nicosia, Cyprus\
E-mail: stathisp@ucy.ac.cy\
Theofanis Sapatinas\
University of Cyprus\
Department of Mathematics and Statistics\
P.O.Box 20537, CY-1678, Nicosia, Cyprus\
E-mail: fanis@ucy.ac.cy

[**Abstract**]{}

The problem of testing equality of the entire second order structure of two independent functional linear processes is considered. A fully functional $L^2$-type test is developed which evaluates, over all frequencies, the Hilbert-Schmidt distance between the estimated spectral density operators of the two processes. The asymptotic behavior of the test statistic is investigated and its limiting distribution under the null hypothesis is derived. Furthermore, a novel frequency domain bootstrap method is developed which approximates more accurately the distribution of the test statistic under the null than the large sample Gaussian approximation obtained. Asymptotic validity of the bootstrap procedure is established and consistency of the bootstrap-based test under the alternative is proved. Numerical simulations show that, even for small samples, the bootstrap-based test has very good size and power behavior. An application to meteorological functional time series is also presented.

[*2010 Mathematics Subject Classification*]{}: Primary 62M10; Secondary 62M15, 60G10\
\
[*Keywords and Phrases*]{}: Bootstrap, Functional Time Series, Hypothesis Testing, Spectral Density Operator\
\
[*Short title*]{}: Testing Equality of Spectral Density Operatorsversion:

Introduction {#S1}
============

Functional time series form a natural class of functional data occurred in many applications such as daily curves of financial transactions, daily images of geophysical and environmental data and daily curves of temperature measurements. Such curves or images are viewed as functions in appropriate spaces since an observed intensity is available at each point on a line segment, a portion of a plane or a volume. Moreover, and most importantly, such functional time series exhibit temporal dependence while ignoring this dependence may result in misleading conclusions and inferential procedures.

Comparing characteristics of two or more groups of functional data forms an important problem of statistical inference with a variety of applications. For instance, comparing the mean functions between independent groups of independent and identically distributed (i.i.d.) functional data or of functional time series has attracted considerable interest in the literature, see, e.g., Benko [*et al.*]{} (2009), Zhang [*et al.*]{} (2010), Horváth and Kokoszka (2012, Chapter 5), Horváth [*et al.*]{} (2013) and Paparoditis and Sapatinas (2016). In contrast to comparing mean functions, the problem of comparing the entire second order structure of two independent functional time series has been much less investigated. Notice that for i.i.d. functional data this problem simplifies to the problem of testing the equality of (the lag zero) covariance operators, see, e.g., Panaretos [*et al.*]{} (2010), Fremdt [*et al.*]{} (2012) (finite-dimensional projections), Pigoli [*et al.*]{} (2014) (distance measures) and Paparoditis and Sapatinas (2016) (fully functional). The same problem of testing the equality of the (lag-zero) covariance operators of two sets of independent functional time series has also been investigated by Zhang and Shao (2015) using a finite-dimensional projections approach.

However, the comparison of the entire second order structure of independent functional time series, is a much more involved problem due to the temporal dependence of the set of random elements considered. In describing the second order structure of functional time series, the spectral density operator, introduced in the functional set-up by Panaretos and Tavakoli (2013), is a very useful tool since it summarizes in a nice way the entire autocovariance structure of the underlying functional time series; see also Hörmann [*et al.*]{} (2015) and van Delft and Eichler (2018). It seems, therefore, very appealing to develop a spectral approach for testing equality of the second order properties of two functional time series. A testing procedure into this direction has been recently proposed by Tavakoli and Panaretos (2016) where projections on finite dimensional spaces of the differences of the estimated spectral density operators of the two functional time series have been used. Although projection-based tests have the advantage to lead to manageable limiting distributions, and can be very powerful when the deviations from the null are captured by the finite-dimensional space projected, such tests have no power for alternatives which are orthogonal to the projection space. It is, therefore, important to develop a fully functional test for the testing problem at hand. Notice that in the finite-dimensional case, i.e., for (univariate or multivariate) real-valued time series, such tests have been developed among others by Eichler (2008) and Dette and Paparoditis (2009).

The contribution of this paper is twofold. First, we develop a fully functional approach for testing equality of the entire second order structure between two independent functional linear processes. The testing procedure proposed, evaluates, for each frequency, the Hilbert-Schmidt norm between the (estimated) spectral density operators of the functional linear process at hand. Integrating these differences over all possible frequencies, leads to a global, $L_2$-type, measure of deviation which is used to test the null hypothesis of interest. We derive the limiting distribution of an appropriately centered version of such a test statistic under the null. Second, and because of the slow convergence of the distribution of the test statistic under the null against its derived limiting Gaussian distribution, we develop a novel frequency domain bootstrap procedure to estimate this distribution. The method works under quite general conditions on the underlying functional linear process and its range of applicability is not restricted to the particular test statistic considered. We show that the bootstrap procedure approximates correctly the distribution of the proposed test statistic under the null. Furthermore, consistency of the bootstrap-based test under the alternative is established. Our theoretical deviations are accomplished by a simulation study which shows a very good behavior of the bootstrap procedure in approximating the distribution of interest and that the test based on the bootstrap critical values has good size and power behavior. Furthermore, a real-life data application is presented which together with some test-diagnostics enables an interesting discussion of the test results obtained.

The remaining of the paper is organized as follows. Section 2 contains the main assumptions on the underlying functional linear processes and states the hypothesis testing problem under study. Section 3 is devoted to the suggested test statistic and its asymptotic behavior while Section 4 discusses a frequency domain bootstrap procedure to estimate the distribution of the test statistic under the null. Asymptotic validity of the bootstrap procedure is established and consistency of the bootstrap-based test under the alternative is also proved. Section 5 contains numerical simulations and an application to meteorological functional time series. Auxiliary results and proofs of the main results are deferred to the Appendix.

Assumptions and the Testing Problem {#S2}
===================================

Suppose that observations $X_1,\dots, X_T$ and $Y_1,\dots, Y_T$ stem from two independent functional linear processes $(X_t)_{t\in\Z}$ and $(Y_t)_{t\in\Z}$ satisfying the following assumption.

[*Assumption 1*]{} : $(X_t)_{t\in\Z}$ and $(Y_t)_{t\in\Z}$ are independent functional linear processes, given by $$\label{eq.lin-proc}
X_t=\sum_{j\in \Z} a_j \varepsilon_{t-j}\quad\text{and}\quad Y_t=\sum_{j\in \Z} b_j e_{t-j},\; t\in \Z,$$ with values in $L^2_\R([0,1], \mu)$, where $\mu$ denotes the Lebesgue measure. The innovations $(\varepsilon_t)_{t \in \Z}$ and $(e_t)_{t \in \Z}$ are two i.i.d. mean zero Gaussian processes with values in $L^2_\R([0,1], \mu)$. The (real-valued) coefficients $(a_j)_{j \in \Z}$ and $(b_j)_{j \in \Z}$ satisfy $\sum_{j \in \Z}|j|^{1/2}(|a_j|+|b_j|)<\infty$.

We are interested in testing for equality of the second order structure between the two functional linear processes given in (\[eq.lin-proc\]). It turns out that a spectral approach to this problem is very appealing. First, note that we can define a spectral density operator in the sense of Panaretos and Tavakoli (2013) in the present set up generalizing the concept of spectral densities for univariate time series and spectral density matrices for multivariate time series. Here and in the sequel, we will abbreviate $L^2_\R([0,1]^d,\mu)$ by $L^2$ if the dimension $d$ becomes clear from the context.

\[l.ex-SDO\] Suppose that $(X_t)_{t\in\Z}$ and $(Y_t)_{t\in \Z}$ are functional linear processes satisfying Assumption 1. Then, for arbitrary $\lambda\in (-\pi,\pi]$, $$f_{X,\lambda}(\cdot,\cdot)=\frac{1}{2\pi}\sum_{t\in\Z} \exp(-i\lambda t) \, r_{X,t}(\cdot,\cdot)\quad\text{and}\quad f_{Y,\lambda}(\cdot,\cdot)=\frac{1}{2\pi}\sum_{t\in\Z} \exp(-i\lambda t)\, r_{Y,t}(\cdot,\cdot)$$ with $r_{X,t}$ and $r_{Y,t}$ denoting the autocovariance kernels of $X$ and $Y$ at lag $t$, respectively, converge absolutely in $L^2$. Moreover, for all $\sigma,\tau \in [0,1]$, $$r_{X,t}(\sigma, \tau)=\int_{(-\pi,\pi]} f_{X,\lambda}(\sigma,\tau)\, e^{it\lambda}\, d\lambda\quad\text{and}\quad
r_{Y,t}(\sigma, \tau)=\int_{(-\pi,\pi]} f_{Y,\lambda} (\sigma, \tau)\, e^{it\lambda}\, d\lambda\quad
\forall t \in \mathbb{Z},$$ where equality holds in $L^2$. The operators ${\mathcal F_{X,\lambda}}$ and ${\mathcal F_{Y,\lambda}}$, induced by right integration of $f_{X,\lambda}$ and $f_{Y,\lambda}$, are self-adjoint, nonnegative definite and it holds $$\mathcal F_{X,\lambda}=\frac{1}{2\pi}\sum_{t\in\Z} e^{-i\lambda t} \mathcal R_{X,t}\quad\text{and}\quad\mathcal F_{Y,\lambda}=\frac{1}{2\pi}\sum_{t\in\Z} e^{-i\lambda t} \mathcal R_{Y,t},$$ where $\mathcal R_{X,t}$ and $\mathcal R_{Y,t}$ denote the autocovariance operators of $X$ and $Y$ at lag $t$, induced by right integration of $r_{X,t}$ and $r_{Y,t}$, respectively. Convergence holds in nuclear norm.

As it can be seen from the proof of Lemma \[l.ex-SDO\], it suffices to assume 4th order integrability of the innovations $(\varepsilon_t)_{t \in \Z}$ and $(e_t)_{t \in \Z}$ rather than Gaussianity here. The kernels $f_{X,\lambda}$ and $f_{Y,\lambda}$ are called the *spectral density kernels* (at frequency $\lambda$) and the operators $\mathcal F_{X,\lambda}$ and $\mathcal F_{Y, \lambda}$ are referred to as the corresponding [*spectral density operators*]{}.

Under the assumptions of Lemma \[l.ex-SDO\], we can now state the hypothesis testing problem of interest as follows $$\begin{aligned}
&{\mathcal H_0}\colon \mathcal F_{X,\lambda}=\mathcal F_{Y, \lambda} \quad\text{ for }\mu\text{-}\text{almost all } \lambda\in (-\pi,\pi], 
 \quad\text{versus}\\
&{\mathcal H_1}\colon \mathcal F_{X, \lambda}\neq\mathcal F_{Y, \lambda} \quad \forall\lambda\in A \text{ for some $A$ with }\mu(A)>0.
\end{aligned}
\label{eq:testM}$$

The Test Statistic and its Asymptotic Behavior {#S3}
==============================================

We first estimate the unknown spectral density operator $\mathcal F_{X, \lambda}$ by an integral operator $\widehat{\mathcal F}_{X,\lambda}$ induced by right-integration with the kernel $$\hat f_{X,\lambda}(\sigma,\tau)=\frac{1}{b T}\sum_{t=-N}^{N} W\left(\frac{\lambda-\lambda_t}{b}\right)\, \widehat p_{X, \lambda_t} (\sigma, \tau), \;\; \mbox{for all} \;\; \sigma,\tau \in [0,1],$$ and, similarly, $\mathcal F_{Y, \lambda}$ by an integral operator $\widehat{\mathcal F}_{Y,\lambda}$ induced by right-integration with the kernel $$\hat f_{Y,\lambda}(\sigma,\tau)=\frac{1}{b T}\sum_{t=-N}^{N} W\left(\frac{\lambda-\lambda_t}{b}\right)\, \widehat p_{Y, \lambda_t} (\sigma, \tau),  \;\; \mbox{for all} \;\; \sigma,\tau \in [0,1].$$ Here, $N=[(T-1)/2]$ and $\lambda_t=2\pi t/T, \;t=-N,\dots, N,$ denote the Fourier frequencies. Further $b=b_T>0$ is an asymptotically vanishing bandwidth and $W$ denotes a weight function. Moreover, as in Panaretos and Tavakoli (2013), $$\widehat p_{X, \lambda}(\sigma,\tau)=\frac{1}{2\pi T}\sum_{s_1,s_2=1}^T X_{s_1}(\sigma)X_{s_2}(\tau)\, \exp(i\lambda( s_1- s_2)),  \;\; \mbox{for all} \;\; \sigma,\tau \in [0,1],$$ and $$\widehat p_{Y, \lambda}(\sigma,\tau)=\frac{1}{2\pi T}\sum_{s_1,s_2=1}^T Y_{s_1}(\sigma)Y_{s_2}(\tau)\, \exp(i\lambda( s_1- s_2)),  \;\; \mbox{for all} \;\; \sigma,\tau \in [0,1],$$ denote the periodogram kernels based on $X_1,\dots, X_T$ and $Y_1,\dots, Y_T$, respectively.

For the hypothesis testing problem (\[eq:testM\]), we propose the following test statistic which evaluates the distance between the estimated spectral density operators, i.e., $${\mathcal U_T}=\int_{-\pi}^\pi\|\widehat{\mathcal F}_{X,\lambda}-\widehat{\mathcal F}_{Y,\lambda}\|_{HS}^2\, d\lambda,
\label{eq:test-stat-main}$$ where $\|\cdot\|_{HS}$ denotes the Hilbert Schmidt norm on $L^2$. The following theorem shows asymptotic normality of the suitably normalized test statistic $\mathcal U_T$ under ${\mathcal H_0}$.

\[t.asy-test\] Suppose that the stretches of observations $X_1,\dots, X_T$ and $Y_1,\dots, Y_T$ stem from the two independent functional linear processes $(X_t)_{t\in\Z}$ and $(Y_t)_{t\in \Z}$, respectively, satisfying Assumption 1. Moreover, assume that

-   $b\sim T^{-\nu }$ for some $\nu\in (1/4,1/2)$,

-   $W$ is bounded, symmetric, positive, and Lipschitz continuous, has bounded support on $(-\pi,\pi]$ and satisfies $\int_{-\pi}^{\pi} W(x)\, dx=2\pi$.

Then, under ${\mathcal H_0}$, $$\sqrt{b} T \, \mathcal U_T-b^{-1/2} \mu_0\stackrel{d}{\longrightarrow} Z\sim{\mathcal N}(0,\theta_0^2),$$ where $$\begin{aligned}
\mu_0&=\frac{1}{\pi}\, \int_{-\pi}^\pi\left\{\int_0^1 f_{X,\lambda}(\tau,\tau) \,d\tau\,\right\}^2 d\lambda\int_{-\pi}^\pi W^2\left(u\right)\, du\\
\theta_0^2&=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \left\{\int_0^1\int_0^1f_{X,\lambda}^2(\sigma, \tau)\,\,d\sigma\, d\tau\right\}^2\, d\lambda.
\end{aligned}$$

Note that the assumptions (i) and (ii) on the weight function $W$ and the bandwidth $b$, respectively, in Theorem \[t.asy-test\] are identical to the assumptions for multivariate time series used in Dette and Paparoditis (2009).

[A careful inspection of the proof of Theorem \[t.asy-test\] shows that the assumption of Gaussianity on the functional innovations $(\varepsilon_t)_{t \in \Z}$ and $(e_t)_{t \in \Z}$ in (\[eq.lin-proc\]) is not essential. In fact, this assumption is solely used to simplify somehow the technical arguments applied in proving asymptotic normality of finite-dimensional distributions associated with the proof of Theorem \[t.asy-test\]. Moreover, this assumption is not required in order to prove convergence of the mean and of the variance of $ \sqrt{b}T \,\mathcal U_T$ or to prove tightness. Consequently, this assumption can be replaced by other assumptions on the stochastic properties of the innovations $(\varepsilon_t)_{t \in \Z}$ and $(e_t)_{t \in \Z}$, which will allow for the use of different technical arguments, for instance arguments based on the convergence of cumulants to the appropriate limits, in order to establish asymptotic normality of finite-dimensional distributions (fidis).]{}

Based on Theorem \[t.asy-test\], the procedure to test hypothesis (\[eq:testM\]) is then defined as follows: Reject ${\mathcal H_0}$ if and only if $$\label{eq.stud-test}
t_{\mathcal U} = \frac{\sqrt{b} T \, \mathcal U_T-b^{-1/2} \widehat{\mu}_0}{\widehat{\theta}_0} \geq z_{1-\alpha},$$ where $ z_{1-\alpha}$ is the upper $1-\alpha$ percentage point of the standard Gaussian distribution and $ \widehat{\mu}_0$ and $\widehat{\theta}_0$ are consistent estimators of $ \mu_0$ and $\theta_0$, respectively. Such estimators can be, for instance, obtained if the unknown spectral density kernel $ f_{X,\lambda}$ is replaced by the pooled estimator $ \widehat{f}_{\lambda}(\tau,\sigma) =\widehat{f}_{X,\lambda}(\tau,\sigma)/2 +  \widehat{f}_{Y,\lambda}(\tau,\sigma)/2$. Notice that, under ${\mathcal H_0}$, $ f_{X,\lambda}=f_{Y,\lambda}=f_{X,\lambda} /2+f_{Y,\lambda} /2$, that is (asymptotically), it makes no difference if $ f_{X,\lambda}$ in $\mu_0$ and $\theta_0$ is replaced by $ \widehat{f}_{X,\lambda}$ (or by $ \widehat{f}_{Y,\lambda}$) instead of the pooled estimator $\widehat{f}_\lambda$. However, under ${\mathcal H_1}$ it matters and, for this reason, we use the pooled estimator $\widehat{f}_\lambda(\tau,\sigma)$ in applying the studentized test statistic $t_{\mathcal U}$ defined in (\[eq.stud-test\]); see also Lemma \[le.power\] in Section 4. Under the assumption that the pooled estimator $\widehat{f}_\lambda$ is uniformly consistent, (see also Assumption 2 below), it is easily seen that, under ${\mathcal H_0}$, $$t_{\mathcal U}  =\frac{\sqrt{b} T \, \mathcal U_T-b^{-1/2} \mu_0}{\theta_0} +o_P(1),$$ i.e., Theorem \[t.asy-test\] implies that the studentized test $t_{\mathcal U} $ is an asymptotically $\alpha$-level test under ${\mathcal H_0}$, for any desired level $\alpha \in (0,1)$.

Bootstrapping The Test Statistic {#sec.boot}
================================

A problem in implementing the above test occurs from the well-known fact that, even in the finite-dimensional case, the convergence of such $L^2$-norm based test statistics towards their limiting distribution is very slow; see, e.g., Härdle and Mammen (1993), Paparoditis (2000) and Dette and Paparoditis (2009). In this case, bootstrap-based approaches may be very effective. In the following, we develop a bootstrap procedure to estimate the distribution of the test statistic $\mathcal U_T$ defined in (\[eq:test-stat-main\]) (and, hence, of the studentized test $ t_{\mathcal U}$ defined in (\[eq.stud-test\])) under ${\mathcal H_0}$, and we prove its asymptotic validity.

We begin by recalling the fact that for any $k\in \N$ and any set of points $ 0\leq s_1 < s_2 < \ldots < s_k\leq 1$ in the interval $[0,1]$, the corresponding $k$-dimensional vector of finite Fourier transforms $$J_{X,\lambda}=\Big(J_{X,\lambda}(s_j) = (2\pi T)^{-1/2}\sum_{t=1}^{T} X_t(s_j) e^{-i t\lambda},\ j=1,2, \ldots, k\Big),$$ satisfies $$\label{eq.fft1}
\left( \begin{array}{c} J_{X,\lambda}(s_1)\\J_{X,\lambda}(s_2)\\ \vdots\\ J_{X,\lambda}(s_k) \end{array}\right) 
\stackrel{d}{\rightarrow} \mathcal N_C\Big( \left( \begin{array}{c} 0\\0\\ \vdots \\ 0 \end{array}\right), \underbrace{\left( \begin{array}{cccc} f_{X,\lambda}(s_1,s_1) &    f_{X,\lambda}(s_1,s_2) &  \ldots &  f_{X,\lambda}(s_1,s_k) \\
 f_{X,\lambda}(s_2,s_1) &    f_{X,\lambda}(s_2,s_2) & \ldots &  f_{X,\lambda}(s_2,s_k) \\ 
  \vdots & \vdots & \ldots & \vdots \\ f_{X,\lambda}(s_k,s_1) &    f_{X,\lambda}(s_k,s_2)& \ldots &  f_{X,\lambda}(s_k,s_k)\end{array} \right)}_{\displaystyle =\Sigma_\lambda} \Big),$$ where $\mathcal N_C$ denotes a circularly-symmetric complex Gaussian distribution with mean zero and complex-valued covariance matrix $ \Sigma_\lambda$. Furthermore, for two different frequencies $0 < \lambda_j \neq \lambda_k <\pi$, the corresponding vectors of finite Fourier transforms $J_{X,\lambda_j}$ and $J_{X,\lambda_k} $ are asymptotically independent; see, e.g., Theorem 5 in Cerovecki and Hörmann (2017). These properties of the finite Fourier transform as well as the fact that $ \widehat{p}_{X,\lambda} (\sigma,\tau) = J_{X,\lambda}(\sigma)\overline{J}_{X,\lambda}(\tau)$, for $\sigma, \tau \in [0,1]$, is the periodogram kernel, motivate the following bootstrap procedure to approximate the distribution of the test statistic $\mathcal U_T$ defined in (\[eq:test-stat-main\]).

1.   For $ \lambda_t=2\pi t/n$, $t=1,2, \ldots, N$, estimate the pooled spectral density operator $ \widehat{\mathcal F}_{\lambda_t}$ by $$\label{eq.sp-pooled}
      \widehat{\mathcal F}_{\lambda_t} = \frac{1}{2}\widehat{\mathcal F}_{X,\lambda_t} +  \frac{1}{2}\widehat{\mathcal F}_{Y,\lambda_t}$$ and denote by $ \widehat{f}_{\lambda_t}(\sigma,\tau)$, for $ \sigma, \tau \in \{s_1, s_2, \ldots, s_k\}$, the corresponding estimated pooled spectral density kernel.

2.   Generate two independent vectors $J^\ast_{X,\lambda_t} $ and $J^\ast_{Y,\lambda_t} $ as $$J^\ast_{X,\lambda_t} \sim \mathcal N_C(0, \widehat{\Sigma}_{\lambda_t})  \ \ \mbox{and} \ \  J^\ast_{Y,\lambda_t} \sim \mathcal N_C(0,  \widehat{\Sigma}_{\lambda_t}),$$ where $ \widehat{\Sigma}_\lambda $ is the matrix obtained by replacing in $ \Sigma_\lambda$ the unknown spectral density kernel $ f_{X,\lambda}$ by its pooled estimator $ \widehat{f}_\lambda$. For $\sigma, \tau \in \{s_1,s_2,\ldots, s_k\}$, let $$p^\ast_{X,\lambda_t}(\sigma,\tau) = J^\ast_{X,\lambda_t}(\sigma)\overline{J}^\ast_{X,\lambda_t}(\tau)  \ \ \mbox{and} \ \ p^\ast_{Y,\lambda_t}(\sigma,\tau) = J^\ast_{Y,\lambda_t}(\sigma)\overline{J}^\ast_{Y,\lambda_t}(\tau)$$ while, for $ t=-1,-2, \ldots, -N$, set $$p^\ast_{X,\lambda_t}(\sigma,\tau)=\overline{p}^\ast_{X,-\lambda_t}(\sigma, \tau) \;\; \mbox{and} \;\;  p^\ast_{Y,\lambda_t}(\sigma,\tau)=\overline{p}^\ast_{Y,-\lambda_t}(\sigma,\tau).$$

3.    For $\sigma, \tau \in \{s_1,s_2,\ldots, s_k\}$, let $$\hat f^\ast_{X,\lambda}(\sigma,\tau)=\frac{1}{b T}\sum_{t=-N}^{N} W\left(\frac{\lambda-\lambda_t}{b}\right)\, \widehat p^\ast_{X, \lambda_t} (\sigma, \tau)$$ and $$\hat f^\ast_{Y,\lambda}(\sigma,\tau)=\frac{1}{b T}\sum_{t=-N}^{N} W\left(\frac{\lambda-\lambda_t}{b}\right)\, \widehat p^\ast_{Y, \lambda_t} (\sigma, \tau).$$

4.    Approximate the distribution of the test statistic $ \mathcal U_T $ defined in (\[eq:test-stat-main\]) by that of the bootstrap test statistic $ \mathcal U^\ast_{T,k} $, where $$\mathcal U^\ast_{T,k}=\frac{2\pi}{Tk^2}\sum_{\substack{l=-N \\l\neq0}}^{N}\sum_{i,j=1}^k \Big| \widehat{f}^\ast_{X,\lambda_l}(s_i,s_j)
    -\widehat{f}^\ast_{Y,\lambda_l}(s_i,s_j) \Big|^2.$$

[The number $k$ of points in the interval $[0,1]$ for which the $k$-dimensional complex-valued random vectors $ J^\ast_{X,\lambda_t}$ and $ J^\ast_{Y,\lambda_t}$ are generated can be set equal to the sampling points $s_j$, $j=1,2, \ldots, k$, at which the observed functional random elements $ X_t(s_j)$ and $ Y_t(s_j)$ are observed in reality. However, and as it is commonly done in functional data analysis, these finite-dimensional vectors can be transformed to functional objects using a basis in $L^2$, for instance, the Fourier basis. In this case, the bootstrap approximation of the test statistic $ \mathcal U_T$ defined in (\[eq:test-stat-main\]) will then be given by $$\mathcal U^\ast_T=\frac{2\pi}{T}\sum_{\substack{l=-N\\l\neq0}}^{N}\int_0^1\int_0^1 \Big| \widehat{f}^\ast_{X,\lambda_l}(\tau,\sigma)
-\widehat{f}^\ast_{Y,\lambda_l}(\tau,\sigma) \Big|^2d\tau d\sigma. 
\label{eq:boot-test-stat}$$ From an asymptotic point of view, and as an inspection of the proof of Theorem \[t.boo-test\] below shows, both bootstrap approximations, $ \mathcal U^\ast_{T,k}$ and $ \mathcal U^\ast_{T} $, will lead to the same result, provided that for $ \mathcal U^\ast_{T,k}$ the number of points $ k$ increases to infinity as the sample size $T$ increases to infinity.]{}

Following the bootstrap procedure described in Steps 1-4, a bootstrap-based test then rejects ${\mathcal H_0}$ if $$t_{\mathcal U} \geq t^\ast_{\mathcal U,1-\alpha},$$ where $ t^\ast_{\mathcal U, 1-\alpha}$ denotes the upper $1-\alpha$ percentage point of the distribution of the bootstrap studentized test $$\label{tst-fanis}
t^\ast_{\mathcal U}=(\sqrt{b} T\,\mathcal U^\ast_T  -b^{-1/2}\widehat{\mu}^\ast_0)/\widehat{\theta}^\ast_0,$$ where $\mathcal U^\ast_T$ is defined in (\[eq:boot-test-stat\]) and $\widehat{\mu}_0^\ast$ and $ \widehat{\theta}^\ast_0$ are obtained by replacing the unknown spectral density kernel $ f_{X,\lambda}$ in the expressions for $\mu_0$ and $\theta_0$ given in Theorem \[t.asy-test\] by its pooled estimator $ \widehat{f}^\ast_{\lambda}(\sigma, \tau) = \widehat{f}^\ast_{X,\lambda}(\sigma, \tau)/2+ \widehat{f}^\ast_{Y,\lambda}(\sigma, \tau)/2 $, for all $\sigma, \tau \in [0,1]$. Notice that this distribution can be evaluated by Monte Carlo.

[It is worth mentioning that, by the definition of $ \widehat{\mu}^\ast_0$ and $ \widehat{\theta}^\ast_0$, the bootstrap studentized test $t^\ast_{\mathcal U}$ imitates correctly also the randomness in $t_{\mathcal U}$ which is introduced by replacing the unknown spectral density kernel $ f_{X,\lambda}$ appearing in $ \mu_0$ and $\theta_0$ by its pooled estimator $ \widehat{f}_\lambda$; see (\[eq.stud-test\]). A computationally simpler alternative will be to ignore this asymptotically negligible effect, that is, to use, instead of $t^\ast_{\mathcal U}$ given in (\[tst-fanis\]), the studentized version $t^+_{\mathcal U}=(\sqrt{b} T\,\mathcal U^\ast_T  -b^{-1/2}\widehat{\mu}_0)/\widehat{\theta}_0$ of the bootstrap-based test.]{}

Before describing the asymptotic behavior of the bootstrap test statistic $ \mathcal U^\ast_T$ defined in (\[eq:boot-test-stat\]), we state the following assumption which clarifies our requirements on the pooled spectral density kernel estimator $ \widehat{f}_\lambda$ used.

 : The pooled spectral density kernel estimator $ \widehat{f}_\lambda $ satisfies $$\sup_{\lambda_t\in [0,\pi]} \Big|\int_0^1 \int_0^1 \big(\widehat{f}_{\lambda_t}(\sigma,\tau) - f_{\lambda_t}(\sigma,\tau)\big)d\sigma d\tau \Big|  = o_P(\sqrt{b}), \quad \mbox{as} \;\; T \rightarrow \infty,$$ where $ f_\lambda$ is the spectral density kernel of the pooled spectral density operator $\mathcal F_\lambda=(1/2) \mathcal F_{X,\lambda} + (1/2) \mathcal F_{Y,\lambda} $.

Notice that the above assumption can be easily verified by using results for uniform consistency of spectral density estimators of univariate time series, since $$\begin{aligned}
\int_0^1\int_0^1 \widehat{f}_{X,\lambda} (\sigma,\tau)d\sigma d\tau & = \frac{1}{Tb}\sum_{t=-N}^N W\Big(\frac{\lambda-\lambda_t}{b}\Big) \int_0^1 \int_0^1 \widehat{p}_{X,\lambda_t}(\sigma,\tau) d\sigma d\tau\end{aligned}$$ can be interpreted as a kernel estimator of the spectral density of the univariate time series $ \int_0^1 X_t(s)ds$, $ t=1,2, \ldots, n$, the periodogram of which at frequency $ \lambda_t$ equals $\int_0^1 \int_0^1 \widehat{p}_{X,\lambda_t}(\sigma,\tau) d\sigma d\tau$. For instance, for the linear functional process $\{X_t, t\in\Z\}$ considered in this paper, $\int_0^1 X_t(s)ds$ is a univariate linear process as well and, under certain conditions, Assumption 2 is satisfied; see Franke and Härdle (1992). However, the same assumption can be fulfilled under different conditions on the integrated process $\int_0^1 X_t(s)ds$; see Wu and Zaffaroni (2015) for a discussion.\
The following theorem establishes the asymptotic validity of the suggested bootstrap procedure.\

\[t.boo-test\] Suppose that the Assumptions 1 and 2 as well as the conditions (i) and (ii) of Theorem \[t.asy-test\] are satisfied. Then, conditional on $X_1,\ldots,X_T,Y_1,\ldots, Y_T$, as $T\rightarrow \infty$, $$\sqrt{b} T \, \mathcal U^\ast_T-b^{-1/2} \widetilde{\mu}_0  \stackrel{d}{\rightarrow}  \mathcal N(0,\widetilde{\theta}_0),$$ in probability, where $$\widetilde{\mu}_0 = \frac{1}{\pi}\, \int_{-\pi}^\pi\left\{\int_0^1 f_{\lambda}(\tau,\tau) \,d\tau\,\right\}^2 d\lambda\int_{-\pi}^\pi W^2\left(u\right)\, du,$$ $$\begin{aligned}
\widetilde{\theta}_0^2&=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \left\{\int_0^1\int_0^1f_{\lambda}^2(\sigma, \tau)\,\,d\sigma\, d\tau\right\}^2\, d\lambda.\end{aligned}$$

Notice that, under $\mathcal H_0$, $ \mu_0=\widetilde{\mu}_0$ and $\theta^2_0=\widetilde{\theta}_0^2$ since $ \mathcal F_{X, \lambda}= \mathcal F_{Y, \lambda} $ (or, respectively, $ f_{X,\lambda}=f_{Y, \lambda}$). Thus, in this case, the asymptotic behavior of the test statistics $\mathcal U_T$ and $\mathcal U^\ast_T$ is identical, that is, the bootstrap procedure estimates consistently the distribution of the test statistic $ \mathcal U_T$ under $\mathcal H_0$. Furthermore, under $\mathcal  H_1$, the following holds true.

\[le.power\] Suppose that the conditions of Theorem \[t.boo-test\] are satisfied. Then, under $\mathcal H_1$, as $ T \rightarrow \infty$, $$\begin{aligned}
t_{\mathcal U} & =  
 \sqrt{b}T \int_{-\pi}^{\pi}\|{\mathcal F}_{X,\lambda} -{\mathcal F}_{Y,\lambda}\|^2_{HS} d\lambda + o_P(\sqrt{b} T) \ \rightarrow +\infty, \quad \mbox{in probability}.\end{aligned}$$

The above result, together with Theorem \[t.boo-test\] and Slusky’s theorem, imply that the power of the studentized test $t_{\mathcal U}$ based on the bootstrap critical values obtained from the distribution of the bootstrap studentized test $t^\ast_{\mathcal U}$ converges to unity as $ T\rightarrow \infty$, i.e., the studentized test $t_{\mathcal U}$ is consistent.

Numerical Results
=================

Choice of the Smoothing Parameter
---------------------------------

Implementing the studentized test $t_{\mathcal U}$ requires the choice of the smoothing bandwidth $b$. For univariate and multivariate time series, this issue has been investigated in the context of a cross-validation type criterion by Beltrão and Bloomfield (1987), Hurvich (1985) and Robinson (1991). However, adaption of the multivariate approach of Robinson (1991) to the spectral density estimator $\widehat{f}_{X,\lambda}(\sigma_r,\tau_s)$, for $r,s \in \{1,2, \ldots, k\}$, faces problems due to the high dimensionality of the periodogram operator involved.

We propose a simple approach to select the bandwidth $b$ used in our testing procedure which is based on the idea to overcome the high-dimensionality of the problem by selecting a single bandwidth based on the “on average" behavior of the pooled estimator $ 
  \widehat{f}_{\lambda}(\sigma_r,\tau_s)$, that is, its behavior over all points $r,s \in \{1,2, \ldots, k\}$ for which the functional objects are observed. To elaborate, define first the following quantities. The averaged periodogram $$\widehat{I}_T(\lambda) = \frac{1}{k^2}\sum_{r=1}^k\sum_{s=1}^k\big\{
  \frac{1}{2}\widehat{p}_{X,\lambda}(\sigma_r,\tau_s) + \frac{1}{2}\widehat{p}_{Y,\lambda}(\sigma_r,\tau_s)\big\}$$ and the averaged spectral density estimator $$\widehat{g}_b(\lambda) =\frac{1}{k^2}\sum_{r=1}^k\sum_{s=1}^k\big\{
  \frac{1}{2}\widehat{f}_{X,\lambda}(\sigma_r,\tau_s) + \frac{1}{2}\widehat{f}_{Y,\lambda}(\sigma_r,\tau_s)\big\}.$$ Notice that $\widehat{I}_T(\lambda)$ can be interpreted as the periodogram at frequency $ \lambda$ of the pooled univariate process $ \{V_t= \frac{1}{2}\int_0^1 X_t(s)ds + \frac{1}{2}\int_0^1 Y_t(s)ds, t \in {\mathbb Z}\}$ while $\widehat{g}_b(\lambda)$ is an estimator of the spectral density $g$ of $ \{V_t,t\in\Z\}$. We then choose the bandwidth $b$ by minimizing the objective function $$CV(b) = \frac{1}{N} \sum_{t=1}^N \big\{ \log(\widehat{g}_{-t}(\lambda_t) + \widehat{I}_T(\lambda_t)/\widehat{g}_{-t}(\lambda_t)\big\},$$ over a grid of values of $b$, where $\widehat{g}_{-t}(\lambda_t) = (Tb)^{-1}\sum_{s\in N_t} W((\lambda_t-\lambda_s)/b) \widehat{I}_T(\lambda_s)$ and $N_t=\{s: -N \leq s \leq N \ \mbox{and} \ s \neq \pm t\}$. That is, $ \widehat{g}_{-t}(\lambda_t)$ is the leave-one-out kernel estimator of $ g(\lambda)$, i.e., the estimator obtained after deleting the $t$-th frequency; see also Robinson (1991).\
Due to the computational complexity of the simulation analysis studied in the next section, the use of this automatic choice of the bandwidth $b$ will only be illustrated in the real-life data example.

Monte-Carlo Simulations
-----------------------

We generated functional time series stemming from the following functional linear processes $$X_t=\varepsilon_t + a_1 \varepsilon_{t-1}+ a_2 \varepsilon_{t-2}, 
\label{mod-1}$$ and $$\quad Y_t=\varepsilon_t + a_1 \varepsilon_{t-1},
\label{mod-2}$$ $ t=1,2,\ldots,T$, where $\varepsilon_t$’s are i.i.d. Brownian Bridges. All curves were approximated using 21 equidistant points in the unit interval and transformed into functional objects using the Fourier basis with 21 basis functions. Three sample sizes $T=50$, $T=100$ and $T=200$ were considered and the bootstrap test was applied using three nominal levels, $\alpha=0.01$, $\alpha=0.05$ and $\alpha=0.10$. All bootstrap calculations were based on $B=1,000$ bootstrap replicates and $R=500$ model repetitions. To investigate the empirical size and power behavior of the bootstrap test, we set $a_1=0.8$ and consider a selection of $a_2$ values, i.e., $a_2 \in \{0.0, 0.2, 0.5, 0.8, 1.0\}$, and various bandwidths $b$. (Notice that $a_2=0$ corresponds to the null hypothesis.)

We first demonstrate the ability of the bootstrap procedure to approximate the distribution of the test statistic under the null. For this, and in order to estimate the exact distribution of the studentized test $t_\mathcal U$ (see (\[eq.stud-test\])), 10,000 replications of the process (\[mod-1\]) and (\[mod-2\]) with $a_2=0$ have been generated, and a kernel density estimate of this exact distribution has been obtained using a Gaussian kernel with bandwidth $h$. The suggested bootstrap procedure is then applied to three randomly selected time series and the bootstrap studentized test $ t^\ast_\mathcal U$ (see (\[tst-fanis\])) has been calculated. Two sample sizes of $ T=50$ and $T=500$ observations have been considered. Figure 5.1 shows the results obtained together with the approximation of the distribution of $t_\mathcal U$ provided by the central limit theorem, i.e., the $\mathcal N(0,1)$ distribution. As it can be seen from this figure, the convergence towards the asymptotic Gaussian distribution is very slow. Even for sample sizes as large as $T=500$, the exact distribution retains its skewness which is not reproduced by the $\mathcal N(0,1)$ distribution. In contrast to this, the bootstrap approximations are very good and the behavior of the exact densities, especially in the critical right hand tale of this distribution, are very accurate.

\[fig.dens\]

![Density plots of the estimated exact standardized distribution of $t_{\mathcal U}$ (red line), the standard Gaussian distibution (black line) and three bootstrap approximations (blue lines). Left panel, $T$=50 ($h=0.2$), right panel $T$=500 ($h=0.04$).](n50h02N.pdf "fig:"){height="6.0cm" width="7.0cm"} ![Density plots of the estimated exact standardized distribution of $t_{\mathcal U}$ (red line), the standard Gaussian distibution (black line) and three bootstrap approximations (blue lines). Left panel, $T$=50 ($h=0.2$), right panel $T$=500 ($h=0.04$).](n500h004.pdf "fig:"){height="6.0cm" width="7.0cm"}

We next investigate the finite sample size and power behavior of the bootstrap studentized test under the aforementioned variety of process parameters and three different sample sizes, $T=50$, $T=100$ and $T=200$. The results obtained are shown in Table \[tab:empsize-power\]. As it is evident from this table, the bootstrap studentized test has a very good empirical size and power behaviour even in the case of $T=50$ observations. The empirical sizes are close to the nominal ones and the empirical power of the test increases to one as the deviations from the null become larger (i.e., larger values of $a_2$) and/or the sample size increases.

  ----- ------- --------------- --------------- --------------- -- --------------- --------------- ---------------
                                     b=0.2                                              b=0.3      
   $T$   $a_2$   $\alpha=0.01$   $\alpha=0.05$   $\alpha=0.10$      $\alpha=0.01$   $\alpha=0.05$   $\alpha=0.10$
   50     0.0        0.014           0.050           0.082              0.012           0.052           0.102
          0.2        0.014           0.062           0.124              0.008           0.064           0.130
          0.5        0.044           0.144           0.234              0.046           0.124           0.226
          0.8        0.074           0.362           0.468              0.098           0.262           0.372
          1.0        0.150           0.372           0.512              0.154           0.336           0.464
                                     b=0.1                                              b=0.2      
         $a_2$   $\alpha=0.01$   $\alpha=0.05$   $\alpha=0.10$      $\alpha=0.01$   $\alpha=0.05$   $\alpha=0.10$
   100    0.0        0.004           0.046           0.100              0.014           0.050           0.084
          0.2        0.012           0.074           0.126              0.022           0.076           0.130
          0.5        0.056           0.190           0.296              0.076           0.222           0.330
          0.8        0.164           0.402           0.554              0.210           0.412           0.574
          1.0        0.330           0.572           0.688              0.366           0.624           0.730
                                    b=0.06                                              b=0.1      
         $a_2$   $\alpha=0.01$   $\alpha=0.05$   $\alpha=0.10$      $\alpha=0.01$   $\alpha=0.05$   $\alpha=0.10$
   200    0.0        0.008           0.055           0.106              0.014           0.048           0.098
          0.2        0.016           0.088           0.144              0.024           0.100           0.186
          0.5        0.126           0.280           0.428              0.174           0.392           0.502
          0.8        0.394           0.646           0.752              0.448           0.710           0.810
          1.0        0.636           0.864           0.926              0.692           0.888           0.952
  ----- ------- --------------- --------------- --------------- -- --------------- --------------- ---------------

  : Empirical size and power of the bootstrap studentized test for the equality of two spectral density operators. The functional time series in each group were generated according to models (\[mod-1\]) and (\[mod-2\]).[]{data-label="tab:empsize-power"}

A Real-Life Data Example
------------------------

We applied the bootstrap studentized test to a data set consisting of temperature measurements recorded in Nicosia, Cyprus, for the winter period, December 2006 to beginning of March 2007 and for the summer period, June 2007 to end of August 2007. It is well-known that the mean temperatures during winter periods are smaller than those of summer periods. Our aim is to test whether there is also a significant difference in the autocovariance structure of the winter and summer periods. The data consists of two samples of curves $\{(X_{t}, Y_t),\;t=1,2,\ldots, 92\}$, where $X_{t}$ represents the temperature of day $t$ for Dec2006-Jan2007-Feb2007-March2007 and $Y_t$ for Jun2007-Jul2007-Aug2007. More precisely, $X_{1}$ represents the temperature of the 1st of December 2006 and $X_{92}$ the temperature of the 2nd of March 2007, whereas $Y_{1}$ represents the temperature of the 1st of June 2007 and $Y_{92}$ the temperature of the 31st of August 2007. The temperature recordings were taken in $15$ minutes intervals, i.e., there are $k=96$ temperature measurements for each day for a total of $T=92$ days in both groups. These measurements were transformed into functional objects using the Fourier basis with 21 basis functions. All curves were rescaled in order to be defined in the unit interval. Figure \[fig:win-sum\] shows the centered temperatures curves of the winter and summer periods, i.e., the curves in each group are transformed by subtracting the corresponding group sample mean functions.

Using the cross-validation algorithm described in Section 5.1, the choice of bandwidth is $b_{CV}=0.075$ (see Figure \[fig.fre-plot\](a)) and the resulting bootstrap $p$-value is equal to 0.030, leading to the rejection of the null hypothesis for almost all commonly used $\alpha$-levels. This implies that dependence properties as measured by autocovariances of the temperature measurements of the winter period differ significantly from those of the summer period.

\[fig.fre-plot\]

In order to understand the reasons leading to this rejection, we decompose the standardized test $ t_{\mathcal U}$ after ignoring the centering sequence $b^{-1/2} \mu_0$ and approximating the integral of the (squared) Hilbert-Schmidt norm by the corresponding Riemann sum over the Fourier frequencies $ \lambda_j=2\pi j/T$, as follows: $$\begin{aligned}
 \label{eq.dec}
\sqrt{b}T \, \mathcal U_T\Big/\theta_0 & \approx  2\pi \sqrt{b}\sum_{j=-N}^{N}  \|\widehat{\mathcal F}_{X,\lambda_j} - \widehat{\mathcal F}_{Y,\lambda_j}\|^2_{HS}\Big/\theta_0 \nonumber \\
& = \sum_{j=-[(T-1)/2]}^{[(T-1)/2]} Q_{T,\lambda_j},\end{aligned}$$ where $$Q_{T,\lambda_j} = 2\pi \sqrt{b} \, \|\widehat{\mathcal F}_{X,\lambda_j} - \widehat{\mathcal F}_{Y,\lambda_j}\|^2_{HS}\Big/\theta_0  \geq 0.$$ Expression (\[eq.dec\]) shows the contributions of the differences $  \|\widehat{\mathcal F}_{X,\lambda_j} - \widehat{\mathcal F}_{Y,\lambda_j}\|^2_{HS} $ for each frequency $ \lambda_j$ to the total value of the test statistic $ \mathcal U_T$. Large values of $ Q_{T,\lambda_j} $ pinpoint to frequency regions from which large contributions to the test statistic $ \mathcal U_T$ occur. Notice that $Q_{T,\lambda_j} $ can be estimated by $ \widehat{Q}_{T,\lambda_j} = 2\pi \sqrt{b}\,  \|\widehat{\mathcal F}_{X,\lambda_j} - 
\widehat{\mathcal F}_{Y,\lambda_j}\|^2_{HS}\Big/\widehat{\theta}_0 $, where $ \widehat{\theta}_0$ is the estimator of $ \theta_0$ discussed in Section 3. A plot of the estimated quantities $\widehat{Q}_{T,\lambda_j} $ against the frequencies $ \lambda_j$, $j=0,1, \ldots, N$, is, therefore, very informative in identifying frequency regions where the differences between the two spectral density operators are large and for interpreting the results of the testing procedure.

Figure \[fig.fre-plot\](b) shows for the real-life temperature data example considered the plot of $ \widehat{Q}_{T, \lambda_j}$ at a log-scale. As it is easily seen from this figure, the large value of the test statistic $ \mathcal U_T$ which leads to a rejection of the null hypothesis is mainly due to large differences between the two spectral density operators at low frequencies. That is, differences in the long term periodicities between the winter and the summer temperature curves seem to be the main reason for rejecting the null hypothesis.

Auxiliary Results and Proofs
============================

[**Proof of Lemma \[l.ex-SDO\]**]{}. Under the assumptions of the lemma, straightforward calculations yield $$\sum_{t\in \Z} \|r_{X,t}\|_2\leq \sum_{t,j\in \Z} a_j a_{j+t} \|E\varepsilon_0^2\|_2^2<\infty.$$ An analogous result holds for $r_{Y,t}$. Hence, the assertions of the lemma is an immediate consequence of Proposition 2.1 in Panaretos and Tavakoli (2013).

The proof of the Theorem \[t.asy-test\] is based on several auxiliary lemmas:

\[l.centering\] Under the assumptions of Theorem \[t.asy-test\], it holds under $\mathcal H_0$, $$\sqrt b T\int_{-\pi}^\pi\|E\widehat{\mathcal F}_{X,\lambda}- {\mathcal F}_{X,\lambda}\|_{HS}^2\, d\lambda+\sqrt b T\int_{-\pi}^\pi\|E\widehat{\mathcal F}_{Y,\lambda}- {\mathcal F}_{Y,\lambda}\|_{HS}^2\, d\lambda=o(1), \quad \mbox{as} \;\; T \rightarrow \infty.$$

For symmetry reasons, we consider $\sqrt b T\int_{-\pi}^\pi\|E\widehat{\mathcal F}_{X,\lambda}- {\mathcal F}_{X,\lambda}\|_{HS}^2$ only. Rewriting the Hilbert Schmidt norm of an integral operator as the double integral over the squared kernel, we obtain $$\begin{aligned}
&\sqrt b T\int_{-\pi}^\pi\|E\widehat{\mathcal F}_{X,\lambda}- {\mathcal F}_{X,\lambda}\|_{HS}^2\, d\lambda\\
&\leq \sqrt b T\int_0^1\int_0^1 |E\varepsilon_0(\sigma)\varepsilon_0(\tau)|^2d\sigma\, d\tau\\
&\quad \times\int_{-\pi}^\pi\left|\frac{1}{2\pi}\sum_{j,h\in\Z}a_j\, a_{j+h}\left[e^{-ih\lambda}-\frac{1}{bT}\sum_{t=-N}^{N}W\left(\frac{\lambda-\lambda_t}{b}\right)e^{-ih\lambda_t}\left(1-\frac{|h|}{T}\right)\I_{|h|<T}\right]\right|^2\, d\lambda.
\end{aligned}$$ Introducing $\widetilde X_t=\sum_{k=1}^\infty a_k \eta_{t-k}$ with real-valued i.i.d. white noise $(\eta_t)_t$ satisfying $E\eta_t^2=1$ and $E\eta_t^4<\infty$, we obtain, for some $C\in (0,\infty)$, $$\sqrt b T\int_{-\pi}^\pi\|E\widehat{\mathcal F}_{X,\lambda}- {\mathcal F}_{X,\lambda}\|_{HS}^2\, d\lambda
= C \, \sqrt b T \int_{-\pi}^\pi\left|E\hat f^{\widetilde X}(\lambda)-f^{\widetilde X}(\lambda) \right|^2\, d\lambda,$$ where $f^{\widetilde X}$ denotes the spectral density of the univariate time series $(\widetilde X_t)_t$ and $$\hat f^{\widetilde X}(\lambda)=\frac{1}{bT}\sum_{t=-N}^{N}W\left(\frac{\lambda-\lambda_t}{b}\right)
\frac{1}{2\pi T}\sum_{t_1,t_2=1}^T \widetilde X_{t_1}\widetilde X_{t_2}\, \exp(i\lambda_t( t_1- t_2)).$$ is the corresponding smoothed periodogram. Using classical results for the smoothed periodogram of univariate linear processes (e.g. following the lines of the proof on page 396 in Taniguchi and Kakizawa (2000)), we obtain that $$\sqrt b T \int_{-\pi}^\pi\left|E\hat f^{\widetilde X}(\lambda)-f^{\widetilde X}(\lambda) \right|^2\, d\lambda=o(1)$$ under our assumptions on the bandwidths.

\[l.cumulant\] Under the assumptions of Theorem \[t.asy-test\], it holds for $\widehat p_{X, \lambda}^c=\widehat p_{X, \lambda}-E\widehat p_{X, \lambda}$ and any $t_1,\dots, t_4\in \{-N,\dots, N\}\backslash\{0\}$ $$\begin{aligned}
&cum(p_{X, \lambda_{t_1}}^c(\sigma_1,\tau_1),p_{X, \lambda_{t_2}}^c(\sigma_2,\tau_2), p_{X, \lambda_{t_3}}^c(\sigma_3,\tau_3), p_{X, \lambda_{t_4}}^c(\sigma_4,\tau_4) )=O(T^{-M/2})\,\prod_{i=1}^4 f_i(\sigma_i)\, g_i(\tau_i)\;
\end{aligned}$$ for $M=\#\{|t_1|,\,|t_2|,\, |t_3|,\,|t_4|\}$ and some $f_i,\; g_i\in L_\R^8([0,1],\mu),\; i=1,\dots 4.$

We have with $(\varepsilon_j\varepsilon_k)^c(\sigma,\tau)=\varepsilon_j(\sigma)\varepsilon_k(\tau)-E(\varepsilon_j(\sigma)\varepsilon_k(\tau))$, $$\label{eq.cum}
\begin{aligned}
&cum(p_{X, \lambda_{t_1}}^c(\sigma_1,\tau_1),p_{X, \lambda_{t_2}}^c(\sigma_2,\tau_2), p_{X, \lambda_{t_3}}^c(\sigma_3,\tau_3), p_{X, \lambda_{t_4}}^c(\sigma_4,\tau_4) )\\
&=\frac{1}{16\pi^4T^4}\sum_{s_1,\dots, s_8=1}^T\sum_{k_1,\dots, k_8\in \Z} a_{k_1}\cdots a_{k_8}\, e^{2\pi i (s_1-s_2)t_1/T}
e^{2\pi i (s_4-s_3)t_2/T}e^{2\pi i (s_6-s_5)t_3/T}e^{2\pi i (s_7-s_8)t_4/T}\\
&\times cum((\varepsilon_{s_1-k_1}\varepsilon_{s_2-k_2})^c(\sigma_1,\tau_1),\,
(\varepsilon_{s_3-k_3}\varepsilon_{s_4-k_4})^c(\sigma_2,\tau_2),\,
(\varepsilon_{s_5-k_5}\varepsilon_{s_6-k_6})^c(\sigma_3,\tau_3),\,
(\varepsilon_{s_7-k_7}\varepsilon_{s_8-k_8})^c(\sigma_4,\tau_4)).
\end{aligned}$$ We now consider terms with different numbers of equal indices of the innovations separately. Obviously, all summands, where either all the indices are equal or there only two different indices, sum up to $O(T^{-2})$. If there are three different indices, then two situations have a non-zero contribution.

1.  2 couples and 1 quad: We consider the case $$s_1-k_1=s_3-k_3=s_5-k_5=s_7-k_7,~ s_2-k_2=s_4-k_4, ~ s_6-k_6= s_8-k_8.$$ The corresponding part on the r.h.s. of can be bounded by $$\begin{aligned}
    &\left(\prod_{i=1}^4 f_i(\sigma_i)\, g_i(\tau_i)\right)\, \frac{1}{T^4} \sum_{s_1,\dots, s_8=1}^Te^{2\pi i(s_1t_1-s_3t_2-s_5 t_3+s_7 t_4)/T}\, e^{2\pi i (s_4t_2-s_2 t_1)/T}\, e^{2\pi i (s_6t_3-s_8 t_4)/T}\\
    &\quad\times\sum_{k_1, k_2, k_6\in \Z}a_{k_1}a_{k_1+s_3-s_1}a_{k_1+s_5-s_1} a_{k_1+s_7-s_1} a_{k_2} a_{k_2+s_4-s_2} a_{k_6}a_{k_6+s_8-s_6} 
    \\
    &\leq\, O(T^{-M/2})\,\prod_{i=1}^4 f_i(\sigma_i)\, g_i(\tau_i) 
    \end{aligned}$$ with the same arguments as used by Brockwell and Davis (1991) to bound their formula (10.3.18). All other situations can either be treated similarly or have zero contribution.

2.  2 triples and 1 couple: We consider one scenario in detail, all others can be handled analogously. Let $$s_1-k_1=s_3-k_3=s_5-k_5,~ s_4-k_4=s_6-k_6=s_8-k_8,~  s_2-k_2=s_7-k_7.$$ The corresponding part on the r.h.s. of can be bounded by $$\begin{aligned}
    &\,\left(\prod_{i=1}^4 f_i(\sigma_i)\, g_i(\tau_i) \right)\, \frac{1}{T^4} \sum_{s_1, \dots, s_8=1}^T
    e^{2\pi i(s_1t_1-s_3t_2-s_5t_3)/T}\, e^{2\pi i (s_4t_2-s_6 t_3-s_8 t_4)/T}\, e^{2\pi i (s_7 t_4-s_2t_1)/T}
    \\&\quad\times \sum_{k_1, k_2, k_4\in \Z} a_{k_1}a_{k_1+s_3-s_1}a_{k_1+s_5-s_1}  a_{k_4} a_{k_4+s_6-s_4}a_{k_4+s_8-s_4}a_{k_2} a_{k_2+s_7-s_2}\\
    &\leq\,  O(T^{-M/2}) \,\prod_{i=1}^4 f_i(\sigma_i)\, g_i(\tau_i) .
    \end{aligned}$$

If there are four different indices, then only the case of four couples has a nonzero distribution. This situation can be handled analogously as the case of 1 quad and 2 couples.

Finally, the cumulant on the r.h.s. of vanishes in the case of more than four different indices.

\[l.bd-thm10.3.1\] Under the assumptions of Theorem \[t.asy-test\], we have for any Fourier frequency $\lambda_t\in[0,\pi],$ $$\widehat p_{X,\lambda_t}(\sigma,\tau)\,=\, \left|\sum_{j=-\infty}^\infty a_j\, e^{-i\lambda_t j}\right|^2 \widehat p_{\varepsilon,\lambda_t}(\sigma,\tau) \,+\, R_{T,\lambda_t}^X(\sigma,\tau),$$ where $\widehat p_{\varepsilon,\lambda}$ denotes the periodogram kernel of the innovation process $(\varepsilon_t)_t$ and\
$ \int_0^1\int_0^1 \sup_{\lambda_t\in[0,\pi]}\,E |R_{T,\lambda_t}^X(\sigma,\tau)|^4\,d\sigma\, d\tau=O(T^{-2})$.

As in Theorem 10.3.1 in Brockwell and Davis (1991), we have $$\frac{1}{\sqrt{2\pi T}}\sum_{s=1}^T X_s(\sigma)\, e^{-is\lambda_t}= \sum_{j=-\infty}^\infty a_j\, e^{-i\lambda_t j}\, \frac{1}{\sqrt{2\pi T}}\sum_{s=1}^T \varepsilon_s(\sigma)\, e^{-is\lambda_t}
+D_{T,\lambda}(\sigma)$$ with $$D_{T,\lambda}(\sigma)=\frac{1}{\sqrt{2\pi T}} \sum_{j=-\infty}^\infty a_j\, e^{-i\lambda_t j}\left\{\sum_{s=1-j}^{T-j} \varepsilon_s(\sigma)\, e^{-is\lambda_t}- \sum_{s=1}^{T} \varepsilon_s(\sigma)\, e^{-is\lambda_t}\right\}.$$ This gives $$\begin{aligned}
|R_{T,\lambda}^X(\sigma,\tau)|\leq &C\, \Bigg\{\left|\frac{1}{\sqrt{2\pi T}}\sum_{s=1}^T \varepsilon_s(\sigma)\, e^{-is\lambda_t}\, D_{T,-\lambda_t}(\tau) \right|
+\left|\frac{1}{\sqrt{2\pi T}}\sum_{s=1}^T \varepsilon_s(\tau)\, e^{is\lambda_t}\, D_{T,\lambda_t}(\sigma)\right| \\
&+ |D_{T,\lambda_t}(\sigma) \, D_{T,-\lambda_t}(\tau)|\Bigg\}
\end{aligned}$$ for some $C\in[0,\infty)$. Hence, it remains to show that $\int_0^1 \sup_{\lambda_t\in[0,\pi]}  E|D_{T,\lambda_t}(\sigma)|^8\, d\sigma=O(T^{-4}).$ This can be done along the lines of the proof of Lemma 1 in Paparoditis (2000) and is omitted here.

\[l.mean\] Suppose that the assumptions of Theorem \[t.asy-test\] hold true. Then $$\sqrt bT\, M_T-b^{-1/2} \, \mu_0=o_P(1),$$ where $$\label{eq.MT}
\begin{aligned}
M_T&=\frac{1}{b^2T^2}\int_{-\pi}^\pi\int_0^1\int_0^1\sum_{t=-N}^N \, W^2\left(\frac{\lambda-\lambda_t}{b}\right)\,
\left|\widehat p_{X, \lambda_t}^c(\sigma,\tau)\,-\, \widehat p_{Y,\lambda_t}^c(\sigma,\tau)\right|^2\, d\sigma\, d\tau\, d\lambda
\end{aligned}$$ and $p_{X,\lambda_t}^c(\sigma,\tau)= p_{X, \lambda_t}(\sigma,\tau)-Ep_{X, \lambda_t}(\sigma,\tau).$

We show that mean and variance of $\sqrt bT\, M_T-b^{-1/2} \, \mu_0$ tend to zero asymptotically. First, we have $$\begin{aligned}
&\frac{1}{b^{3/2} T}\int_{-\pi}^\pi\int_0^1\int_0^1\sum_{t=-N}^N \, W^2\left(\frac{\lambda-\lambda_t}{b}\right)\,
E\left|\widehat p_{X, \lambda_t}^c(\sigma,\tau)\right|^2\, d\sigma\, d\tau\, d\lambda\\
&=\frac{1}{b^{3/2} T}\int_{-\pi}^\pi\sum_{t=-N}^N W^2\left(\frac{\lambda-\lambda_t}{b}\right)\, d\lambda\\
&\quad \times\int_0^1\int_0^1\frac{1}{4\pi^2 T^2}\sum_{s_1,\dots, s_4=1}^T e^{2\pi i(s_1-s_2-s_3+s_4)t/T} \, cov(X_{s_1}(\sigma)X_{s_2}(\tau), X_{s_3}(\sigma) X_{s_4}(\tau))\, d\sigma\, d\tau\\
&=o(1) \,+\, \frac{1}{4\pi^2\,b^{3/2} T^3}\int_{-\pi}^\pi\sum_{t=-N}^N W^2\left(\frac{\lambda-\lambda_t}{b}\right)\, d\lambda\, \sum_{k_1,\dots, k_4\in \Z}a_{k_1}\, a_{k_2} \, a_{k_3}\, a_{k_4}\\
&\quad\times\sum_{\substack{s_1,\dots, s_4=1\\s_1-k_1=s_3-k_3\\s_2-k_2=s_4-k_4}}^T e^{2\pi i(s_1-s_2-s_3+s_4)t/T}\, \int_0^1\int_0^1cov(\varepsilon_{s_1-k_1}(\sigma)\varepsilon_{s_2-k_2}(\tau), \varepsilon_{s_3-k_3}(\sigma) \varepsilon_{s_4-k_4}(\tau))\, d\sigma\, d\tau\\
&\quad +\frac{1}{4\pi^2\,b^{3/2} T^3}\int_{-\pi}^\pi\, \sum_{t=-N}^N \, W^2\left(\frac{\lambda-\lambda_t}{b}\right)\, d\lambda\, \sum_{k_1,\dots, k_4\in \Z}a_{k_1}\, a_{k_2} \, a_{k_3}\, a_{k_4}\\
&\quad\times\sum_{\substack{s_1,\dots, s_4=1\\s_1-k_1=s_4-k_4\\s_2-k_2=s_3-k_3}}^T e^{2\pi i(s_1-s_2-s_3+s_4)t/T}\, \int_0^1\int_0^1cov(\varepsilon_{s_1-k_1}(\sigma)\varepsilon_{s_2-k_2}(\tau), \varepsilon_{s_3-k_3}(\sigma) \varepsilon_{s_4-k_4}(\tau))\, d\sigma\, d\tau\\
&=: o(1)+M_{T,1} +M_{T,2}
\end{aligned}$$ with an obvious notation for $M_{T,1}$ and $M_{T,2}$. We consider the two remaining parts separately. Using the definition of $(\widetilde X_t)_t$ in the proof of Lemma \[l.centering\] and the summability of $\gamma_{\widetilde X}$, we obtain $$\begin{aligned}
M_{T,1}&= \left(\int_0^1 E\varepsilon_1^2(\sigma) d\sigma\right)^2
\frac{1}{4\pi^2b^{3/2}T^3}\int_{-\pi}^\pi \sum_{t=-N}^N W^2\left(\frac{\lambda-\lambda_t}{b}\right)\, d\lambda\, 
\\
&\quad
\times\sum_{k_1\in \Z}a_{k_1}\sum_{\substack{s_1,s_3=1}}^T\, a_{k_1+s_3-s_1} e^{2\pi i(s_1-s_3)t/T}\, \sum_{k_2\in \Z}\, a_{k_2} \sum_{\substack{s_2, s_4=1}}^T\,  a_{k_2+s_4-s_2} e^{2\pi i(s_4-s_2)t/T}\, \\
&=\left(\int_0^1 E\varepsilon_1^2(\sigma) d\sigma\right)^2
\frac{1}{4\pi^2b^{1/2}T^3}\int_{-\pi/b}^{\pi/b}\sum_{t=-N}^N W^2\left(\lambda-\frac{\lambda_t}{b}\right)\, d\lambda\, 
\\
&\quad
\times \sum_{\substack{s_1,s_3=1}}^T\, \gamma_{\widetilde X}(s_3-s_1) e^{2\pi i(s_1-s_3)t/T}\, \sum_{\substack{s_2, s_4=1}}^T\,  \gamma_{\widetilde X}(s_4-s_2) e^{2\pi i(s_4-s_2)t/T}.
\end{aligned}$$ Hence, with similar arguments as in the proof of Lemma 3 in Paparoditis (2000), $$\begin{aligned}
M_{T,1}-\frac{1}{2b^{1/2}} \mu_0
&=\frac{1}{4\pi^2b^{1/2}}\,\left(\int_0^1 E\varepsilon_1^2(\sigma) d\sigma\right)^2
\Bigg\{\frac{1}{T^3}\int_{-\pi/b}^{\pi/b}\sum_{t=-N}^N W^2\left(\lambda-\frac{\lambda_t}{b}\right)\, d\lambda\, 
\\
&\quad
\times \sum_{\substack{s_1,s_3=1}}^T\, \gamma_{\widetilde X}(s_3-s_1) e^{2\pi i(s_1-s_3)t/T}\, \sum_{\substack{s_2, s_4=1}}^T\,  \gamma_{\widetilde X}(s_4-s_2) e^{2\pi i(s_4-s_2)t/T}
\\
&\quad -\,
  \int_{-\pi}^\pi W^2(u)\, du \,\frac{1}{2\pi}\,\int_{-\pi}^\pi\left(\sum_{h\in \Z}\gamma_{\widetilde X}(h)e^{ ih\lambda}\right)^2\, d\lambda\Bigg\}\\
&\leq o(1)\,+ \, \frac{1}{4\pi^2}\,\left(\int_0^1 E\varepsilon_1^2(\sigma) d\sigma\right)^2\int_{-\pi}^\pi W^2(u)\, du\\
&\quad
\frac{1}{b^{1/2}}\Bigg\{\frac{1}{T^3}\sum_{t=-N}^N
\left(\sum_{\substack{s_1,s_3=1}}^T\, \gamma_{\widetilde X}(s_3-s_1) e^{2\pi i(s_1-s_3)t/T}\right)^2
\\
&\quad -\,
   \,\frac{1}{2\pi 
    }\,\int_{-\pi}^\pi\left(\sum_{h\in \Z}\gamma_{\widetilde X}(h)e^{ ih\lambda}\right)^2\, d\lambda\Bigg\}\\
&\ninfty 0.
\end{aligned}$$ Following the lines of the proof of Theorem 10.3.2 in Brockwell and Davis (1991) (in particular using their methods to bound (10.3.18)), we obtain $$\begin{aligned}
M_{T,2}&= \int_0^1\int_0^1 \left(E(\varepsilon_1(\sigma)\varepsilon_1(\tau))\right)^2 d\sigma\,d\tau
\frac{1}{4\pi^2b^{3/2}T^3}\int_{-\pi}^\pi \sum_{t=-N}^N W^2\left(\frac{\lambda-\lambda_t}{b}\right)\, d\lambda\, 
\\
&\quad
\times\sum_{k_1\in \Z}a_{k_1}\sum_{\substack{s_1,s_4=1}}^T\, a_{k_1+s_4-s_1} e^{2\pi i(s_1+s_4)t/T}\, \sum_{k_2\in \Z}\, a_{k_2} \sum_{\substack{s_2, s_3=1}}^T\,  a_{k_2+s_3-s_2} e^{-2\pi i(s_3+ s_2)t/T}\\
&= O\left((bT)^{-1/2}\right).
\end{aligned}$$ In complete analogy to the above calculations, we get, from $f_{X,\lambda}=f_{Y,\lambda}$, $$\begin{aligned}
&\frac{1}{b^{3/2} T}\int_{-\pi}^\pi\int_0^1\int_0^1\sum_{t=-N}^N \, W^2\left(\frac{\lambda-\lambda_t}{b}\right)\,
E\left|\widehat p_{Y, \lambda_t}^c(\sigma,\tau)\right|^2\, d\sigma\, d\tau\, d\lambda-\frac{1}{2b^{1/2}} \mu_0=0.
\end{aligned}$$ Hence, from independence of the time series $(X_t)_t$ and $(Y_t)_t$ we get $
E M_T\ninfty 0.
$

We turn now to $var(M_T)$. First, note that by independence of the time series $(X_t)_t$ and $(Y_t)_t$, and in complete analogy to the proof of Theorem 10.3.2 in Brockwell and Davis (1991), that for $t_1\neq t_2$ $$\begin{aligned}
&cov\left(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\,{\widehat p_{Y, -\lambda_{t_1}}^c(\sigma,\tau)}\,,\,
\widehat p_{X, \lambda_{t_2}}^c(\sigma,\tau)\,{\widehat p_{Y, -\lambda_{t_2}}^c(\sigma,\tau)}\right)\\
&=cov(\widehat p_{X, \lambda_{t_1}}(\sigma,\tau)\,,\, \widehat p_{X, \lambda_{t_2}}(\sigma,\tau) )\,
cov({\widehat p_{Y,- \lambda_{t_1}}(\sigma,\tau)}\,,\, {\widehat p_{Y, -\lambda_{t_2}}(\sigma,\tau)})\\
&=f(\sigma,\tau,\bar\sigma,\bar\tau)\, O(T^{-2})
\end{aligned}$$ for some $f\in L^2$, where $O(T^{-2})$ does not depend on $\sigma,\tau,\bar\sigma,\bar\tau$. All other mixed covariances of this type can be treated similarly. Since additionally $var(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau))=f(\sigma,\tau)\, O(1)$ for some $f\in L^2$, $$\begin{aligned}
&var(M_T) \\=&\frac{1}{b^3T^2}\sum_{t_1,t_2=-N}^N\int_{-\pi}^\pi W^2\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\, d\lambda \int_{-\pi}^\pi W^2\left(\frac{\bar\lambda-\lambda_{t_2}}{b}\right)d\bar\lambda\\
&\times \int_0^1\int_0^1\int_0^1\int_0^1
cov\left(|\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\,-\, \widehat p_{Y,\lambda_{t_1}}^c(\sigma,\tau)|^2,|\widehat p_{X, \lambda_{t_2}}^c(\bar \sigma,\bar\tau)\,-\, \widehat p_{Y, \lambda_{t_2}}^c(\bar\sigma,\bar\tau)|^2\right)d\sigma\,d\tau\,d\bar\sigma\,d\bar\tau\\
=&o(1)\,+\, \frac{1}{b^3T^2}\sum_{t_1,t_2=-N}^N\int_{-\pi}^\pi W^2\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\, d\lambda \int_{-\pi}^\pi W^2\left(\frac{\bar\lambda-\lambda_{t_2}}{b}\right)d\bar\lambda\\
&\times \int_0^1\int_0^1\int_0^1\int_0^1
cov\left(|\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)|^2,|\widehat p_{X, \lambda_{t_2}}^c(\bar \sigma,\bar\tau)|^2\right)\\
&\qquad
+cov\left(|\widehat p_{Y, \lambda_{t_1}}^c(\sigma,\tau)|^2,|\widehat p_{Y, \lambda_{t_2}}^c(\bar \sigma,\bar\tau)|^2\right)
d\sigma\,d\tau\,d\bar\sigma\,d\bar\tau.
\end{aligned}$$ We consider only the part depending on $X_1,\dots, X_T$ in more detail, since the remaining part (depending on $Y_1,\dots, Y_T$) behaves similarly. From Lemma \[l.cumulant\] and the calculations showing $EM_T\to 0$ we get $$\begin{aligned}
&\frac{1}{b^3T^2}\sum_{t_1,t_2=-N}^N\int_{-\pi}^\pi W^2\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\, d\lambda \int_{-\pi}^\pi W^2\left(\frac{\bar\lambda-\lambda_{t_2}}{b}\right)d\bar\lambda\\
&\quad\times \int_0^1\int_0^1\int_0^1\int_0^1
cov\left(|\widehat p_{X,\lambda_{t_1}}^c(\sigma,\tau)|^2,|\widehat p_{X, \lambda_{t_2}}^c(\bar \sigma,\bar\tau)|^2\right)
d\sigma\,d\tau\,d\bar\sigma\,d\bar\tau\\
&=\frac{1}{b^3T^2}\sum_{t_1,t_2=-N}^N\int_{-\pi}^\pi W^2\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\, d\lambda \int_{-\pi}^\pi W^2\left(\frac{\bar\lambda-\lambda_{t_2}}{b}\right)d\bar\lambda\\
&\quad\times \int_0^1\int_0^1\int_0^1\int_0^1
cum\left(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau),\widehat p_{X, \lambda_{t_2}}^c(\bar \sigma,\bar\tau),\widehat p_{X, -\lambda_{t_1}}^c(\sigma,\tau),\widehat p_{X, -\lambda_{t_2}}^c(\bar \sigma,\bar\tau)\right)\\
&\quad + cov\left(\widehat p_{X, \lambda_{t_1}}(\sigma,\tau),\widehat p_{X, -\lambda_{t_2}}(\bar \sigma,\bar\tau)\right) \,cov\left(\widehat p_{X, \lambda_{t_2}}(\bar \sigma,\bar\tau),\widehat p_{X, \lambda_{t_1}}(\sigma,\tau)\right)\\
&\quad + cov\left(\widehat p_{X, \lambda_{t_1}}(\sigma,\tau),\widehat p_{X, \lambda_{t_2}}^c(\bar \sigma,\bar\tau)\right)\,cov\left(\widehat p_{X, \lambda_{t_2}}(\bar \sigma,\bar\tau),\,\widehat p_{X, \lambda_{t_1}}(\sigma,\tau)\right)
d\sigma\,d\tau\,d\bar\sigma\,d\bar\tau\\
& = o(1).
\end{aligned}$$

\[l.var\] Suppose that the assumptions of Theorem \[t.asy-test\] hold true. Then, $$var(\sqrt b T\, L_T)\ninfty\theta_0$$ for $\theta_0$ defined in Theorem \[t.asy-test\], where $$\label{eq.LT}
\begin{aligned}
L_T&= \frac{1}{b^{2}T^2}\int_{-\pi}^\pi\int_0^1\int_0^1\sum_{t_1,t_2=-N,\, t_1\neq t_2}^N \, W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\,W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,
\\
&\quad\times\left(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\,-\, \widehat p_{Y, \lambda_{t_1}}^c(\sigma,\tau)\right)\,\left(\widehat p_{X, -\lambda_{t_1}}^c(\sigma,\tau)\,-\, \widehat p_{Y, -\lambda_{t_2}}^c(\sigma,\tau)\right)\, d\sigma\, d\tau\, d\lambda
\end{aligned}$$ with $\widehat p_{X, \lambda}^c=\widehat p_{X, \lambda}-E\widehat p_{X, \lambda}$ and $\widehat p_{Y, \lambda}^c=\widehat p_{Y, \lambda}-E\widehat p_{Y, \lambda}$.

By independence of the time series $(X_t)_t$ and $(Y_t)_t$, we get $$\label{eq.vt}
\begin{aligned}
var((\sqrt b &T\, L_T)=
\frac{1}{b^{3}T^2}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\sum_{\substack{t_3,t_4=-N \\ t_3\neq t_4}}^N\int_{-\pi}^\pi \, W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,d\lambda\\
&\quad\times\;\int_{-\pi}^\pi W\left(\frac{\bar\lambda-\lambda_{t_3}}{b}\right)W\left(\frac{\bar\lambda-\lambda_{t_4}}{b}\right)d\bar\lambda\,
\\
&\quad\times\int_0^1\int_0^1\int_0^1\int_0^1 cov\Bigg(\left[\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\,-\, \widehat p_{Y, \lambda_{t_1}}^c(\sigma,\tau)\right]\,\left[\widehat p_{X, -\lambda_{t_2}}^c(\sigma,\tau)\,-\, \widehat p_{Y, -\lambda_{t_2}}^c(\sigma,\tau)\right],\,
\\
&\quad
\left[\widehat p_{X, \lambda_{t_3}}^c(\sigma,\tau)\,-\, \widehat p_{Y, \lambda_{t_3}}^c(\sigma,\tau)\right]\,\left[\widehat p_{X, -\lambda_{t_4}}^c(\bar\sigma,\bar\tau)\,-\, \widehat p_{Y, -\lambda_{t_4}}^c(\bar\sigma,\bar\tau)\right]\Bigg) \, d\sigma\, d\tau\,d\bar\sigma\, d\bar\tau \\
&=\frac{1}{b^{3}T^2}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\sum_{\substack{t_3,t_4=-N \\ t_3\neq t_4}}^N\int_{-\pi}^\pi \, W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,d\lambda\;
\\
&\quad\times\int_{-\pi}^\pi W\left(\frac{\bar\lambda-\lambda_{t_3}}{b}\right)W\left(\frac{\bar\lambda-\lambda_{t_4}}{b}\right) d\bar\lambda\,
\\
&\quad\times\int_0^1\int_0^1\int_0^1\int_0^1 cov\Bigg(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\widehat p_{X, -\lambda_{t_2}}^c(\sigma,\tau),\,
\widehat p_{X,\lambda_{t_3}}^c(\bar\sigma,\bar\tau)\widehat p_{X, -\lambda_{t_4}}^c(\bar\sigma,\bar\tau)\Bigg)\\
&\quad+cov\Bigg(\widehat p_{Y, \lambda_{t_1}}^c(\sigma,\tau)\widehat p_{Y, -\lambda_{t_2}}^c(\sigma,\tau),\,
\widehat p_{Y, \lambda_{t_3}}^c(\bar\sigma,\bar\tau)\widehat p_{Y, -\lambda_{t_4}}^c(\bar\sigma,\bar\tau)\Bigg)\\
&\quad+4cov\Bigg(\widehat p_{X, \lambda_{t_1}}(\sigma,\tau)\,,\,\widehat p_{X, \lambda_{t_3}}(\bar\sigma,\bar\tau)\Bigg)
\,cov\Bigg(\widehat p_{Y, -\lambda_{t_2}}(\sigma,\tau),\,\widehat p_{Y, -\lambda_{t_4}}(\bar\sigma,\bar\tau)\Bigg) \, d\sigma\, d\tau\,d\bar\sigma\, d\bar\tau .
\end{aligned}$$ We consider the first summand on the r.h.s. in more detail. Invoking Lemma \[l.cumulant\], we obtain $$\begin{aligned}
&\frac{1}{b^{3}T^2}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\sum_{\substack{t_3,t_4=-N \\ t_3\neq t_4}}^N\int_{-\pi}^\pi \, W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,d\lambda\;
\\
&\quad\times\int_{-\pi}^\pi W\left(\frac{\lambda-\lambda_{t_3}}{b}\right)W\left(\frac{\bar\lambda-\lambda_{t_4}}{b}\right) d\bar\lambda\,
\\
&\quad\times\int_0^1\int_0^1\int_0^1\int_0^1 cov\Bigg(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\widehat p_{X, -\lambda_{t_2}}^c(\sigma,\tau),\,
\widehat p_{X, \lambda_{t_3}}^c(\bar\sigma,\bar\tau)\widehat p_{X, \lambda_{t_4}}^c(\bar\sigma,\bar\tau)\Bigg) \, d\sigma\, d\tau\,d\bar\sigma\, d\bar\tau\\
&=  o(1)\,+\,\frac{1}{bT^2}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\sum_{\substack{t_3,t_4=-N \\ t_3\neq t_4}}^N\int_{-\pi/b}^{\pi/b} \,W\left(\lambda-\frac{\lambda_{t_1}}{b}\right)W\left(\lambda-\frac{\lambda_{t_2}}{b}\right)\,d\lambda\;\\
&\quad\times\int_{-\pi/b}^{\pi/b} W\left(\bar\lambda-\frac{\lambda_{t_3}}{b}\right)W\left(\bar\lambda-\frac{\lambda_{t_4}}{b}\right) d\bar\lambda\,
\\
&\quad\times\int_0^1\int_0^1\int_0^1\int_0^1 cov\Big(\widehat p_{X, \lambda_{t_1}}(\sigma,\tau),\, \widehat p_{X, \lambda_{t_3}}(\bar\sigma,\bar\tau)\Big)\, cov\Big(\widehat p_{X, -\lambda_{t_2}}(\sigma,\tau),\, \widehat p_{X, -\lambda_{t_4}}(\bar\sigma,\bar\tau)\Big)\\
&\quad + cov\Big(\widehat p_{X,\lambda_{t_1}}(\sigma,\tau),\, \widehat p_{X, -\lambda_{t_4}}(\bar\sigma,\bar\tau)\Big)\, cov\Big(\widehat p_{X, -\lambda_{t_2}}(\sigma,\tau),\, \widehat p_{X, \lambda_{t_3}}(\bar\sigma,\bar\tau)\Big) \, d\sigma\, d\tau\,d\bar\sigma\, d\bar\tau .\\
\end{aligned}$$ Similar to the proof of Lemma \[l.mean\], the latter can be bounded by $$\label{eq.step2}
\begin{aligned}
&o(1) \,+\, \frac{{2}}{bT^2}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\left\{\int_{-\pi/b}^{\pi/b} \, W\left(\lambda-\frac{\lambda_{t_1}}{b}\right)W\left(\lambda-\frac{\lambda_{t_2}}{b}\right)\,d\lambda\right\}^2\,
\\
&\quad\times\int_0^1\int_0^1\int_0^1\int_0^1 cov\Big(\widehat p_{X, \lambda_{t_1}}(\sigma,\tau),\, \widehat p_{X, \lambda_{t_1}}(\bar\sigma,\bar\tau)\Big)\, cov\Big(\widehat p_{X, -\lambda_{t_2}}(\sigma,\tau),\, \widehat p_{X, -\lambda_{t_2}}(\bar\sigma,\bar\tau)\Big) \, d\sigma\, d\tau\,d\bar\sigma\, d\bar\tau\\
&\ninfty
\frac{1}{2\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \int_0^1\int_0^1\int_0^1\int_0^1f_{X,\lambda}^2(\sigma,\bar\sigma)\, f_{X,\lambda}^2(\tau,\bar \tau)\,\,d\sigma\,d\bar\sigma\, d\tau\,d\bar\tau d\lambda.
\end{aligned}$$ Proceeding similarly for the remaining summands of , we obtain that $var((\sqrt b T\, L_T)\ninfty\theta_0^2$.

We can now turn to the proof of the theorem itself:

First note that, under $\mathcal H_0$, $${\mathcal U_T}=\int_{0}^\pi\|\widehat{\mathcal F}_{X,\lambda}- \mathcal F_{X,\lambda}+\mathcal F_{Y,\lambda}-\widehat{\mathcal F}_{Y,\lambda}\|_{HS}^2\, d\lambda.$$ In view of Lemma \[l.centering\], it remains to verify that $$\widetilde{\mathcal U}_T:=\sqrt b T\int_{-\pi}^\pi\|\widehat{\mathcal F}_{X,\lambda}-E\widehat{\mathcal F}_{X,\lambda}+ E\widehat{\mathcal F}_{Y,\lambda}-\widehat{\mathcal F}_{Y,\lambda}\|_{HS}^2\, d\lambda- b^{-1/2} \mu_0\stackrel{d}{\longrightarrow} Z.$$ To this end, we split up $$\widetilde{\mathcal U}_T\,:=\,\sqrt b T( M_T\,+\, L_T)\,-\,  b^{-1/2} \mu_0,$$ where $M_T$ and $L_T$ are defined as in Lemma \[l.mean\] and Lemma \[l.var\], respectively. In view of Lemma \[l.mean\], it remains to show that $\sqrt b T L_T\stackrel{d}{\longrightarrow} Z.$ We proceed in several steps. First we derive a CLT for the fidis $$\label{eq.fidis}
\begin{aligned}
\widetilde L_T(
{\sigma_l,\tau_l})&:=\sum_{t_1,t_2=-N,\, t_1\neq t_2}^N\frac{1}{b^{3/2}T}\int_{-\pi}^\pi \, W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\,W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,\,d\lambda
\\
&\quad\times\left(\widehat p_{X, \lambda_{t_1}}^c(\sigma_l,\tau_l)\,-\, \widehat p_{Y, 2\pi t_1/T}^c(\sigma_l,\tau_l)\right)\,\left(\widehat p_{X, -\lambda_{t_2}}^c(\sigma_l,\tau_l)\,-\, \widehat p_{Y, -\lambda_{t_2}}^c(\sigma_l,\tau_l)\right)\\
\end{aligned}$$ for $l=1,\dots, L,\; L\in\N$, in Steps 1 and 2 and prove tightness in $L^2$ in Step 3. The continuous mapping theorem then yields the desired result.\
[ *(A simplifying approximation of $\widetilde L_T$)*]{} We will show that $$\begin{aligned}
|\widetilde L_T(\sigma,\tau)-\widehat L_T(\sigma,\tau)|&= \,o_P(1)\, f(\sigma,\tau)
\end{aligned}$$ for some $f\in L^2$ with $$\label{eq.H}
\begin{aligned}
\widehat L_T(\sigma,\tau)&=
\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N \,\frac{1}{b^{3/2}T}\, \int_{-\pi}^\pi W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\,W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,d\lambda\,
\\
&\quad\times\left(\left|\sum_{j=-\infty}^\infty a_j\, e^{-i j\lambda_{t_1}}\right|^2\widehat p_{\varepsilon, \lambda_{t_1}}^c\sigma,\tau)\,-\, \left|\sum_{j=-\infty}^\infty b_j\, e^{-i j\lambda_{t_1}}\right|^2\widehat p_{e, \lambda_{t_1}}^c(\sigma,\tau)\right)\\
&\quad\times\left(\left|\sum_{j=-\infty}^\infty a_j\, e^{-i j\lambda_{t_2}}\right|^2\widehat p_{\varepsilon, -\lambda_{t_2}}^c(\sigma,\tau)\,-\, \left|\sum_{j=-\infty}^\infty b_j\, e^{-ij \lambda_{t_2}}\right|^2\widehat p_{e, -\lambda_{t_2}}^c(\sigma,\tau)\right)\\
&=:\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N  H_{t_1,t_2}(\sigma,\tau).
\end{aligned}$$ From Lemma \[l.bd-thm10.3.1\], we obtain with $R_{T,\lambda}^{X,c}=R_{T,\lambda}^{X}-ER_{T,\lambda}^{X}$ (and similarly for $R^{Y,c}_{T,\lambda}$) $$\begin{aligned}
&|\widetilde L_T(\sigma,\tau)-\widehat L_T(\sigma,\tau)|\\
&=o_P(1) \, \widetilde f(\sigma,\tau)+\frac{1}{b^{3/2}T}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N \,\int_{-\pi}^\pi W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\,W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,\, d\lambda
\\
&\quad\times\Bigg\{\left(\left|\sum_{j=-\infty}^\infty a_j\, e^{-i j \lambda_{t_1}}\right|^2\widehat p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\,-\, \left|\sum_{j=-\infty}^\infty b_j\, -e^{-i\pi j \lambda_{t_1}}\right|^2\widehat p_{e, \lambda_{t_1}}^c(\sigma,\tau)\right)\,\\
&\times (R_{T,-\lambda_{t_2}}^{X,c}(\sigma,\tau)\,+\, R_{T,-\lambda_{t_2}}^{Y,c}(\sigma,\tau))+(R_{T,\lambda_{t_1}}^{X,c}(\sigma,\tau)\,+\, R_{T,\lambda_{t_1}}^{Y,c}(\sigma,\tau))\\
&\times
\left(\left|\sum_{j=-\infty}^\infty a_j\, e^{-i j \lambda_{t_2}}\right|^2\widehat p_{\varepsilon, -\lambda_{t_2}}^c(\sigma,\tau)\,-\, \left|\sum_{j=-\infty}^\infty b_j\, e^{-i  j \lambda_{t_2}}\right|^2\widehat p_{e, -\lambda_{t_1}}^c(\sigma,\tau)\right)\,\Bigg\}
\end{aligned}$$ for some $\widetilde f\in L^2$. In order to show that the r.h.s. vanishes asymptotically, we consider exemplarily $$\begin{aligned}
&\frac{1}{b^{3/2}T}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N \,\int_{-\pi}^\pi W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\,W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\, d\lambda \\
& \times\,\left|\sum_{j=-\infty}^\infty a_j\, e^{-ij\lambda_{t_1}}\right|^2\widehat p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\, R_{T,-\lambda_{t_2}}^{X,c}(\sigma,\tau)
\end{aligned}$$ and show asymptotic negligibility of its absolute second moment $$\begin{aligned}
&\frac{1}{bT^2}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N \sum_{\substack{t_3,t_4=-N \\ t_3\neq t_4}}^N 
\int_{-\pi/b}^{\pi/b}\, W\left(\lambda-\frac{\lambda_{t_1}}{b}\right)\,W\left(\lambda-\frac{\lambda_{t_2}}{b}\right)d\lambda\\
&\times 
\int_{-\pi/b}^{\pi/b}\, W\left(\bar\lambda-\frac{\lambda_{t_3}}{b}\right)\,W\left(\bar\lambda-\frac{\lambda_{t_4}}{b}\right)d\bar\lambda\,\left|\sum_{j=-\infty}^\infty a_j\, e^{-ij\lambda_{t_1}}\right|^2\,\left|\sum_{j=-\infty}^\infty a_j\, e^{-ij\lambda_{t_3}}\right|^2 \\
&\times E\left(\widehat p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\,\widehat p_{\varepsilon, -\lambda_{t_3}}^c(\sigma,\tau)\, R_{T,-\lambda_{t_2}}^{X,c}(\sigma,\tau)\,R_{T,\lambda_{t_4}}^{X,c}(\sigma,\tau)\right).
\end{aligned}$$ If $M:=\#\{|t_1|,\,|t_2|,\, |t_3|,\,|t_4|\}<4$ this can be easily deduced from $T^2\, \sup_{\lambda\in[0,\pi]}\,E |R_{T,\lambda_t}^X(\sigma,\tau)|^4\in L^2$ and $T\, E(\widehat p_{\varepsilon, \lambda_s}^c(\sigma,\tau)\,\widehat p_{\varepsilon, \lambda_t}^c(\sigma,\tau))\,\in L^2$ for $\lambda_{s}\neq\lambda_t\in [0,\pi]$ using a cumulant type of argument again. Hence, it suffices to show that $$T^2\,\left|E\left(\widehat p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\,\widehat p_{\varepsilon, -\lambda_{t_3}}^c(\sigma,\tau)\, R_{T,-\lambda_{t_2}}^{X,c}(\sigma,\tau)\,R_{T,\lambda_{t_4}}^{X,c}(\sigma,\tau)\right)\, \right|\in L^2 
$$ if $M=4$. We obtain $$\begin{aligned}
&\left|E\left(\widehat p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\,\widehat p_{\varepsilon, -\lambda_{t_3}}^c(\sigma,\tau)\, R_{T,-\lambda_{t_2}}^{X,c}(\sigma,\tau)\,R_{T,\lambda_{t_4}}^{X,c}(\sigma,\tau)\right)\,\right|\\
&=\frac{1}{T^4}\Bigg|E\sum_{s_1,\dots, s_4=1}^T (\varepsilon_{s_1}\varepsilon_{s_2})^c(\sigma,\tau)(\varepsilon_{s_3}\varepsilon_{s_4})^c(\sigma,\tau)e^{i2\pi(s_1-s_2)t_1/T}e^{-i2\pi(s_3-s_4)t_3/T}\\
&\quad\times
\sum_{k_1,k_2=-\infty}^\infty a_{k_1} a_{k_2} \times \Bigg\{\sum_{s_6=1}^{T} \\
& \left[\sum_{s_5=1-k_1}^{T-k_1} e^{i(s_5+k_1-s_6 -k_2)2\pi t_2/T}(\varepsilon_{s_5}\varepsilon_{s_6})^c(\sigma,\tau) - 
\sum_{s_5=1}^{T} e^{i(s_5+k_1-s_6 -k_2)2\pi t_2/T}(\varepsilon_{s_5}\varepsilon_{s_6})^c(\sigma,\tau)\right] + \sum_{s_5=1-k_1}^{T-k_1}\\
&\left[\sum_{s_6=1-k_2}^{T-k_2} e^{i(s_5+k_1-s_6 -k_2)2\pi t_2/T}(\varepsilon_{s_5}\varepsilon_{s_6})^c(\sigma,\tau)- \sum_{s_6=1}^{T} e^{i(s_5+k_1-s_6 -k_2)2\pi t_2/T}(\varepsilon_{s_5}\varepsilon_{s_6})^c(\sigma,\tau)\right]\Bigg\}\\
&\quad\times
\sum_{k_3,k_4=-\infty}^\infty a_{k_3} a_{k_4}\Bigg\{\sum_{s_8=1}^{T} \\
& \left[\sum_{s_7=1-k_3}^{T-k_3} e^{i(s_7+k_3-s_8 -k_4)2\pi t_4/T}(\varepsilon_{s_7}\varepsilon_{s_8})^c(\sigma,\tau)- \sum_{s_7=1}^{T} e^{i(s_7+k_3-s_8 -k_4)2\pi t_4/T}(\varepsilon_{s_7}\varepsilon_{s_8})^c(\sigma,\tau)\right]+\sum_{s_7=1-k_3}^{T-k_3}\\
&\left[\sum_{s_8=1-k_4}^{T-k_4} e^{i(s_7+k_3-s_8 -k_4)2\pi t_4/T}(\varepsilon_{s_7}\varepsilon_{s_8})^c(\sigma,\tau)- \sum_{s_8=1}^{T} e^{i(s_7+k_3-s_8 -k_4)2\pi t_4/T}(\varepsilon_{s_7}\varepsilon_{s_8})^c(\sigma,\tau)\right]\Bigg\}
\Bigg|.\end{aligned}$$

For symmetry reasons, it suffices to investigate $$\begin{aligned}
&\frac{1}{T^4}\Bigg|E\sum_{s_1,\dots, s_4=1}^T (\varepsilon_{s_1}\varepsilon_{s_2})^c(\sigma,\tau)(\varepsilon_{s_3}\varepsilon_{s_4})^c(\sigma,\tau)e^{i2\pi(s_1-s_2)t_1/T}e^{-i2\pi(s_3-s_4)t_3/T} \nonumber \\
&\quad\times
\sum_{k_1,k_2=-\infty}^\infty a_{k_1} a_{k_2}\sum_{s_6=1}^{T}\Big[\sum_{s_5=1-k_1}^{T-k_1} e^{i(s_5+k_1-s_6 -k_2)2\pi t_2/T}(\varepsilon_{s_5}\varepsilon_{s_6})^c(\sigma,\tau) \nonumber \\
& \ \ \ \ - \sum_{s_5=1}^{T} e^{i(s_5+k_1-s_6 -k_2)2\pi t_2/T}(\varepsilon_{s_5}\varepsilon_{s_6})^c(\sigma,\tau)\Big] \nonumber\\
&\quad\times
\sum_{k_3,k_4=-\infty}^\infty a_{k_3} a_{k_4}\sum_{s_8=1}^{T}\Big[\sum_{s_7=1-k_3}^{T-k_3} e^{i(s_7+k_3-s_8 -k_4)2\pi t_4/T}(\varepsilon_{s_7}\varepsilon_{s_8})^c(\sigma,\tau) \nonumber \\
& \ \ \ \ - \sum_{s_7=1}^{T} e^{i(s_7+k_3-s_8 -k_4)2\pi t_4/T}(\varepsilon_{s_7}\varepsilon_{s_8})^c(\sigma,\tau)\Big]\Bigg|.
\label{eq.hilf}\end{aligned}$$ Now, we discuss terms with different numbers of equal indices $s_1,\dots, s_8$ separately. Straightforward calculations show that the desired order $O(T^{-2})$ is reached if the number of different indices is not equal to 3 or 4. For the remaining case of 3 and 4 different indices, we consider exemplarily the situation $s_1=s_3,~s_2=s_4$, $s_5=s_7$, $s_6=s_8$. Using that $\sum_{s=1}^T e^{2\pi i \lambda_s}=0$ for any Fourier frequency $\lambda\notin T\Z$, the corresponding summands of can be bounded by $$\begin{aligned}
&\frac{\hat f(\sigma,\tau)}{T^4}\sum_{k_1,\dots, k_4=-\infty}^\infty |a_{k_1}\cdots a_{k_4}|
\sum_{s_5=-\infty}^\infty|\I(s_5\in \{1-k_1,\dots,T-k_1\})-\I(s_5\in \{1,\dots,T\})|\\
&\times\Bigg|\sum_{\substack{s_1,s_2=1 \\ s_1\neq s_2 \\ s_1,s_2\neq s_5}}^T
e^{i2\pi(s_1-s_2)(t_1-t_3)/T}\sum_{\substack{s_6=1 \\ s_6\notin\{s_1,s_2, s_5\}}}^T e^{-i2\pi s_6(t_2+t_4)/T}
\Bigg|\\
&\leq \frac{\hat f(\sigma,\tau)}{T^2}\sum_{k_1,\dots, k_4=-\infty}^\infty |a_{k_1}\cdots a_{k_4}|\min\{T, k_1\}\\
&=O(T^{-2})\, \hat f(\sigma,\tau)
\end{aligned}$$ for some $\hat f\in L^2$.\
[ *(CLT for fidis)*]{} Due to the previous step and the Cramér-Wold device, it suffices to verify $$\widehat W_T:=\sum_{l=1}^L d_l\,\widehat L_T(\sigma_l,\tau_l)\stackrel{d}{\longrightarrow} \hat Z\sim \mathcal N(0, V_{\widehat W})$$ with $$V_{\widehat W}=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\,\sum_{l_1,l_2=1}^L d_{l_1}\, d_{l_2} \int_{-\pi}^\pi f_{X,\lambda}^2(\sigma_{l_1}, \tau_{l_1})f_{X,\lambda}^2(\sigma_{l_2}, \tau_{l_2}) \, d\lambda$$ for any $d_1,\dots, d_L\in \R$ in order to show the desired weak convergence result for $(\widetilde L_T(\sigma_1,\tau_1),$ $ \dots, $ $ \widetilde L_T(\sigma_L,\tau_L)) $ defined in . To this end, first note that $$\widehat L_T(\sigma, \tau)=\sum_{\substack{t_1,t_2=1 \\ t_1\neq t_2}}^N \widetilde H_{t_1,t_2}(\sigma, \tau)+\sum_{t_1=-N}^N H_{t_1,0}(\sigma, \tau)+\sum_{t_2=-N}^N H_{0,t_2}(\sigma, \tau)$$ with $H_{t_1,t_2}$ defined in and $
\widetilde H_{t_1,t_2}:=H_{t_1,t_2}+ H_{-t_1,t_2}+ H_{t_1,-t_2}+ H_{-t_1,-t_2}. 
$ Straightforward calculations yield that $\sum_{t_1=-N}^N H_{t_1,0}+\sum_{t_2=-N}^N H_{0,t_2}=o_P(1)$ in $L^2$. Now, we apply Theorem 2.1 of de Jong (1987) to $$\widetilde W_T=\sum_{\substack{t_1, t_2=1 \\ t_1\neq t_2}}^N \sum_{l=1}^L d_l\, \widetilde H_{t_1,t_2}(\sigma_l,\tau_l)$$ with $$\mathbb X_t=(p_{\varepsilon, \lambda_t}^c(\sigma_1,\tau_1),\dots,p_{\varepsilon, \lambda_t}^c(\sigma_L,\tau_L),p_{e, \lambda_t}^c(\sigma_1,\tau_1),\dots, p_{e,\lambda_t}^c(\sigma_L,\tau_L))$$ in their notation. First, note that for fixed $\sigma,\tau\in [0,1]$ and under the assumption of normally distributed innovations, the corresponding periodograms evaluated at different Fourier frequencies are stochastically independent which implies independence of $\mathbb X_{1},\dots, \mathbb X_N$. Moreover, this yields $E(\widetilde H_{t_1,t_2}\mid \mathbb X_{t_1})=E(\widetilde H_{t_1,t_2}\mid \mathbb X_{t_2})=0~a.s.$ for $t_1\neq t_2$ which implies that $\widetilde W_T$ is clean (see Definition 2.1 in de Jong (1987)). It remains to check conditions (a) and (b) of Theorem 2.1 of de Jong (1987). It can be shown in analogy to Lemma \[l.var\] that $var(\widetilde W_T)$ converges to the finite constant $V_{\widetilde W}$. Subsequently, we only consider the non-trivial case of $V_{\widetilde W}>0$. For condition (a), it remains to verify that $$\max_{t_1\in\{1,\dots, N\}}\sum_{\substack{t_2=1 \\ t_2\neq t_1}}^N var\left(\sum_{l=1}^L d_l\,\widetilde H_{t_1,t_2}(\sigma_l,\tau_l)\right)=o(1).$$ This is an immediate consequence of $var(H_{t_1,t_2}(\sigma,\tau))=0$ for $|t_1-t_2|>bT$ and $$\begin{aligned}
&var(H_{t_1,t_2}(\sigma,\tau))\\
&\leq O\left(\frac{1}{b\,T^2}\right)\,\big( E|p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\, p_{\varepsilon, \lambda_{t_2}}^c(\sigma,\tau)|^2+ 
 E|p_{e, \lambda_{t_1}}^c(\sigma,\tau)\, p_{e, \lambda_{t_2}}^c(\sigma,\tau)|^2\\
&\quad+
 E|p_{\varepsilon, \lambda_{t_1}}^c(\sigma,\tau)\, p_{e, \lambda_{t_2}}^c(\sigma,\tau)|^2+
 E|p_{e, \lambda_{t_1}}^c(\sigma,\tau)\, p_{\varepsilon, \lambda_{t_2}}^c(\sigma,\tau)|^2\big) \\
&=o\left(\frac{1}{bT}\right)
\end{aligned}$$ for $|t_1-t_2|\leq bT$ and any $\sigma,\tau\in [0,1]$. Finally, we have to check assumption (b) of Theorem 2.1 of de Jong (1987), i.e. $$E\widetilde W_T^4\ninfty 3 V_{\widetilde W}^2.$$ To this end, we argue that $E\widetilde W_T^2\ninfty V_{\widetilde W}$ and that the forth-order cumulant of $\widetilde W_T$ vanishes asymptotically.\
[ *(Tightness)*]{} Tightness can be proven along the lines of Lemma 16 of Cerovecki and Hörmann (2017) if we can verify validity of their Lemma 14 with $$p_j^{(T)}=E|\langle \widetilde L_T,\, e_j\rangle|^2 \geq 0$$ where $(e_j)_{j\in \N}$ be an arbitrary complete ONB in $L_2^\R([0,1]^2,\mu)$. As in Lemma \[l.var\], we obtain $$\begin{aligned}
\lim_{T\to \infty} p_j^{(T)}
&=\lim_{T\to\infty} E\Bigg|\int_0^1\int_0^1\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\frac{1}{b^{3/2}T}\int_{-\pi}^\pi \, W\left(\frac{\lambda-\lambda_{t_1}}{b}\right)\,W\left(\frac{\lambda-\lambda_{t_2}}{b}\right)\,\,d\lambda
\\
&\quad\times\left(\widehat p_{X, \lambda_{t_1}}^c(\sigma,\tau)\,-\, \widehat p_{Y, \lambda_{t_1}}^c(\sigma,\tau)\right)\,\left(\widehat p_{X, -\lambda_{t_2}}^c(\sigma,\tau)\,-\, \widehat p_{Y, -\lambda_{t_2}}^c(\sigma,\tau)\right)\, e_j(\sigma, \tau)\, d\sigma\, d\tau \Bigg|^2\\
&=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \left\{\int_0^1\int_0^1f_{X,\lambda}^2(\sigma, \tau)\,e_j(\sigma,\tau)\,d\sigma\, d\tau\right\}^2\, d\lambda\\
&=: p_j^{(0)}\geq 0.
\end{aligned}$$ Thus, conditions (a) and (b) in Lemma 14 of Cerovecki and Hörmann (2017) are satisfied. Moreover, Parseval’s identity gives $$\begin{aligned}
\sum_{j=1}^\infty p_j^{(0)}
&=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi\sum_{j=1}^\infty \langle f^2_{X,\lambda},\, e_j\rangle^2 \, d\lambda\\
&=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \int_0^1\int_0^1 f_{X,\lambda}^4(\sigma,\tau)\,d\sigma\, d\tau\, d\lambda<\infty,
\end{aligned}$$ which yields their condition (c). Their assumption (d) is an immediate consequence of $$\begin{aligned}
\lim_{T\to\infty}\sum_{j=1}^\infty p_j^{(T)}&=\lim_{T\to\infty}E\left(\int_0^1\int_0^1 \widetilde L_T^2(\sigma,\tau)\, d\sigma\, d\tau\right)^2\\
&=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \int_0^1\int_0^1 f_{X,\lambda}^4(\sigma,\tau)\,d\sigma\, d\tau\, d\lambda,
\end{aligned}$$ which can be verified along the lines of the proof of Lemma \[l.var\]. Moreover, it follows from that $\sum_j p_j^{(T)}<\infty \;\forall T$, i.e., condition (e) of Lemma 14 in Cerovecki and Hörmann (2017) is satisfied, which in turn implies the desired tightness result.

We first note that to simplify notation and calculations, all indices in the sums considered, run in the set $ \{-N, -N+1, \ldots, -1,1,\ldots, N-1, N\}$, where $ N=[(T-1)/2]$; i.e., we ignore the asymptotically negligible zero frequency in the consideration of $ \widehat{f}^\ast_{X,\lambda_t}$ and $\widehat{f}^\ast_{Y,\lambda_t} $ in Step 3 of the bootstrap algorithm. Using the notation $$d^*_{t}(\sigma,\tau)= J^\ast_{X,\lambda_t}(\sigma) \overline{J}^\ast_{X,\lambda_t}(\tau)- J^\ast_{Y,\lambda_t}(\sigma) \overline{J}^\ast_{Y,\lambda_t}(\tau),$$ the bootstrap test statistic $$\mathcal U^\ast_T = \frac{2\pi }{T}\sum_{l=-N}^N \int_0^1\int_0^1 \Big|  \widehat{f}^\ast_{X,\lambda_l}(\sigma,\tau)
-\widehat{f}^\ast_{Y,\lambda_l}(\sigma,\tau) \Big|^2 d\sigma d\tau$$ can be decomposed as $$\begin{aligned}
\mathcal U^\ast_T & = \frac{2\pi}{T^3b^2}\sum_{t_1=-N}^N\sum_{t_2=-N}^N \sum_{l=-N}^N W\Big(\frac{\lambda_l-\lambda_{t_1}}{b} \Big)
W\Big(\frac{\lambda_l-\lambda_{t_2}}{b} \Big) \int_0^1\int_0^1 d^*_{t_1}(\sigma,\tau)\overline{d^*_{t_2}}(\sigma,\tau)d\sigma d\tau \\
& =  \frac{2\pi}{T^3b^2}\sum_{l=-N}^N \sum_{t=-N}^N W^2\Big(\frac{\lambda_l-\lambda_{t}}{b} \Big)
 \int_0^1\int_0^1 |d^{*}_{t}(\sigma,\tau)|^2d\sigma d\tau\\
 & \ \ \ \ +  \frac{2\pi}{T^3b^2}\sum_{\substack{t_1,t_2=-N \\t_1\neq t_2}}^N \sum_{l=-N}^N W\Big(\frac{\lambda_l-\lambda_{t_1}}{b} \Big)
W\Big(\frac{\lambda_l-\lambda_{t_2}}{b} \Big) \int_0^1\int_0^1 d^*_{t_1}(\sigma,\tau)\overline{d^*_{t_2}}(\sigma,\tau)d\sigma d\tau\\
& = M^\ast_T + L^\ast_T,  \end{aligned}$$ with an obvious notation for $ M_T^\ast$ and $ L^\ast_T$. We show that $$\label{eq.M}
\sqrt{b} T M_T^\ast \,-\, b^{-1/2}\,\widetilde{\mu}_0 \stackrel{P}{\rightarrow} 0,$$ and $$\label{eq.L}
\sqrt{b} T  L_T^\ast \stackrel{d}{\rightarrow} \mathcal N(0,\widetilde{\theta}_0^2).$$ Consider (\[eq.M\]) and notice that $$\begin{aligned}
E^\ast(\sqrt{b} T M^\ast_T) & = \frac{4\pi}{T^2b^{3/2}} \sum_{t=-N}^{N}\sum_{l=-N}^{N}W^2\Big(\frac{\lambda_l-\lambda_{t}}{b} \Big)\Big(\int_0^1 \widehat{f}_{\lambda_t}(\tau,\tau) d\tau \Big)^2\\
& = \frac{4\pi}{T^2b^{3/2}} \sum_{t=-N}^{N}\sum_{l=-N}^{N}W^2\Big(\frac{\lambda_l-\lambda_{t}}{b} \Big)\Big(\int_0^1 f_{\lambda_t}(\tau,\tau) d\tau \Big)^2\\
& \ \ \  +  O_P\Big(b^{-1/2}\sup_{\lambda_t \in [0, \pi]} \Big| \int_0^1 \widehat{f}_{\lambda_t}(\tau,\tau) d\tau - \int_0^1 f_{\lambda_t}(\tau,\tau) d\tau \Big| \Big)\end{aligned}$$ and, therefore, $E^\ast(\sqrt{b} T M^\ast_T) - b^{-1/2}\widetilde{\mu}_0 \stackrel{P}{\rightarrow} 0$ under Assumption 2. Furthermore, because of the independence of $ d^*_{t}(\sigma, \tau)$ at different frequencies we get $$\begin{aligned}
Var^\ast(\sqrt{b} T M^\ast_T)  & =\frac{1}{T^4b^3}\sum_{t=-N}^N\sum_{l_1=-N}^{N}\sum_{l_2 =-N}^N W^2\Big(\frac{\lambda_{l_1}-\lambda_t}{b}\Big) W^2\Big(\frac{\lambda_{l_2}-\lambda_t}{b}\Big)\\
&\ \ \ \ \times  \int_0^1\cdots \int_0^1 Var\big(|d^{*}_t(\sigma_1,\tau_1)|^2\cdot |d^{*}_t(\sigma_2,\tau_2)|^2 \big) d\sigma_1 \cdots d\tau_2\\
& = O_P(T^{-1}b^{-1}) \rightarrow 0.\end{aligned}$$

Consider next (\[eq.L\]). Notice first that $ E^\ast( L^\ast_T)=0$ and that $$\begin{aligned}
Var^\ast(\sqrt{b}T L^\ast_T) & = \frac{1}{T^4b^3}\sum_{\substack{t_1,t_2=-N \\ t_1\neq t_2}}^N\sum_{\substack{t_3,t_4=-N \\ t_3\neq t_4}}^N\sum_{l_1,l_2 =-N}^N W\Big(\frac{\lambda_{l_1}-\lambda_{t_1}}{b}\Big)W\Big(\frac{\lambda_{l_1}-\lambda_{t_2}}{b}\Big)\\
&\ \  \times W\Big(\frac{\lambda_{l_2}-\lambda_{t_3}}{b}\Big)W\Big(\frac{\lambda_{l_2}-\lambda_{t_4}}{b}\Big)\int_0^1\cdots \int_0^1 \\
& \ \ \times \Big\{E^\ast\big(d^\ast_{t_1}(\sigma_1,\tau_1)\overline{d^\ast_{t_3}}(\sigma_2,\tau_2)\big) 
E^\ast\big(d^\ast_{t_2}(\sigma_1,\tau_1)\overline{d^\ast_{t_4}}(\sigma_2,\tau_2)\big)\\
& \ \ \ \ \ \ + E^\ast\big(d^\ast_{t_1}(\sigma_1,\tau_1)\overline{d^\ast_{t_4}}(\sigma_2,\tau_2)\big) 
E^\ast\big(d^\ast_{t_2}(\sigma_1,\tau_1)\overline{d^\ast_{t_3}}(\sigma_2,\tau_2)\big)\\
& \ \ \ \ \ \ + cum\big(d^\ast_{t_1}(\sigma_1,\tau_1), \overline{d^\ast_{t_2}}(\sigma_1,\tau_1),d^\ast_{t_3}(\sigma_2,\tau_2),
\overline{d^\ast_{t_4}}(\sigma_2,\tau_2) \big)\Big\}d\sigma_1 \cdots  d\tau_2 \\
& = V_{1,T}^\ast + V_{2,T}^\ast + V_{3,T}^\ast,\end{aligned}$$ with an obvious notation for $V_{i,T}^\ast$, $i=1,2,3$. Using the independence of the two processes as well as the fact that $E^\ast \big(d^\ast_{t_1}(\sigma_1,\tau_1)
\overline{d^\ast_{t_2}}(\sigma_1,\tau_1)\big)=0$ for $ |t_1|\neq |t_2|$ and verifying by straightforward calculations that $$E^\ast \big(d^\ast_{t}(\sigma_1,\tau_1)
\overline{d^\ast_{t}}(\sigma_2,\tau_2)\big)=2\widehat{f}_{\lambda_t}(\sigma_1,\sigma_2)\overline{\widehat{f}}_{\lambda_t}(\tau_1,\tau_2),$$ it follows that $$\begin{aligned}
V^\ast_{1,T} & = \frac{4}{T^4b^3}\sum_{t_1=-N}^{N}\sum_{t_2=-N}^N \sum_{l_1=-N}^N\sum_{l_2=-N}^NW\Big(\frac{\lambda_{l_1}-\lambda_{t_1}}{b}\Big)W\Big(\frac{\lambda_{l_2}-\lambda_{t_1}}{b}\Big)
W\Big(\frac{\lambda_{l_1}-\lambda_{t_2}}{b}\Big)W\Big(\frac{\lambda_{l_2}-\lambda_{t_2}}{b}\Big)\\
& \ \ \times \int_0^1\cdots \int_0^1\widehat{f}_{\lambda_{t_1}}(\sigma_1,\sigma_2)\overline{\widehat{f}}_{\lambda_{t_1}}(\tau_1,\tau_2)
\widehat{f}_{\lambda_{t_2}}(\sigma_1,\sigma_2)\overline{\widehat{f}}_{\lambda_{t_2}}(\tau_1,\tau_2)d\sigma_1 \cdots d\tau_2,\end{aligned}$$ which converges in probability to $$\frac{1}{\pi^2}\int_{-2\pi}^{2\pi}\Big\{\int_{-\pi}^\pi W(u)W(u-x)du\Big\}^2dx \int_{-\pi}^\pi\Big\{ \int_0^1\int_0^1 f^2_\lambda (\sigma,\tau)d\sigma d\tau\Big\}^2 d\lambda,$$ by standard arguments regarding convergence of the Riemann sums involved and using the uniform consistency of the estimators $ \widehat{f}_\lambda$. By the same arguments, the term $V_{2,T}^\ast$ converges to the same limit in probability. Furthermore, using the notation $$w(i,j,k,l) = W\Big(\frac{\lambda_{l_i}-\lambda_{k}}{b}\Big)W\Big(\frac{\lambda_{l_j}-\lambda_{k}}{b}\Big)
W\Big(\frac{\lambda_{l_i}-\lambda_{l}}{b}\Big)W\Big(\frac{\lambda_{l_j}-\lambda_{l}}{b}\Big) ,$$ and the fact that $ d^\ast_{-t}(\tau,\sigma) = d^\ast_t(\sigma, \tau)$, we get by the independence of the random variables $d^\ast_{t}(\tau,\sigma) $ and $d^\ast_{s}(\tau,\sigma) $ for frequencies $ |t|\neq |s|$, that $$\begin{aligned}
V_{3,T}^\ast & = \frac{1}{T^4b^3}\int_0^1\cdots \int_0^1  \sum_{\substack{t_1, t_2=-N \\t_1\neq t_2}}^N\sum_{l_1=-N}^N\sum_{l_2=-N}^N
w( l_1,l_2,t_1,t_2) \\
& \ \ \times cum^\ast\Big(d^\ast_{t_1}(\tau_1,\sigma_1), d^\ast_{t_1}(\tau_2,\sigma_2) , d^\ast_{t_2}(\tau_1,\sigma_1), d^\ast_{t_2}(\tau_2,\sigma_2\Big)d\tau_1 \cdots d\sigma_2\\
& = \frac{1}{T^4b^3}\int_0^1\cdots \int_0^1  \sum_{\substack{t_1, t_2=1 \\ t_1\neq t_2}}^N\sum_{l_1=-N}^N\sum_{l_2=-N}^N\\
& \ \ \times\Big\{w( l_1,l_2,-t_1,-t_2) cum^\ast\big(d^\ast_{t_1}(\sigma_1, \tau_1), d^\ast_{t_1}(\sigma_2,\tau_2) , \bar d^\ast_{t_2}(\sigma_1,\tau_1), \bar d^\ast_{t_2}(\sigma_2,\tau_2)\big)\\
& \ \ \ \ \ \ + w( l_1,l_2,-t_1,t_2) cum^\ast\big(d^\ast_{t_1}(\sigma_1, \tau_1), d^\ast_{t_1}(\sigma_2,\tau_2) , \bar d^\ast_{t_2}(\tau_1,\sigma_1), \bar d^\ast_{t_2}(\tau_2,\sigma_2)\big)\\
& \ \ \ \ \ \ + w( l_1,l_2,t_1,-t_2) cum^\ast\big(d^\ast_{t_1}(\tau_1, \sigma_1), d^\ast_{t_1}(\tau_2,\sigma_2) , \bar d^\ast_{t_2}(\sigma_1,\tau_1),\bar  d^\ast_{t_2}(\sigma_2,\tau_2)\big)\\
& \ \ \ \ \ \ + w( l_1,l_2,t_1,t_2) cum^\ast\big(d^\ast_{t_1}(\tau_1, \sigma_1), d^\ast_{t_1}(\tau_2,\sigma_2) , \bar d^\ast_{t_2}(\tau_1,\sigma_1), \bar d^\ast_{t_2}(\tau_2,\sigma_2)\big)\Big\} d\tau_1 \cdots d\sigma_2\\ 
& =0,\end{aligned}$$ due to the independence of the bootstrap periodogram kernels at different Fourier frequencies. Thus, we conclude that $ Var^\ast(\sqrt{b}T L^\ast_T) \stackrel{P}{\rightarrow} \widetilde{\theta}_0^2$. We next establish the weak convergence of $\sqrt{b} T L^\ast_T $ to $ \mathcal N(0,\widetilde{\theta}_0^2)$. We do this in two steps.\
[ *(Convergence of fidis)*]{} Notice that $$\sqrt{b}T L^\ast_T  = \int_0^1\int_0^1 \sum_{1\leq  t_1<t_2 \leq N}  H^\ast_{t_1,t_2}(\sigma,\tau)  d\sigma d\tau,$$ where $$\label{eq.Hstar}
H^\ast_{t_1,t_2}(\sigma,\tau) = 2\big\{h^\ast_{t_1,t_2}(\sigma,\tau)+ h^\ast_{-t_1,t_2}(\sigma,\tau)+h^\ast_{t_1,-t_2}(\sigma,\tau)+ h^\ast_{-t_1,-t_2}(\sigma,\tau)\big\},$$ and $$h^\ast_{t,s}(\sigma,\tau) = \frac{2\pi}{T^2b^{3/2}} \sum_{\substack{l=-N \\l\neq 0}}^N W \Big(\frac{\lambda_l-\lambda_{t}}{b}\Big) 
W\Big( \frac{\lambda_l-\lambda_{s}}{b}\Big) d_{t}^\ast(\sigma,\tau)\overline{d_{s}^\ast}(\sigma,\tau).$$ Since the quadratic form $\sum_{1\leq t_1 < t_2 \leq N}   H^\ast_{t_1,t_2}(\sigma,\tau)$ is clean, (see Definition 2.1 in de Jong (1987)), we prove convergence of the finite-dimensional vector $ \big(\sum_{1 \leq t_1 < t_2 \leq  N } H^\ast_{t_1,t_2}(\sigma_l,\tau_l), \ l=1,2, \ldots, L\big)$, by verifying Conditions (a) and (b) of Theorem 2.1 of de Jong (1987) for the linear combinations $  \sum_{l=1}^L c_l  \sum_{1\leq t_1 < t_2 \leq N}  H^\ast_{t_1,t_2}(\sigma_l,\tau_l)$ for any $c_l \in \mathbb{R}$.

Evaluating $ \sigma_{t_1,t_2}^2 = E^\ast(\sum_{l=1}^L c_l H^\ast_{t_1,t_2}(\sigma_l,\tau_l))^2$ using (\[eq.Hstar\]), yields the expression $$4\sum_{l_1,l_2=1}^{L} c_{l_1}\,c_{l_2}\,\sum_{r_1\in\{-t_1,t_1\}}\sum_{s_1\in\{-t_2,t_2\}}\sum_{r_2\in\{-t_1,t_1\}}\sum_{s_2\in\{-t_2,t_2\}} cov(h^\ast_{r_1,s_1}(\sigma_{l_1},\tau_{l_1}),h^\ast_{r_2,s_2}(\sigma_{l_2},\tau_{l_2})),$$ where $ t_1 \neq t_2 $. Taking into account the independence of the random variables involved for $ t_1\neq t_2$, the covariance terms in the above sum are very similar with a typical one given, for instance for $r_1=t_1,s_1=t_2,r_2=-t_1, s_2=-t_2 $, by $$\begin{aligned}
\frac{1}{T^4b^3}\sum_{j_1}\sum_{j_2} & W\Big(\frac{\lambda_{j_1}-\lambda_{t_1}}{b}\Big)W\Big(\frac{\lambda_{j_1}-\lambda_{t_2}}{b}\Big)W\Big( \frac{\lambda_{j_2}+\lambda_{t_1}}{b}\Big) W\Big(\frac{\lambda_{j_2}+\lambda_{t_2}}{b} \Big) \\
& \ \  \times \sum_{l_1,l_2=1}^{L} c_{l_1}c_{l_2}\widehat{f}_{\lambda_{t_1}}(\sigma_{l_1},\tau_{l_1})\overline{ \widehat{f}}_{\lambda_{-t_1}}(\sigma_{l_2},\tau_{l_2})  \widehat{f}_{\lambda_{t_2}}(\sigma_{l_1},\tau_{l_1})  \widehat{f}_{\lambda_{-t_2}}(\sigma_{l_2},\tau_{l_2}) \\
&= O_P(T^{-2}b^{-1}), \end{aligned}$$ where the $O_P(T^{-2}b^{-1})$ term is uniform in $ t_1$ and $t_2$. Taking into account that\
$ 0< \sigma^2_T=E^\ast(  \sum_{l=1}^L c_l  \sum_{1\leq t_1 < t_2 \leq N}  H^\ast_{t_1,t_2}(\sigma_l,\tau_l) )^2=O_P(1)$, which follows from the calculations of $Var^\ast(\sqrt{b} T L^\ast_T) $, we get that $$\frac{1}{\sigma^2_T} \max_{1\leq t_1 \leq  N}\sum_{1\leq t_2\leq N} \sigma^2_{t_1,t_2} =  O_P(T^{-1}b^{-1}) \rightarrow 0,$$ as $ T \rightarrow \infty$, which shows that Condition (a) of Theorem 2.1 of de Jong (1987) is satisfied.\
Consider Condition (b). From (\[eq.Hstar\]), the fourth moment of $ \sum_{l=1}^L c_l  \sum_{1\leq t_1 < t_2 \leq N}  H^\ast_{t_1,t_2}(\sigma_l,\tau_l)$ equals $$\begin{aligned}
16& \sum_{l_1,l_2,l_3,l_4=1}^L c_{l_1}c_{l_2}c_{l_3}c_{l_4}\sum_{1\leq t_1 <t_2 \leq N}  \sum_{1\leq r_1 <r_2 \leq N}\sum_{1\leq s_1 <s_2 \leq N}\sum_{1\leq m_1 <m_2 \leq N}
\sum_{j_1 \in\{ -t_1,t_1\} \atop  j_2\in\{-t_2,t_2\}} \sum_{k_1 \in\{ -r_1,r_1\} \atop k_2\in\{-r_2,r_2\}} \\
&\times \sum_{n_1 \in\{ -s_1,s_1\} \atop  n_2\in\{-s_2,s_2\}} 
\sum_{v_1 \in\{ -m_1,m_1\} \atop v_2\in\{-m_2,m_2\}}E^\ast\Big( h^\ast_{j_1,j_2}(\sigma_{l_1},\tau_{l_1})h^\ast_{k_1,k_2}(\sigma_{l_2},\tau_{l_2})h^\ast_{n_1,n_2}(\sigma_{l_3},\tau_{l_3})h^\ast_{v_1,v_2}(\sigma_{l_4},\tau_{l_4})\Big)\end{aligned}$$ where only for the following four cases the expectation term is different from zero: 1) $(j_1,j_2)=(k_1,k_2) \neq (n_1,n_2)=(v_1,v_2)$, 2) $(j_1,j_2)=(n_1,n_2) \neq (k_1,k_2)=(v_1,v_2)$, 3) $(j_1,j_2)=(v_1,v_2) \neq (k_1,k_2)=(n_1,n_2)$ and 4) $(j_1,j_2)=(k_1,k_2) = (n_1,n_2)=(v_1,v_2)$ and where the notation $ (i,j)=(l,k)$ means $ i=l$ and $j=k$. Straightforward calculations show that case 4) vanishes asymptotically while cases 1), 2) and 3) converge to the same limit as $ \sigma^4_T$ converges, from which we conclude that $  E[\sum_{l=1}^L c_l  \sum_{1\leq t_1 < t_2 \leq N}  H^\ast_{t_1,t_2}(\sigma_l,\tau_l)]^4/\sigma^4_T \rightarrow 3$, as $ T \rightarrow \infty$.\
[ *(Tightness)*]{} We show tightness in probability of the sequence $\{\sqrt{b}T L^\ast_T\}_{T\in \N}$ by verifying conditions (a)-(e) of Lemma 14 of Cerovecki and Hörmann (2017) for the sequence $ p_{j,T}^{\ast} = E^\ast|\langle \sqrt{b}T L^\ast_T, e_j\rangle |^2$, where $ (e_j)_{j\in \N}$ is an arbitrary complete orthonormal basis in $ L^\R_2([0,1]^2,\mu)$. As in the proof of $Var^\ast(\sqrt{b}TL^\ast_T)$, we derive that $ \lim_{T\rightarrow \infty}p_{j,T}^{\ast}=p_j^{0}$ in probability, where $$p_j^0=\frac{2}{\pi^2}\int_{-2\pi}^{2\pi}\left\{\int_{-\pi}^\pi W(u)W(u-x)\, du\right\}^2dx\, \int_{-\pi}^\pi \left\{\int_0^1\int_0^1f_{\lambda}^2(\sigma, \tau)e_j(\sigma,\tau)\,\,d\sigma\, d\tau\right\}^2\, d\lambda \geq 0.$$ In other words, conditions (a)-(b) of Lemma 14 in Cerovecki and Hörmann (2017) are satisfied. condition (c) is verified exactly as in proving that $ \sum_{j=1}^\infty p_j^{(0)}=p <\infty$ in Step 5 of the proof of Theorem \[t.asy-test\], while condition (d) follows because, again as in the proof of $Var^\ast(\sqrt{b}TL^\ast_T)$, we get $$\begin{aligned}
\lim_{T\rightarrow \infty}\sum_{j=1}^\infty p^{\ast}_{j,T} & =\lim_{T\rightarrow\infty}E^\ast\Big(\int_0^1\int_0^1 \sum_{1\leq t_<t_2\leq N}H^\ast_{t_1,t_2}(\sigma,\tau) d\sigma d\tau \Big)^2= p,\end{aligned}$$ in probability. Finally, by the construction of bootstrap pseudo-random elements $H^\ast_{t_1,t_2}$ and the convergence above, we get that $ \sum_{j=1}^\infty p_{j,T}^\ast $ is, bounded in probability, uniformly in $T$.

By Assumption 2, we have $ b^{-1/2}\big(\widehat{\mu}_0 -\widetilde{\mu}_0) \stackrel{P}{\rightarrow} 0$ and $ \widehat{\theta}_0\stackrel{P}{\rightarrow} \widetilde{\theta}_0$. Furthermore, $$\begin{aligned}
 \label{eq.power1}
\sqrt{b}T\, \mathcal U_T -b^{-1/2}\widetilde{\mu}_0 &= 
\sqrt{b}T \int_{-\pi}^{\pi}\| E\widehat{\mathcal F}_{X,\lambda} - E\widehat{\mathcal F}_{Y,\lambda} \|^2_{HS}d\lambda \nonumber \\
& \ \ \ \ + \sqrt{b}T \int_{-\pi}^{\pi}\| \big(\widehat{\mathcal F}_{X,\lambda}- E\widehat{\mathcal F}_{X,\lambda}\big) -\big(\widehat{\mathcal F}_{Y,\lambda} - E\widehat{\mathcal F}_{Y,\lambda}\big) \|^2_{HS}d\lambda -b^{-1/2}\widetilde{\mu}_0 \nonumber  \\
& \ \ \ \  + 2 \sqrt{b}T \int_{-\pi}^\pi \langle  \big(\widehat{\mathcal F}_{X,\lambda}- E\widehat{\mathcal F}_{X,\lambda}\big) -\big(\widehat{\mathcal F}_{Y,\lambda} - E\widehat{\mathcal F}_{Y,\lambda}\big), E\widehat{\mathcal F}_{X,\lambda} - E\widehat{\mathcal F}_{Y,\lambda}  \rangle_{HS} d\lambda.\end{aligned}$$ The assertion of the proposition follows since $ \sqrt{b}T \int_{-\pi}^{\pi}\| E\widehat{\mathcal F}_{X,\lambda} - E\widehat{\mathcal F}_{Y,\lambda} \|^2_{HS}d\lambda=\sqrt{b}T \int_{-\pi}^{\pi}\| \mathcal F_{X,\lambda} - \mathcal F_{Y,\lambda} \|^2_{HS}d\lambda + o(1)$ (see Lemma \[l.centering\]), $ \sqrt{b}T \int_{-\pi}^{\pi}\| \big(\widehat{\mathcal F}_{X,\lambda}- E\widehat{\mathcal F}_{X,\lambda}\big) -\big(\widehat{\mathcal F}_{Y,\lambda} - E\widehat{\mathcal F}_{Y,\lambda}\big) \|^2_{HS}d\lambda -b^{-1/2}\widetilde{\mu}_0 =O_P(1)$, while the last term in (\[eq.power1\]) is $ O_{P}(\sqrt{T})$ since $\sqrt{b} T$ times the expression of the inner product is bounded by $$\sqrt{T}\Big\{\|\sqrt{Tb}(\widehat{\mathcal F}_{X,\lambda}- E\widehat{\mathcal F}_{X,\lambda} )\|_{HS}  + \|\sqrt{Tb}(\widehat{\mathcal F}_{Y,\lambda}- E\widehat{\mathcal F}_{Y,\lambda}) \|_{HS}\Big\}\Big\{\| E\widehat{\mathcal F}_{X,\lambda} - E\widehat{\mathcal F}_{Y,\lambda} \|_{HS}\Big\}=O_P(\sqrt{T}).$$ Notice that the last equality follows because the sequences $ \sqrt{Tb}(\widehat{\mathcal F}_{X,\lambda}- E\widehat{\mathcal F}_{X,\lambda})$ and $\sqrt{Tb}(\widehat{\mathcal F}_{Y,\lambda}- E\widehat{\mathcal F}_{Y,\lambda}) $ are bounded in probability, uniformly in $\lambda$, since as in Theorem 3.7 of Panaretos and Tavakoli (2013), both sequences converge weakly to Gaussian elements in $ L^2([0,1]^2,\C)$ with covariance kernels that can be bounded uniformly in $\lambda$.

This research is partly funded by the Volkswagen Foundation (Professorinnen für Niedersachsen des Niedersächsischen Vorab) and by a University of Cyprus Research Grant.

Beltrão, K.I. & Bloomfield, P. (1987). Determining the bandwidth of a kernel spectrum estimate. [*Journal of Time Series Analysis*]{}, Vol. [**8**]{}, 21–38.

Benko, M., Härdle, W. & Kneip, A. (2009). Common functional principal components. [*Annals of Statistics*]{}, Vol. [**37**]{}, 1–34.

Brockwell, P.J. & Davis, R. A. (1991). [*Time Series: Theory and Methods*]{}. New York: Springer.

Cerovecki, C. & Hörmann, S. (2017). On the CLT for discrete Fourier transforms of functional time series. [*Journal of Multivariate Analysis*]{}, Vol. [**154**]{}, 281–295.

de Jong, P. (1987). A central limit theorem for generalized quadratic forms. [*Probability Theory and Related Fields*]{}, Vol. [**75**]{}, 261–277.

Dette, H. & Paparoditis, E. (2009). Bootstrapping frequency domain tests in multivariate time series with an application to comparing spectral densities. [*Journal of the Royal Statistical Society*]{}, Series B, Vol. [**71**]{}, 831–857.

Eichler, M. (2008). Testing nonparametric and semiparametric hypotheses in vector stationary processes. [*Journal of Multivariate Analysis*]{}, Vol. [**99**]{}, 968–1009.

Franke, J. & Härdle, W. (1992). On bootstrapping kernel spectral estimates. [*Annals of Statistics*]{}, Vol. [**20**]{}, 121–145.

Fremdt, S., Steinebach, J.G., Horváth, L. & Kokoszka, P. (2012). Testing the equality of covariance operators in functional samples. [*Scandinavian Journal of Statistics*]{}, Vol. [**40**]{}, 38–152.

Härdle, W. & Mammen, E. (1993). Comparing nonparametric versus parametric regression fits. [*The Annals of Statistics*]{}, Vol. [**21**]{}, 1926–1947. Hörmann, S., Kidziński, L. & Hallin, M. (2015). Dynamic functional principal components. [*Journal of the Royal Statistical Society*]{}, Series B, Vol. [**77**]{}, 319–348.

Horváth, L. & Kokoszka, P. (2012). [*Inference for Functional Data with Applications*]{}. New York: Springer-Verlag.

Horváth, L., Kokoszka, P. & Reed, R. (2013). Estimation of the mean of functional time series and a two-sample problem. [*Journal of the Royal Statistical Society*]{}, Series B, Vol. [**75**]{}, 103-122.

Hurvich, C.M. (1985). Data riven choice of a spectrum estimate: extending the applicability of cross-validation methods. [*Journal of the American Statistical Association*]{}, Vol. [**80**]{}, 933–940.

Panaretos, V.M., Kraus, D. & Maddocks, J.H. (2010). Second-order comparison of Gaussian random functions and the geometry of DNA minicircles. [*Journal of the American Statistical Association*]{}, Vol. [**105**]{}, 670–682.

Panaretos, V.M. & Tavakoli, S. (2013). Fourier analysis of stationary time series in function space. [*The Annals of Statistics*]{}, Vol. [**41**]{}, 568–603.

Paparoditis, E. (2000). Spectral density based goodness-of-fit tests for time series models. [*Scandinavian Journal of Statistics*]{}, Vol. [**27**]{}, 143–176.

Paparoditis, E. & Sapatinas, T. (2016). Bootstrap-based testing of equality of mean functions or equality of covariance operators for functional data. [*Biometrika*]{}, Vol. [**103**]{}, 727–733,

Pigoli, D., Aston, J.A.D., Dryden, I.L. & Secchi, P. (2014). Distances and inference for covariance operators. [*Biometrika*]{}, Vol. [**101**]{}, 409–422.

Robinson, P.M. (1991). Automatic frequency domain inference on semiparametric and nonparametric models. [*Econometrica*]{}, Vol. [**59**]{}, 1329–1363

Taniguchi, M. & Kakizawa, Y. (2000). [*Asymptotic Theory of Statistical Inference for Time Series.*]{} New York: Springer.

Tavakoli, S. & Panaretos, V.M. (2016). Detecting and localizing differences in functional time series dynamics: a case study in molecular biophysics. [*Journal of the American Statistical Association*]{}, Vol. [**111**]{}, 1020–1035.

van Delft, A. & Eichler, M. (2018). Locally stationary functional time series. [*Electronic Journal of Statistics*]{}, Vol. [**12**]{}, 107–170.

Wu, W.B. & Zaffaroni, P. (2015). Uniform convergence of multivariate spectral density estimates. [*Preprint*]{} (arXiv:1505.03659).

Zhang, X. & Shao, X. (2015). Two sample inference for the second-order property of temporally dependent functional data. [*Bernoulli*]{}, Vol. [**21**]{}, 909–929.

Zhang, C., Peng, H. & Zhang, J.-T. (2010). Two samples tests for functional data. [*Communications in Statistics - Theory and Methods*]{}, Vol. [**39**]{}, 559–578.
