John Aston
email: j.aston@statslab.cam.uk

Dominique Dehay
email: dominique.dehay@univ-rennes2.fr

Anna Dudek
email: aedudek@agh.edu.pl

Jean-Marc Freyermuth
email: jean-marc.freyermuth@univ-amu.fr

Denes Szucs

Lincoln Colling
email: l.colling@sussex.ac.uk

Spectrum inference for replicated spatial locally time-harmonizable time series

Keywords: Primary 62G05, 62G20; secondary 62M15 Harmonizable spatiotemporal processes, Nonparametric spectral analysis, Circular Block Bootstrap, Functional connectivity, ElectroEncephaloGraphy

In this paper, we develop tools for statistical inference on replicated realizations of spatiotemporal processes that are locally time-harmonizable. Our method estimates both the rescaled spatial timevarying Loève-spectrum and the spatial time-varying dual-frequency coherence function under realistic modeling assumptions. We construct confidence intervals for these parameters of interest using the Circular Block Bootstrap method and prove its consistency. We illustrate the application of our methodology on a dataset arising from an experiment in neuropsychology. From EEG recordings, our method allows studying the dynamic functional connectivity within the brain associated to visual working memory performance.

Introduction

The paper is concerned with a class of spatiotemporal processes that are locally time-harmonizable, that is, they possess a local two-dimensional spectrum. In order to introduce such a class of processes, we need first to recall some basic facts concerning second order harmonizable processes that are due to [START_REF] Loève | Fonctions aléatoire du second ordre[END_REF]. A centered P -variate discrete-time process {X} = {X t , t ∈ Z} with the finite second order moments is called harmonizable if it admits a Cramér's representation of the following form:

X t = π -π e itω dZ (ω) , (1.1)
where X t = (X 1,t , . . . , X P,t ) ′ and the spectral process {Z(ω)} = {Z(ω) = (Z 1 (ω), . . . , Z P (ω)) ′ , ω ∈ (-π, π]} is a zero-mean stochastic process. Here and hereafter the symbol (•) ′ denotes the transpose of a vector.

The simplest example of processes admitting the representation (1.1) are stationary sequences. In that case, the process {Z(ω)} has orthogonal and cross-orthogonal increments (see e.g., [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]). However, for the class of harmonizable processes, {Z(ω)} has correlated increments. They form a broad class of processes that includes many nonstationary ones, such as periodically correlated time series (see e.g. [START_REF] Hurd | Periodically Correlated Random Sequences: Spectral Theory and Practice[END_REF]). A very important feature of harmonizable processes is that their covariance is a Fourier transform of a finite measure.

Harmonizable processes are particularly useful in modeling real-world data when the main interest is frequency domain analysis. They are widely applied, for instance in signal theory, communications and mechanics (see e.g. [START_REF] Napolitano | Cyclostationarity: New trends and applications[END_REF][START_REF] Gardner | Cyclostationarity : half a century of research[END_REF][START_REF] Setoodeh | Fundamentals of Cognitive Radio[END_REF][START_REF] Serpedin | Bibliography on cyclostationarity[END_REF]). Our approach was initially motivated by the analysis of ElectroEncephaloGraphy (EEG) data. However, it can be applied to other types of problems with possibly some minor modifications.

In recent years, numerous studies of brain signals have explored networks of functional connections to reveal subtle mechanisms of brain activity. Essentially, this involves measuring the relationships in the activity of different brain regions. The analysis is often performed using coherence, which is a frequency domain equivalent of correlation. It takes value in [0, 1], a value close to 1 indicates a strong synchronization. More specifically, for a P -variate harmonizable process {X}, its Loève spectrum f = f pq 1≤p,q≤P is a P × P -matrix defined as follows Cov(dZ(ω 1 ), dZ(ω 2 )) = f (ω 1 , ω 2 )dω 1 dω 2 .

(1.

2)

The dual-frequency coherence between a pair of processes ({X p,t }, {X q,t }) and a pair of frequencies (ω 1 , ω 2 ) is given by:

ρ pq (ω 1 , ω 2 ) := |Cov (dZ p (ω 1 ) , dZ q (ω 2 ))| 2
Var (dZ p (ω 1 )) Var (dZ q (ω 2 ))

= |f pq (ω 1 , ω 2 )| 2 f pp (ω 1 , ω 1 ) f qq (ω 2 , ω 2 ) . (1.3)
The dual-frequency coherence (1.3) allows capturing dependencies at two different frequencies. [START_REF] Soedjak | Consistent estimation of the bispectral density function of a harmonizable process[END_REF] developed inference tools for the Loève spectrum for such a model, but the Loève spectrum of the harmonizable process is constant in time and the model does not consider any spatial localization. Therefore, it cannot sufficiently capture the complexity of the brain mechanisms. Consequently, these results needed to be extended accordingly. This was achieved thanks to the recent and important contribution of [START_REF] Gorrostieta | Time-dependent dual-frequency coherence in multivariate non-stationary time series[END_REF]. The authors follow the approach of [START_REF] Dahlhaus | Locally stationary processes[END_REF] to introduce multivariate locally-harmonizable processes. They describe a windowed Fourier based estimation procedure for the time-varying dual-frequency coherence. They derive exact confidence intervals for testing if the coherence differs from zero under i.i.d. Gaussian assumptions, and also obtain asymptotic confidence intervals.

In this paper, we extend the existing results in several ways. First, we introduce new inference tools that take into account both time and space (i.e. spatial location). We define the rescaled spatiotemporal local Loève spectrum and the spatiotemporal coherence. In other words, we measure the time-evolving squared correlation coefficient at different frequencies between any pairs of spatial locations. Our approach uses spatial correlations to improve the estimation of these quantities by exploiting spatial location information in the spirit of the [START_REF] Ombao | Spatio-spectral analysis of brain signals[END_REF] method. Second, we consider more realistic modeling assumptions. Third, in order to construct confidence intervals for the spatiotemporal coherence, we adapt the Circular Block Bootstrap (CBB) method and show its consistency.

The paper is organized as follows. In Section 2 we introduce a spatial locally time-harmonizable process model along with an appropriate estimation procedure under realistic model assumptions. In Section 3 we discuss asymptotic properties of our estimators. Moreover, we show consistency of the CBB approach. Finally, in Section 4 we illustrate the application of our method on a real data set. All proofs and additional information on the real data can be found in Section 6.

Rescaled spatiotemporal spectrum estimation

In this section, we generalize some of the ideas presented by [START_REF] Ombao | Spatio-spectral analysis of brain signals[END_REF] and [START_REF] Gorrostieta | Time-dependent dual-frequency coherence in multivariate non-stationary time series[END_REF]. For the sake of clarity, we start by introducing the notion of spatial time-harmonizable process and the corresponding Loève spectrum. Next, we introduce the spatial time-varying local Loève spectrum for a general spatial process. Then, we describe our modeling assumptions, in particular the spatiotemporal rescaling. They ensure notably that the quantities of interest lie on a bounded spatiotemporal domain and satisfy some smoothness conditions. We construct a rescaled spectrum estimator that is based on replicated observations of the process, and give its asymptotic properties. Finally, we adapt the CBB method to construct bootstrap confidence intervals and we prove the bootstrap consistency.

Spatiotemporal Loève spectrum

Let X S = X S t , t ∈ Z := X S s,t , t ∈ Z, s ∈ {1, . . . , S 1 } × {1, . . . , S 2 } , S := (S 1 , S 2 ) ∈ N * 2 , be a family of spatial time-harmonizable processes, i.e.,

X S t = π -π e iωt dZ S (ω) , such that Cov dZ S (ω 1 ), dZ S (ω 2 ) = f S (ω 1 , ω 2 ) dω 1 dω 2 , where f S := f S s 1 ,s 2 s 1 ,s 2 ∈{1,...,S1}×{1,...,S2} is the Loève spectrum. Then C S (t 1 , t 2 ) := Cov X S t1 , X S t2 = π -π π -π f S (ω 1 , ω 2 ) e i(ω1t1-ω2t2) dω 1 dω 2 .
Here C S (t 1 , t 2 ) and f S (ω 1 , ω 2 ) are S 1 × S 2 × S 1 × S 2 -matrices, N * := {1, 2, . . . }. Notice that the process is defined for each spatial locations s 1 , s 2 and each time points (t 1 , t 2 ). A sufficient condition for time-harmonizability and the existence of a two-dimensional spectral density for the discrete-time spatial second order random processes X S is given by the following condition

(t1,t2)∈Z 2 C S (t 1 , t 2 ) < ∞,
where |•| is a matrix norm. Then the Loève spectrum is a continuous function and it coincides with

f S (ω 1 , ω 2 ) = 1 4π 2 (t1,t2)∈Z 2 C S (t 1 , t 2 ) e -i(ω1t1-ω2t2) .
Remark that the above definition does not include the stationary case as the Loève spectrum is twodimensional while the spectrum of a stationary process is one-dimensional.

Localized Loève spectrum

For the purpose of our application, the notion of harmonizable processes is not sufficient. Therefore, in this section, we generalize the previous considerations by introducing the notion of spectrum for a spatial second order process {X S } that is not necessarily time-harmonizable. We also introduce its estimator. We define the (spatiotemporal) localized Loève spectrum of the process {X S } as

( f S t1,t2 (ω 1 , ω 2 ) := 1 4π 2 t1+N -1 k1=t1-N t2+N -1 k2=t2-N C S (k 1 , k 2 )e -i(ω1k1-ω2k2) (2.1)
using a local rectangular time window centered at (t 1 , t 2 ) with size 2N . For any t 1 and t 2 we obtain that

C S (t 1 , t 2 ) = π 2 N 2 N -1 j1=-N N -1 j2=-N ( f S t1,t2 ω N j1 , ω N j2 e i ω N j 1 t1-ω N j 2 t2 , (2.2)
where w N j := jπ N , j = -N, . . . , N -1 are the Fourier frequencies of the local rectangular time window. When {X S t } is a family of spatial time-harmonizable processes with spectrum f S (ω 1 , ω 2 ) one can easily verify that

( f S t1,t2 (ω 1 , ω 2 ) = 1 4π 2 π -π π -π D N (ω ′ 1 )D N (ω ′ 2 )e i ω ′ 1 t1-ω ′ 2 t2 f S ω 1 + ω ′ 1 , ω 2 + ω ′ 2 dω ′ 1 dω ′ 2 ,
where D N (0) = 2N and D N (ω) = 2i sin(ωN )

e iω -1
otherwise. Furthermore, if

(t1,t2)∈Z 2 C S (t 1 , t 2 ) < ∞,
for any k 1 and k 2 , then lim

N →∞ ( f S t1,t2 (ω 1 , ω 2 ) = f S (ω 1 , ω 2 ).

The observations

In the following, we consider replicates {X S,r }, r ∈ N * , of a spatial zero-mean second order process {X S }. This means that the processes {X S,r } have the same distribution as {X S }. Here the process is not necessarily time-harmonizable. From now on, we assume that the replicates are dependent, more precisely, that the family of processes {X S,r }, r ∈ N * , is nonstationary with respect to t and stationary with respect to r. Consequently, we denote

C S (t 1 , t 2 ) := Cov X S,1 t1 , X S,1 t2 = Cov X S,r t1 , X S,r t2 (2.3)
for any positive integer r, and

( f S t1,t2 (ω 1 , ω 2 ) := 1 4π 2 t1+N -1 k1=t1-N t2+N -1 k2=t2-N C S (k 1 , k 2 ) e -i(ω1k1-ω2k2) . (2.4) Then C S (t 1 , t 2 ) = π 2 N 2 N -1 j1=-N N -1 j2=-N ( f S t1,t2 ω N j1 , ω N j2 e i ω N j 1 t1-ω N j 2 t2 ,
where w N j := jπ N , j = -N, . . . , N -1. The r-th replicate is observed at time instants 0, . . . , T -1 and at S 1 × S 2 different spatial locations. Hence

{X S,r } = {X S,r t , t = 0, . . . , T -1} = X S,r s,t , t = 0, . . . , T -1, s ∈ {1, . . . , S 1 } × {1, . . . , S 2 } ,
where S = (S 1 , S 2 ) ∈ N * 2 . For the sake of simplicity we set X S,r t = 0 S1×S2 (the null S 1 × S 2 -matrix) for t / ∈ {0, . . . , T -1}.

In the following, we study the asymptotic behavior of the localized Loève spectrum ( f S t1,t2 (ω 1 , ω 2 ). For that purpose we introduce the rescaled spatiotemporal spectrum and we construct its estimator. All asymptotic results are obtained as S 1 , S 2 , T, R go to ∞. The time window size 2N can be fixed or going to ∞.

Assumptions

To obtain the asymptotic results we assume the following conditions.

(L) Rescaling conditions. There exists a function f : [0, 1] 6 × (-π, π] 2 → C and positive constants L and Q such that

f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) -f u 3 ,u 4 ,τ3,τ4 (ω 1 , ω 2 ) ≤ L (∥u 1 -u 3 ∥ + ∥u 2 -u 4 ∥ + |τ 1 -τ 3 | + |τ 2 -τ 4 |) (2.5) for any u 1 , u 2 , u 3 , u 4 , ∈ [0, 1] 2 , τ 1 , τ 2 ∈ [0, 1] and ω 1 , ω 2 ∈ (-π, π] and ( f S s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) -f s1 ,s 2 , ẗ1, ẗ2 (ω 1 , ω 2 ) ≤ Q 1 S 1 + 1 S 2 + 1 T , (2.6)
where

s i := (s i,1 , s i,2 ), si := (s i,1 /S 1 , s i,2 /S 2 ), ẗi := t i /T, i = 1, 2 for N ≤ t 1 , t 2 ≤ T -N . Inequality (2.6
) is assumed to be true for all S 1 , S 2 and T large enough, and for n fixed or sufficiently large, as the case may be.

Hereafter, f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) is called the rescaled Loève spectrum (SR) The replications {X S,r }, r ∈ N * have the same distribution and are stationary with respect to r. (MR) Mixing property for the replicates: The family {X S,r }, r ∈ N * , S ∈ N * 2 , is α-mixing with respect to r and such that one of the following two conditions holds: (i) sup t,S X S,1 t

< C almost surely for some finite C > 0 and where F r (S) := σ X S,q s,t : q ≤ r, t ∈ Z and all locations s and F r+κ (S) := σ X S,q s,t : q ≥ r + κ, t ∈ Z and all locations s .

κ α X (κ) < ∞, (ii) sup t,S E X S,1 t 4+δ < ∞ and κ α X (κ) δ/(4+δ) < ∞
In order to state the asymptotic covariance of our estimator we consider an additional rescaling assumption, which is a generalization of the condition (L) to the four-dimensional spectrum. Denote t := (t 1 , t 2 , t 3 , t 4 ) ∈ Z 4 , τ := (τ 1 , τ 2 , τ 3 , τ 4 ) ∈ [0, 1] 4 and ω := (ω 1 , ω 2 , ω 3 , ω 4 ) ∈ (-π, π] 4 . Moreover, for s j ∈ N * 2 and u j ∈ [0, 1] 2 , j = 1, . . . , 4, let s := (s 1 , s 2 , s 3 , s 4 ) ∈ N * 8 and u := (u 1 , u 2 , u 3 , u 4 ) ∈ [0, 1] 8 . Under the stationary condition (SR) the covariance C S,κ s,t := Cov X κ+r s 1 ,t1 X κ+r s 2 ,t2 , X r s 3 ,t3 X r s 4 ,t4

does not depend on r ≥ max{0, -κ}, for any κ ∈ Z. Then define

( f S,κ s,t (ω) := 1 16π 4 t1+N -1 k1=t1-N t2+N -1 k2=t2-N t3+N -1 k3=t3-N t4+N -1 k4=t4-N C S,κ s,k e -i(ω1k1-ω2k2-ω3k3+ω4k4) .
The rescaling assumption is as follows.

(LR) Rescaling condition for the replicates. There exist functions f κ : [0, 1] 12 × (-π, π] 4 → C, κ ∈ N * , and some positive constants L and Q such that for each

u i ∈ [0, 1] 8 , τ i ∈ [0, 1] 4 , ω ∈ (-π, π] 4 , i = 1, 2 and each κ ∈ N * , f κ u1,τ1 (ω) -f κ u2,τ2 (ω) ≤ L 4 j=1 ∥u j,1 -u j,2 ∥ + |τ j,1 -τ j,2 | (2.7) and ( f S,κ s,t (ω) -f κ s, ẗ(ω) ≤ Q 1 S 1 + 1 S 2 + 1 T , (2.8)
where s j = (s j,1 , s j,2 ), sj = (s j,1 /S 1 , s j,2 /S 2 ), ẗj = t j /T, j = 1, 2, 3, 4 and for

N ≤ t 1 , t 2 , t 3 , t 4 ≤ T -N . Furthermore, assume that κ∈Z f κ u,τ (ω M l ) < ∞, (2.9)
for the Fourier frequencies ω M l = (ω M l1 , . . . , ω M l4 ) and ω M li = liπ M , l i = -M, . . . , M -1}, i = 1, . . . , 4, where the integer M > 0 is fixed and N = nM . Inequality (2.8) is assumed to be true for all S 1 , S 2 and T sufficiently large, and for n fixed or sufficiently large, as the case may be. 1. Under conditions (L) and (LR) the functions f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) and f κ u,τ (ω) are L-Lipschitz-continuous in space and time components uniformly with respect to the frequencies ω 1 , . . . , ω 4 and the shift κ between replicates. 2. Identifiability. In condition (L), relation (2.6) is assumed to be true for all S 1 , S 2 and T sufficiently large.

Hence, if f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) exists then it is unique. Similarly, under condition (LR) the function f κ u,τ (ω) is unique. 3. When we assume that N → ∞, then the rescaled Loève spectrum f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) does not depend on N . Of course, if N is assumed to be fixed then f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) may depend on N . 4. Example for condition (L). Let {X S } be a spatial time-harmonizable process with the Loève spectrum of the form

f S s 1 ,s 2 (ω 1 , ω 2 ) = A S (s 1 , s 2 ) ϕ(ω 1 , ω 2 )
, where the function ϕ(ω 1 , ω 2 ) is bounded, say |ϕ(ω 1 , ω 2 )| ≤ c, c > 0, and

A S (s 1 , s 2 ) -A(s 1 , s2 ) ≤ Q c 1 S 1 + 1 S 2 for some (L/c)-Lipschitz-continuous function A : [0, 1] 4 → C. Then assumption (L) is fulfilled with f s1 ,s 2 , ẗ1, ẗ2 (ω 1 , ω 2 ) = A(s 1 , s2 ) ϕ(ω 1 , ω 2 )
. (See also [START_REF] Ombao | Spatio-spectral analysis of brain signals[END_REF]) 5. The α-mixing function α X is a weak dependence measure. Hence, replicated processes {X S,r1 } and {X S,r2 } that are close to each other, i.e. such that the distance κ := |r 1 -r 2 | between replications is small, can be dependent, while when κ is large, they are almost independent. The replicates are Mdependent, M ≥ 1, if and only if α X (κ) = 0 for any κ ≥ M . This generalizes the modeling assumptions in [START_REF] Gorrostieta | Time-dependent dual-frequency coherence in multivariate non-stationary time series[END_REF], where the replicates are assumed to be independent, that is α X (κ) = 0 for any κ ̸ = 0. For properties and examples of other dependence measures, we refer the reader to [START_REF] Doukhan | Mixing: properties and examples[END_REF]. 6. Gaussian framework. For pedagogical purposes, we present the results for a Gaussian process in Section 6.1 in the appendix. In this case, we do not need the mixing condition, and we replace the condition (LR) by (LGR). Then we give an expression for the four-dimensional rescaled spectrum f κ u,τ ω in terms of the two-dimensional rescaled spectrum. See relation (6.3).

In the following, we provide the results in two cases: N fixed and N → ∞, that is n fixed and n → ∞ for M fixed with N = nM . The integer M being defined below according to the frequency resolution. The case N → ∞ denotes that we consider infinitely many time points around each instant t.

Estimator of the rescaled Loève spectrum

In this section we introduce an estimation procedure for the rescaled Loève spectrum f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ). For that purpose, we first define two kernel functions that we use for rescaling in space and time. To simplify the presentation, let us consider two non-negative functions w, W : R → [0, ∞) and two positive numbers h and ℏ. We define

w u (s) := 1 S 1 S 2 h 2 w u 1 -s 1 /S 1 h w u 2 -s 2 /S 2 h , and
W τ (t) := 1 T ℏ W τ -t/T ℏ where u = (u 1 , u 2 ) ∈ [0, 1] 2 , s = (s 1 , s 2 ) ∈ {1, . . . ,
, k ′ ∈ N * , -1 = υ 1 < • • • < υ k = 1 and -1 = τ 1 < • • • < τ k ′ = 1 such that w(•) and W (•) are Lipschitz-continuous on each interval (υ j , υ j+1 ), 1 ≤ j ≤ k -1 and (τ j ′ , τ j ′ +1 ), 1 ≤ j ′ ≤ k ′ -1, respectively.
Note that under the condition (KS) the kernel functions w(•) and W (•) are bounded. It holds for instance for rectangular and triangular kernels. Now we define the dual-frequency periodogram of the r-th replicate for the spatial locations s 1 , s 2 and the instants t 1 , t 2 at frequencies ω 1 , ω 2 and over a time window of size 2N as

I r s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) := 1 4π 2 d r s 1 ,t1 (ω 1 ) d r s 2 ,t2 (ω 2 ), where d r s,t (ω) := t+N -1 k=t-N X S,r s,k e -iωk = N -1 k=-N X S,r s,k+t e -iω(k+t)
is the discrete Fourier transform of the r-th replicate for the spatial location s around the instant t. Recall that we set X r s,k = 0 for k / ∈ {0, . . . , T -1}. Then the estimator of the local Loève spectrum is defined as the average of the dual-frequency periodograms of replicates i.e.,

fs 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) := 1 R R r=1 I r s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) .
Finally, the estimator of the rescaled Loève spectrum

f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) is given by f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) := t1 t2 s 1 s2 W τ1 (t 1 )W τ2 (t 2 )w u 1 (s 1 )w u 2 (s 2 ) fs 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) . (2.10)
The rescaled coherence is defined as

ρ u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) := f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) 2 f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) , (2.11)
and its estimator is given by

ρ u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) := f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) 2 f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) . (2.12)
Due to the limitation of the frequency resolution capacity in the real life experiment, in the sequel we consider the convergence of the estimator

f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) for a finite number of Fourier frequencies ω M l := lπ M , -M ≤ l ≤ M -1,
where M > 1 is some fixed integer. Furthermore, in order to ensure the identifiability of the frequencies, we take the window size 2N equal to an integer multiple of M : N = nM . This choice of modeling allows us to derive a more accessible asymptotic theory presented in Section 3. Moreover, it is motivated by our real data application for which we typically consider a finite number of frequency bands of interest. In particular, we consider the sets of Fourier frequencies Ω

i := ω M j : L i ≤ j ≤ L i + l i -1 , for some l i ≥ 1, i = 1, 2.
Then the estimator is computed as an average over the frequencies

f u 1 ,u 2 ,τ1,τ2 (Ω 1 , Ω 2 ) := 1 l 1 l 2 L1+l1-1 j1=L1 L2+l2-1 j2=L2 f u 1 ,u 2 ,τ1,τ2 ω M j1 , ω M j2 .
(2.13)

Main results

Below we state some asymptotic properties of our estimation procedure like convergence in quadratic mean and asymptotic normality. All the proofs are deferred to Section 6.2 in the appendix. From now on, any complex number z is treated as a vector of its real and imaginary parts, i.e., z = (ℜz, ℑz) ′ .

Theorem 3.1. (Convergence in quadratic mean) Let M ≥ 1, u 1 , u 2 ∈ (0, 1) 2 and τ 1 , τ 2 ∈ (0, 1) be fixed. Assume that the assumptions (L) and (MR) hold. Then

lim R→∞ f u 1 ,u 2 ,τ1,τ2 ω M l1 , ω M l2 = f u 1 ,u 2 ,τ1,τ2 ω M l1 , ω M l2 in quadratic mean, for the Fourier frequencies ω M li = liπ M , l i = -M, . . . , M -1, i = 1, 2, provided that N = nM as well as h 2 S 1 , h 2 S 2 , ℏ 2 T → ∞ and n 2 (h + ℏ), n 4 R -1 → 0 as T, S 1 , S 2 , R → ∞, h, ℏ → 0 independently of the behavior of n ≥ 1.
Below we express the asymptotic covariance matrix of the estimator

f u 1 ,u 2 ,τ1,τ2 . Proposition 3.1. Let M ≥ 1, u i ∈ (0, 1) 2 , τ i ∈ (0, 1
) and the Fourier frequencies ω M li be fixed, i = 1, 2, 3, 4. Assume that the assumptions (SR) and (LR) hold. Then

lim R→∞ R Cov f u 1 ,u 2 ,τ1,τ2 ω M l1 , ω M l2 , f u 3 ,u 4 ,τ3,τ4 ω M l3 , ω M l4 = κ∈Z f κ u,τ ω M l , provided that N = nM as well as h 2 S 1 , h 2 S 2 , ℏ 2 T → ∞ and n 4 (h + ℏ) → 0 as T, S 1 , S 2 , R → ∞, h, ℏ → 0 independently of the behavior of n ≥ 1.
Before we formulate the multivariate central limit theorem we introduce some additional notation. Let

f := f u 1,1 ,u 2,1 ,τ1,1,τ2,1 (ω 1,1 , ω 2,1 ) ′ , . . . , f u 1,k ,u 2,k ,τ 1,k ,τ 2,k (ω 1,k , ω 2,k ) ′ ′ , and
f := f u 1,1 ,u 2,1 ,τ1,1,τ2,1 (ω 1,1 , ω 2,1 ) ′ , . . . , f u 1,k ,u 2,k ,τ 1,k ,τ 2,k (ω 1,k , ω 2,k ) ′ ′ ,
where k is some positive integer, u i,j ∈ (0, 1) 2 , τ i,j ∈ (0, 1), ω i,j = ω M li,j = li,j π M , i = 1, 2 and j = 1, . . . , k.

Now we state the asymptotic normality of the estimator.

Theorem 3.2. Assume that the assumptions (L), (SR), (MR) and (LR) hold. Then

lim R→∞ L √ R f -f = N 2k (0, Σ 2k ) , provided that (i) either N = nM is a constant, T, S 1 , S 2 , R → ∞, h, ℏ → 0 with Rh -4 (S -2 1 + S -2 2 ), Rℏ -4 T -2 , R(h 2 + ℏ 2 ) → 0; (ii) or N = nM → ∞, T, S 1 , S 2 , R → ∞, h, ℏ → 0 with Rh -4 (S -2 1 + S -2 2 ), Rℏ -4 T -2 , Rn 4 (h 2 + ℏ 2 ) → 0, and the additional condition n T ℏ T -1 t=0 X S,1 s i ,t ≤ C (3.1)
almost surely, or

n T ℏ T -1 t=0 E X S,1 s i ,t 4+δ 1/(4+δ) ≤ C (3.2)
for some finite C > 0 which does not depend on the locations. The elements of the covariance (2k × 2k)-matrix Σ 2k can be calculated from Proposition 3.1.

Remark 3.1. When N → ∞, conditions (3.1) and (3.2) can be replaced by more subtle assumptions. For the sake of clarity, this technical remark is detailed in the annexe. See conditions (ii) in Proposition6.1 in the annexe and to the subsequent remarks.

Theorem 3.2 is crucial to study the behavior of f u 1 ,u 2 ,τ1,τ2 (Ω 1 , Ω 2 ) given by the equation (2.13).

Corollary 3.1. Under conditions of Theorem 3.2, the estimator ρ of the rescaled spatiotemporal coherence ρ, defined respectively by (2.12) is asymptotically normal i.e.,

lim R→∞ L √ R ρu 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) -ρ u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) = N (0, γ 2 ), (3.3)
where the Fourier frequencies

ω i = ω M li with l i ∈ {M, . . . , M -1}, i = 1, 2 and provided that f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 )× f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) ̸ = 0.
Here γ 2 = (∇ (ω1,ω2)

u 1 ,u 2 ,τ1,τ2 )Σ 6 (∇ (ω1,ω2) u 1 ,u 2 ,τ1,τ2
) ′ , where ∇ denotes the gradient operator. The covariance 6 × 6 matrix

Σ 6 is given in Theorem 3.2 for k = 3, τ 1,1 = τ 2,1 = τ 1,3 = τ 1 , τ 1,2 = τ 2,2 = τ 2,3 = τ 2 , u 1,1 = u 2,1 = u 1,3 = u 1 , u 1,2 = u 2,2 = u 2,3 = u 2 , ω 1,1 = ω 2,1 = ω 1,3 = ω 1 , and ω 1,2 = ω 2,2 = ω 2,3 = ω 2 . Moreover, ∇ (ω1,ω2) u 1 ,u 2 ,τ1,τ2 = -f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) 2 f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) 2 f u 2 ,u 2 ,τ2,τ2, (ω 2 , ω 2 ) , 0 , -f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) 2 f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) 2 , 0 , 2ℜf u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) , 2ℑf u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) ′ .

Bootstrap approach

Using Corollary 3.1 one may construct confidence interval for the spatiotemporal dual-frequency coherence ρ u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ). However, since the asymptotic variance γ 2 depends on unknown parameters, it is in practice very difficult to estimate. Thus, we present below a bootstrap approach that allows to obtain consistent confidence intervals for ρ τ1,τ2,u 1 ,u 2 (ω 1 , ω 2 ).

Let us recall that we have R replicates {X (r) } = {X S,r s,t , t ∈ Z, s ∈ {1, . . . , S 1 } × {1, . . . , S 2 }}, r = 1, . . . , R. The process X S,r s,t is stationary in r and nonstationary in t. We will bootstrap our observations in replicates not in time. For that purpose we use the CBB (see [START_REF] Politis | A circular block-resampling procedure for stationary data[END_REF]). The CBB is a modification of the Moving Block Bootstrap method [START_REF] Künsch | The jackknife and the bootstrap for general stationary observations[END_REF][START_REF] Liu | Moving block jackknife and bootstrap capture weak dependence[END_REF], which allows to reduce bias of the bootstrap estimator. Below we present how to adapt the CBB algorithm to our problem. Let B i , i = 1, . . . , R be the block of replicates from our sample X (1) , . . . , X (R) , that starts with replicate X (i) and has the length b ∈ N, i.e.

B i := X (i) , . . . , X (i+b-1) .

If i + b -1 > R then the missing part of the block is taken from the beginning of the sample and we get

B i = X (i) , . . . , X (R) , X (1) , . . . , X (b-R+i-1) for i = R -b + 2, . . . , R.

CBB algorithm

1. Choose a block size b < R. Then our sample X (1) , . . . , X (R) can be divided into l blocks of length b and the remaining part is of length r, i.e. R = lb + r, R = 0, . . . , b -1. 2. From the set {B 1 , . . . , B R } choose randomly with replacement l + 1 blocks. 3. Join the selected l + 1 blocks (B * 1 , . . . , B * l+1 ) and take the first R observations to get the bootstrap sample X * (1) , . . . , X * (R) of the same length as the original one.

We apply the CBB to get bootstrap estimators of f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) and ρ u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) and finally to be able to construct confidence intervals for these characteristics. We use the bootstrap algorithm described above. The bootstrap version of f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) is given by

f * u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) := t1 t2 s 1 s2 W τ1 (t 1 )W τ2 (t 2 )w u 1 (s 1 )w u 2 (s 2 ) f * s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) , (3.4)
where

f * s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) := 1 R R r=1 I * ,r s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) , I * ,r s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) := 1 4π 2 d * ,r s 1 ,t1 (ω 1 ) d * ,r s 2 ,t2 (ω 2 ),
and

d * ,r s,t (ω) := t+N -1 k=t-N X * ,r s,k e -iωk .
Below we state the consistency of our bootstrap approach for the spatial time-varying dual-frequency coherence function. The bootstrap estimator of the spatial coherence is defined as

ρ * u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) := f * u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) 2 f * u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) f * u 2 ,u 2 τ2,τ2 (ω 2 , ω 2 )
. Theorem 3.3. Under conditions of Theorem 3.2 and assuming that b

-1 + R -1 b = o(1) the CBB is consistent i.e., sup x∈R P *    √ R   ρ * u 1 ,u 2 ,τ 1 ,τ 2 (ω1, ω2) - E * ℜ f * u 1 ,u 2 ,τ 1 ,τ 2 (ω1, ω2) 2 + E * ℑ f * u 1 ,u 2 ,τ 1 ,τ 2 (ω1, ω2) 2 E * f * u 1 ,u 1 ,τ 1 ,τ 1 (ω1, ω1) E * f * u 2 ,u 2 ,τ 2 ,τ 2 (ω2, ω2)    ≤ x    -P √ R ρu 1 ,u 2 ,τ 1 ,τ 2 (ω1, ω2) -ρu 1 ,u 2 ,τ 1 ,τ 2 (ω1, ω2) ≤ x p -→ 0 as R -→ ∞ (3.5)
for Fourier frequencies

ω i = ω M li with l i ∈ {M, . . . , M -1}, i = 1, 2. Centering of ρ * u 1 ,u 2 ,τ1,τ2, (ω 1 , ω 2 ) may seem surprising. One could expect to use simply E * (ρ * u 1 ,u 2 ,τ1,τ2, (ω 1 , ω 2 )
). But in fact the spatial time-varying dual-frequency coherence function is a function of the rescaled spatiotemporal Loève spectrum and therefore to show convergence (3.5), one needs first to obtain bootstrap consistency for f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ), then to generalize this result to a multidimensional case and finally to apply the delta method (see Propositions 6.2 and 6.3 in the annexe). While applying block bootstrap a natural question that appears concerns the choice of the block length. In the case of stationary sequences this problem is well investigated (see [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF]). It is well known that for the CBB the optimal block length obtained by minimization of the mean squared error of the bootstrap estimator is b = O(R 1/3 ) (see Theorem 5.4 in [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF]).

Real data application

We illustrate the application of our method on a dataset derived from an experiment in neuropsychology. It aims at improving our understanding of the brain mechanisms involved in Visual Working Memory performance. After a brief description of the scientific context and data, we demonstrate the usage of our methodology by providing a visualization of the estimated spatiotemporal dual-frequency coherence and an estimation of the dual-frequency functional connectivity networks.

Scientific context

Working Memory (WM) is an essential cognitive resource because it is strongly correlated with general cognitive abilities. Its function is to maintain access to relevant information during a brief time-span, which enables a person to perform activities such as navigation, communication, problem solving. . . Over the past 20 years there has been an explosion of more specific research on Visual Working Memory (VWM). Following [START_REF] Luck | Visual working memory capacity: from psychophysics and neurobiology to individual differences[END_REF], Visual Working Memory is an "active maintenance of visual information to serve the needs of ongoing tasks". There are key issues at stake in describing and identifying sources of VWM limitation and variability, particularly from the perspective of brain connectivity [START_REF] Fougnie | Variability in the quality of visual working memory[END_REF]. Brain connectivity describes how localized activity can be statistically dependent from one part of the brain to another. In the neuroscience community, this is referred to as functional connectivity [START_REF] Friston | Functional and effective connectivity: A review[END_REF]. In our data example, the study of these brain mechanisms is based on the analysis of EEG signals. In brief, electrical currents generated in the brain by ensembles of neurons firing in a synchronized manner propagate through the cerebral cortex to the scalp, where they are recorded by spatially localized EEG electrodes. These electrodes measure electric potentials over time, which represent the oscillations of the brain waves. Hence, the study of functional connectivity can be addressed using coherence analysis. It has already been proved useful in order to reveal interesting facts about Working Memory [START_REF] Sauseng | Fronto-parietal eeg coherence in theta and upper alpha reflect central executive functions of working memory[END_REF]. A challenging aspect is that these dynamic functional connections may involve brain waves oscillatory components of different frequencies [START_REF] Gorrostieta | Exploring dependence between brain signals in a monkey during learning[END_REF][START_REF] Pascual-Marqui | The dual frequency rv-coupling coefficient: a novel measure for quantifying cross-frequency information transactions in the brain[END_REF]. In other terms, bursts of high frequencies in some area of the brain could occur preferentially during specific phases of low frequency activity in other areas. It is worth noticing that electrical currents at the scalp surface are spatiotemporal phenomenon sampled at the specific localization of the electrodes. We showed that our method is an appropriate tool for modeling such a phenomenon because it consistently estimates the corresponding spatio-spectral characteristics. In fact, neuroscientists are interested in studying certain specific frequency bands that relate to different brain states and that can be interpreted in a meaningful way. More specifically, in the sequel, we consider the so-called theta, alpha and beta frequency bands ( [START_REF] Napolitano | Cyclostationarity: New trends and applications[END_REF][START_REF] Soedjak | Consistent estimation of the bispectral density function of a harmonizable process[END_REF]-Hz, [START_REF] Soedjak | Consistent estimation of the bispectral density function of a harmonizable process[END_REF][START_REF] Doukhan | Mixing: properties and examples[END_REF]-Hz and [START_REF] Doukhan | Mixing: properties and examples[END_REF][START_REF] Sauseng | Fronto-parietal eeg coherence in theta and upper alpha reflect central executive functions of working memory[END_REF]-Hz, respectively) denoted as Ω θ , Ω α and Ω β .

Experimental details

Our real data comes from an experiment that consists of the following consecutive steps (an illustration is provided in Section 6.3 in the appendix:

• Memory set -Memorize: the subject is placed in front of a computer screen. An arrow appears on the screen and the subject has 2 seconds to memorize its orientation and color.

-Retain: a blank screen appears for 0.3 seconds, then, for the next 0.1 seconds, multiple arrows appear to knock out the immediate memory. Finally, a blank screen appears again for 0.9 seconds.

• Memory test: using a joystick, the subject has 1.7 seconds to reproduce the orientation and the color of the arrow.

Notice that the subject answers about the color he remembers by selecting it from a color scale wrapped on a circle (see Section 6.3 in the Annexe). Henceforth, we compute the VWM errors for both orientation and color as angles between the truth and the subject's answers. This results as a set of two-dimensional VWM error measures denoted hereafter as {y (r) ∈ [0, 2π) 2 ; r = 1, . . . , R}. While the subject is performing these tasks, EEG traces are recorded using a Hydrocel GSN equipment with 129 electrodes that are placed on the subject's scalp at specific spatial locations. These electrodes record the electric potential (in micro-volts) over time with a sampling rate of 500 Hz. The subject performs this experiment R = 2400 times. In the following, we denote the set of replicated spatially localized EEG traces as {X S,r t , s ∈ M, r = 1, . . . , R}, where M is the set of electrode coordinates in the two-dimensional plane.

Remark 4.1. The EEG electrodes are spatially localized in 3d space over a template of the human head. Standard practice is to use projected coordinates on the 2d plane. All the information and code to obtain the 2d layout associated with the Hydrocel GSN can be obtained from [START_REF] Oostenveld | Fieldtrip: Open source software for advanced analysis of meg, eeg, and invasive electrophysiological data[END_REF].

Remark 4.2. We developed our method based on realistic modeling assumptions for such real data applications:

1. Since we are interested in studying EEG connectivity and there is empirical evidence for correlations between oscillatory components of brain waves at different frequencies (see [START_REF] Pascual-Marqui | The dual frequency rv-coupling coefficient: a novel measure for quantifying cross-frequency information transactions in the brain[END_REF]), we considered modeling these data as some kind of harmonizable processes. 2. Along the experiment, EEGs corresponds to the electrical activity of sequence of different brain states, rapidly changing from one state to another. For example, the brain states related to visual information acquisition, memorization, joystick usage. . . Piecewise stationary models have been proved useful in such regime/state switching situation [START_REF] Kumar | Classification of seizure and seizure-free eeg signals using multi-level local patterns[END_REF][START_REF] Schröder | Fresped: Frequency-specific change-point detection in epileptic seizure multichannel eeg data[END_REF]. 3. EEG signals represent a sample of a process that is inherently spatial, which justifies a spatial approach. 4. The test subject repeats many times the same experiment. This experiment has a precisely timed performance of different tasks. This is taken into account by our model considering the same distribution of replicates. We additionally introduce short-term dependencies between replicates to account for fatigue and the effect of training. 5. [START_REF] Elul | Gaussian behavior of the electroencephalogram: changes during performance of mental task[END_REF] shows that the Gaussian behavior of EEG is violated most of the time during mental tasks. Therefore, we do not assume Gaussianity in the main results. 6. We use the assumptions of uniform Lipschitz continuity which we find to be mathematically convenient, while at the same time not violating fundamental properties of our real data.

Statistical analysis

To illustrate the application of our method to real data, we proceed in three steps. We have replicated time series associated with two covariates: the orientation and the color errors. Since our estimators are computed on replicated observations, we need to cluster our replicates into meaningful subgroups according to these variables. All replicated time series in a given subgroup will be used to estimate the corresponding spectral quantities. Therefore, the first step of the analysis consists of unsupervised clustering of replicates according to the WVM scores. The second step consists of visualizing the data in order to compare the corresponding spatial timevarying dual-frequency coherence functions within each cluster, and finally the third step is to compare the dual-frequency connectivity networks.

Step 1: clustering with toroidal mixture

Figure 1 shows the bidimensional angular errors for all replicates. Note that both orientation and color errors are well centered around (0, 0), meaning that on average the subject has an unbiased assessment of angle and color. We observe a seemingly more precise quality of memorization for colors than for orientations. Our first step is to model the joint distribution of errors. Using the R package "BAMBI" [START_REF] Chakraborty | BAMBI: Bivariate Angular Mixture Models[END_REF][START_REF] Chakraborty | Bambi: An r package for fitting bivariate angular mixture models[END_REF] and considering the weighted AIC criterion, our best fit is obtained using a two-component mixture of bivariate von Mises distributions. It gives a satisfactory clustering, as shown in Figure 1. The first subgroup of replicates (colored in red) can be interpreted as 'poor' memorization scores, the second subgroup (colored in blue) as 'good' memorization scores. The first subgroup contains approximately 10% of the total number of replicates.

Step 2: estimation of the spatial time-varying dual-frequency coherence

We can now proceed to the estimation of the spatial time-varying dual-frequency coherence functions for each cluster based on formula (2.13). In Remark 4.3. In this analysis and after we consider a size of the Fourier window of 0.5 seconds. The time window was chosen as a rule of thumb. It is an actual research question in this context to choose a proper length for the time window. It should be chosen small enough to avoid bias due to the nonstationarity and large enough to get a suitable frequency resolution.

Step 3: estimation of the dual-frequency functional connectivity networks

Neuroscientists are interested in interpreting significant and sufficiently large coherence values. Hereafter, we consider that coherence values passing above 0.3 are of neurophysiological interest. We use our bootstrap approach to check whether the coherence values are above this reference value by constructing 95% left-sided bootstrap confidence intervals following Section 3. This is done for each pairs of spatial locations (here restricted to a subset of spatial locations of electrodes of interest) and for each time blocks. The block length for the CBB is taken as the integer part of the cubic root of the number of replicates. Next, we construct adjacency matrices of dual-frequency connectivity that refer to different spatial locations at given time points. From these matrices, we construct a dynamic visualization of the network. The resulting networks of (Ω α , Ω β ) dual-frequency functional connectivity at a time point of interest is shown in Figure 3. The graphs present the connectivity estimated from the set of replicates related to poor (left side) and good (right side) scores. Blue lines are drawn between spatial locations for which the lower limit of the bootstrap confidence interval for dual-frequency coherence passes over the predefined threshold value of 0.3.

Conclusions on real data analysis

VWM involves sophisticated functional connections within different areas of the brain, in particular the visual cortex and the prefrontal cortex appear to play fundamental roles [START_REF] Grimault | Oscillatory activity in parietal and dorsolateral prefrontal cortex during rentention in visual short-term memory: additive effects of spatial attention and memory load[END_REF][START_REF] Li | Visual working memory load-related changes in neural activity and functional connectivity[END_REF][START_REF] Barton | Visual working memory in human cortex[END_REF][START_REF] Dai | Eeg cortical connectivity analysis of working memory reveals topological reorganization in theta and alpha bands[END_REF], the latter being involved in encoding task-relevant information ( [START_REF] Lara | The role of the prefrontal cortex in working memory: a mini review[END_REF][START_REF] Funahasi | Working memory in the prefrontal cortex[END_REF]). Interestingly, by examining the dynamics of dual-frequency connectivity during the experiment, our method reveals that brain mechanisms associated with good memorization show a significant correlation between oscillatory components of moderate (alpha) and high frequencies (beta) within the prefrontal cortex during the "memory set stage" (see, Section 4.2). This is illustrated in Figure 3 which captures the connectivity at a specific time moment during the memory stage. This observation appears to be consistent with the current state of understanding of these brain mechanisms, and it provides novel insight that this connectivity is noticeably between the alpha and beta frequency bands.

Conclusions

In this paper, we introduce spectral analysis for a novel model for replicated spatiotemporal processes that are locally time-harmonizable. We propose a consistent estimation procedure for the rescaled spatial time-varying Loève spectrum and the spatial time varying dual frequency coherence. We model dependency across replicated observations and we proved the consistency of the circular block bootstrap. This method allows to obtain valid confidence interval for inference. As an application example, we consider the analysis of replicated measurements of EEG signals in a neuropsychology experiment. We demonstrated the ability of our method to provide a novel way to visualize topographic maps of EEG voltage and to describe the dynamic dual-frequency functional connectivity.

Appendix

In the first section (Section 6.1) of this appendix, we provide results in the Gaussian framework as mentioned in the Remark (2.1). Then, in the Section 6.2 we prove all the results presented in Section 3 and 6.1. In Section 6.3 we provide additional information about the real data experiment.

Gaussian framework

In this section the spatiotemporal random family X S,r : r = 1, . . . , R is assumed to be Gaussian for any and define

S
f S,κ t1,t2 (ω 1 , ω 2 ) := 1 4π 2 t1+N -1 k1=t1-N t2+N -1 k2=t2-N C S,κ (k 1 , k 2 ) e -i(ω1k1-ω2k2) . (6.2) Then C S,κ (t 1 , t 2 ) = π 2 N 2 N -1 j1=-N N -1 j2=-N f S,κ t1,t2 ω N j1 , ω N j2 e i ω N j 1 t1-ω N j 2 t2 ,
where w N j := πj N , j = -N, . . . , N -1 are Fourier frequencies. Moreover,

C S (t 1 , t 2 ) = C S,0 (t 1 , t 2 ) and f S t1,t2 (ω 1 , ω 2 ) = f S,0 t1,t2 (ω 1 , ω 2 ).
In this Gaussian framework we replace conditions (LR) and (MR) by the following condition on f S,κ s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) defined by relation (6.2).

(LGR) There exists some positive constant values L, Q > 0 and a family of functions

f κ : [0, 1] 6 × (-π, π] 2 → C, κ ∈ Z, such that f κ u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) -f κ u 3 ,u 4 ,τ3,τ4 (ω 1 , ω 2 ) ≤ L (∥u 1 -u 3 ∥ + ∥u 2 -u 4 ∥ + |τ 1 -τ 3 | + |τ 2 -τ 4 |) for any u 1 , u 2 , u 3 , u 4 , ∈ [0, 1] 2 , τ 1 , τ 2 ∈ [0, 1] and ω 1 , ω 2 ∈ (-π, π]. f S,κ s 1 ,s 2 ,t1,t2 (ω 1 , ω 2 ) -f κ s1 ,s 2 , ẗ1, ẗ2 (ω 1 , ω 2 ) ≤ Q 1 S 1 + 1 S 2 + 1 T ,
where

s i = (s i,1 , s i,2 ), si = (s i,1 /S 1 , s i,2 /S 2 ), ẗi = t i /T, i = 1, 2 for N ≤ t 1 , t 2 ≤ T -N . In addition, assume that κ∈Z f κ u 1 ,u 2 ,τ1,τ2 ω M l1 , ω M l2 2 < ∞ and lim R→∞ R -1/2 R κ=-R f κ u 1 ,u 2 ,τ1,τ2 ω M l1 , ω M l2 = 0 for ω M li = 2πli M , l i = -M, . . . , M -1, i = 1, 2.
Notice that the first part of the condition (LGR) is a generalization of the condition (L) for replicates in the considered Gaussian case. Since the replicates are not necessarily independent, the replicate-κ-shifted rescaled spatiotemporal Loève spectrum f κ is not necessarily null, and consequently this fact is reflected in the additional superscript κ.

If the condition (LGR) is satisfied, then the condition (LR) is also satisfied with

f κ u,t (ω) = f κ u 1 ,u 3 ,τ1,τ3 (ω 1 , ω 3 )f κ u 2 ,u 4 ,τ2,τ4 (-ω 2 , -ω 4 ) + f κ u 1 ,u 4 ,τ1,τ4 (ω 1 , -ω 4 )f κ u 2 ,u 3 ,τ2,τ3 (-ω 2 , ω 3
). ( 6.3) Theorem 6.1. Assume that conditions (GR), (SR) and (LGR) are fulfilled. Then the conclusions of Theorem 3.1 and of Theorem 3.2 hold.

Proofs

We start with properties of the kernel that are used later in the document to prove properties of our estimators.

Some properties of kernels

Let us recall that the kernel w(•) is bounded and piecewise Lipschitz and supp(w(•) ⊂ [-1, 1]. Thus, there are

k = k w ∈ N, -1 = υ 1 < • • • < υ k = 1 such that w(•) is Lipschitz on each interval (υ j , υ j+1
). This includes the rectangular kernel as well as the triangular kernel. Then

1 Sh S s=1 w u -s/S h = 1 + O 1 Sh 2
uniformly with respect to u such that h ≤ u ≤ 1 -h, provided that h < 0.5, S ∈ N * , and Sh 2 bounded away from 0. We deduce that

s w u (s) = 1 + O 1 S 1 h 2 1 + O 1 S 2 h 2 = 1 + 1 h 2 O 1 S 1 + 1 S 2
uniformly with respect to h ≤ u 1 , u 2 ≤ 1 -h, provided that h < 0.5, S 1 , S 2 ∈ N * , S 1 h 2 as well as S 2 h 2 bounded away from 0.

Proof. First we can write

S s=1 1 Sh w u -s/S h = 1 h 1 0 w u -v h dv + 1 h S s=1 s S s-1 S w u -s/S h -w u -v h dv.

Next notice that for

h ≤ u ≤ 1 -h 1 h 1 0 w u -v h dv = 1.
Moreover, there exists a constant c > 0 such that

w u -s/S h -w u -v h ≤ c Sh
for any h > 0, any S ∈ N * , any s = 1, . . . , S, except a finite number that is bounded by k w and for any v such that (s -1)/S < v < s/S. The constant c and the bound k w do not depend on u, v, h and S. Furthermore in any case we have

w u -s/S h -w u -v h ≤ 2 sup x w(x) < ∞.
Hence we get that

1 h S s=1 s S s-1 S w u -s/S h -w u -v h dv ≤ c Sh 2 + k w × 2 sup x w(x) Sh = O 1 Sh 2 .
This completes the proof of the equalities.

Notation

From now on, for the sake of simplicity, when there is no possibility of confusion, we denote

f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) by f 1,2 (ω 1 , ω 2 ), and f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) by f 1,2 (ω 1 , ω 2 ). Moreover, let us denote d r i (ω) = d r u i ,τi (ω) = t s W τi (t)w u i (s) t+N -1 k=t-N X S,r s,k e -iωk ,
where i = 1, 2. Then the estimator f R 1,2 (ω 1 , ω 2 ) can be equivalently expressed as

f 1,2 (ω 1 , ω 2 ) = 1 R R r=1 I r u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) = 1 4π 2 R R r=1 d r 1 (ω 1 ) d r 2 (ω 2 ), (6.4)
where the space and time smoothed periodogram I r u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) is defined by From Lemma 6.1 and Lemma 6.2 below we easily deduce Theorem 3.1.

I r 1,2 (ω 1 , ω 2 ) = I r u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) := 1 4π 2 d r 1 (ω 1 ) d r 2 (ω 2 ). ( 6
Lemma 6.1. (Limit of the expectation) Let u 1 , u 2 ∈ (0, 1) 2 , τ 1 , τ 2 ∈ (0, 1), ω 1 , ω 2 ∈ [-π, π) and M ≥ 1 fixed. Under condition (L), the expectation

E f 1,2 (ω 1 , ω 2 ) does not depend on R. Moreover for ω M li = liπ M , l i = -M, . . . , M -1, i = 1, 2, we have E f 1,2 ω M l1 , ω M l2 -→ f u 1 ,u 2 ,τ1,τ2 (ω M l1 , ω M l2 )
provided that N = nM as well as

h 2 S 1 , h 2 S 2 , ℏ 2 T → ∞ and n 2 (h + ℏ) → 0 as T, S 1 , S 2 → ∞, h, ℏ → 0 independently of the behavior of n ≥ 1 and R ≥ 1.
Proof. From the definitions of f 1,2 (ω 1 , ω 2 ), d r 1 (ω 1 ) and d r 2 (ω 2 ), we have

E f 1,2 (ω 1 , ω 2 ) = 1 4π 2 R R r=1 E d r 1 (ω 1 ) d r 2 (ω 2 ) = 1 4π 2 R R r=1 t1 t2 s 1 s2 W τ1 (t 1 )W τ2 (t 2 )w u 1 (s 1 )w u 2 (s 2 ) × t1+N -1 k1=t1-N t2+N -1 k2=t2-N E X S,r s 1 ,k1 X S,r s 2 ,k2 e -i(ω1k1-ω2k2) .
From assumption (SR) and relation (2.4) we have

E X S,r s 1 ,k1 X S,r s 2 ,k2 = π 2 N 2 N -1 j1=-N N -1 j2=-N f S s 1 ,s 2 ,k1,k2 ω N j1 , ω N j2 e i(ω N j 1 k1-ω N j 2 k2) ,
where ω N j = jπ N and we deduce that E f 1,2 (ω 1 , ω 2 ) does not depend on R as well as that

E f 1,2 (ω 1 , ω 2 ) = 1 4N 2 t1 t2 s 1 s2 W τ1 (t 1 )W τ2 (t 2 )w u 1 (s 1 )w u 2 (s 2 ) × t1+N -1 k1=t1-N t2+N -1 k2=t2-N N -1 j1=-N N -1 j2=-N f S s 1 ,s 2 ,k1,k2 ω N j1 , ω N j2 e i(ω N j 1 k1-ω N j 2
k2) e -i(ω1k1-ω2k2) .

However, the assumption (L) ensures that

f S s 1 ,s 2 ,k1,k2 ω N j1 , ω N j2 = f s1 ,s 2 , k1, k2 ω N j1 , ω N j2 + O S -1 1 + S -1 2 + T -1 = f u 1 ,u 2 ,τ1,τ2 (ω N j1 , ω N j2 ) + O S -1 1 + S -1 2 + N T -1 + h + ℏ for |τ i -ki | ≤ ℏ + N T and ∥u i -si ∥ ≤ h, i = 1, 2.
Since the supports of the rescaling kernels w(•) and

W (•) are contained in [-1,1], we obtain that E f 1,2 (ω 1 , ω 2 ) = t1 t2 s 1 s2 W τ1 (t 1 )W τ2 (t 2 )w u 1 (s 1 )w u 2 (s 2 ) × 1 4N 2 t1+N -1 k1=t1-N t2+N -1 k2=t2-N N -1 j1=-N N -1 j2=-N f u 1 ,u 2 ,τ1,τ2 (ω N j1 , ω N j2 ) + O S -1 1 + S -1 2 + N T -1 + h + ℏ e i(ω N j 1 -ω1)k1 e -i(ω N j 2 -ω2)k2 .
Notice that the O(•) does not depend on t 1 , t 2 , s 1 , s 2 , k 1 , k 2 , j 1 , j 2 . Using the fact that

N -1 k=-N e iπjk N
= 0 for j ̸ = 0 mod 2N and N = nM with n, M ∈ N * , M being fixed, we deduce that

E f 1,2 ω M l1 , ω M l2 = 1 + O S -1 1 h -2 + S -1 2 h -2 + T -1 ℏ -2 × f u 1 ,u 2 ,τ1,τ2 ω M l1 , ω M l2 + N 2 O S -1 1 + S -1 2 + N T -1 + h + ℏ for 2N T + ℏ ≤ τ 1 , τ 2 ≤ 1 -ℏ -2N -1
T and l 1 , l 2 = -M, . . . , M -1, which concludes the proof of the lemma, noticing that

n 2 S -1 i = n -2 × (n 2 h) 2 × h -2 S -1 i and n 3 T -1 = n -1 × (n 2 ℏ) 2 × ℏ -2 T -1 . .
Then we can easily determine the rate of convergence. Corollary 6.1. (Rate of convergence for the bias) Let u ∈ (0, 1) 2 , τ i ∈ (0, 1), ω M li = liπ M , l i = -M, . . . , M -1, i = 1, 2 and M ≥ 1 fixed. Under condition (L), we have lim

R→∞ √ R E f 1,2 ω M l1 , ω M l2 -f u 1 ,u 2 ,τ1,τ2 ω M li , ω M li = 0 provided that N = nM as well as Rh -4 (S -2 1 +S -2 2 ), Rℏ -4 T -2 , Rn 4 (h 2 +ℏ 2 ) → 0 as R, T, S 1 , S 2 → ∞, h, ℏ → 0, independently of the behavior of n > 1.
This is a direct consequence of the proof of Lemma 6.1, noticing that

Rn 4 S -2 i = Rh -4 S -2 i × Rn 4 h 2 × R -1 h 2 and Rn 6 T -2 = (Rn) -2 × (Rn 4 ℏ 2 ) 2 × Rℏ -4 T -2 .
Lemma 6.2 (Bound for the variance). Assume that the mixing assumption (MR) is satisfied. Then

Var f 1,2 (ω 1 , ω 2 ) ≤ c N 4 R where c is some positive constant independent of R, T, S, N, τ 1 , τ 2 , u 1 u 2 , ω 1 , ω 2 . Proof. First Var f 1,2 (ω 1 , ω 2 ) = 1 R 2 R r1=1 R r2=1 Cov I r1 1,2 (ω 1 , ω 2 ), I r2 1,2 (ω 1 , ω 2 ) .
From assumption (MR)(ii) we have

16π 4 Cov I r1 1,2 (ω 1 , ω 2 ), I r2 1,2 (ω 1 , ω 2 ) = Cov d r1 1 (ω 1 ) d r1 2 (ω 2 ), d r2 1 (ω 1 ) d r2 2 (ω 2 ) ≤ c s 1 s 2 s 3 s 4 W τ1 (t 1 )W τ2 (t 2 )W τ3 (t 3 )W τ4 (t 4 ) ×w u 1 (s 1 )w u 2 (s 2 )w u 3 (s 3 )w u 4 (s 4 ) N 4 α δ 4+δ X (|r 1 -r 2 |).
Then, due to the properties of the window kernels w(•) and W (•), we conclude that

Var f 1,2 (ω 1 , ω 2 ) ≤ c R 2 R r1=1 R r2=1 N 4 α δ 4+δ X (|r 1 -r 2 |) ≤ c N 4 R R-1 κ=-R+1 1 - |κ| R α δ 4+δ X (|κ|).
The lemma is proved under condition (MR)(ii). Following the same reasoning one may easily prove the lemma under condition (MR)(i).

Proof of Proposition 3.1

Proof. First notice that

R 2 Cov f 1,2 (ω 1 , ω 2 ) , f 3,4 (ω 3 , ω 4 ) = 1 16π 4 R r1=1 R r2=1 Cov d r1 1 (ω 1 ) d r1 2 (ω 2 ), d r2 3 (ω 3 ) d r2 4 (ω 4 ) . (6.6)
Let κ = r 1 -r 2 . Then, from the stationarity with respect to the replicates (condition (SR)) we have

Cov d r1 1 (ω 1 ) d r1 2 (ω 2 ), d r2 3 (ω 3 ) d r2 4 (ω 4 ) = t1 • • t4 s 1 • • s 4 4 j=1 W τj (t j )w u j (s j ) × t1+N -1 k1=t1-N • • t4+N -1 k4=t4-N C S,κ s (k) e -i(ω1k1-ω2k2-ω3k3+ω4k4) = π 4 N 4 t1 • • t4 s 1 • • s 4 4 j=1 W τj (t j )w u j (s j ) t1+N -1 k1=t1-N • • t4+N -1 k4=t4-N N -1 l1=-N • • N -1 l4=-N f S,κ s,k ω N l × e i(ω N l 1 -ω1)k1 e -i(ω N l 2 -ω2)k1 e -i(ω N l 3 -ω3)k3 e i(ω N l 4 -ω4)k4
,

where k = (k 1 , k 2 , k 3 , k 4 ), ω N l = (ω N l1 , . . . , ω N l4 ), ω N li = liπ N .
Since condition (LR) is fulfilled and N = nM with n, M ∈ N * , we deduce that

16π 4 Cov d r1 1 ω M p1 d r1 2 ω M p2 , d p2 3 ω M p3 d r2 4 ω M l4 = 1 + O S -1 1 h -2 + S -1 2 h -2 + T -1 ℏ -2 × f κ u,τ ω M p + N 4 O S -1 1 + S -1 2 + N T -1 + h + ℏ
for the Fourier frequencies ω M pi := piπ M , p i = -M, . . . , M -1, i = 1, . . . , 4. Notice that

n 4 S -1 i = (n 2 h) 2 × h -2 S -1 i and n 5 T -1 = n -3 × (n 4 ℏ) 2 × ℏ -2 T -1 .
Then the use of relation (2.9) in the condition (LR) ends the proof of the proposition.

Remark Since the observed process

{X S t } is real-valued, we have f 1,2 (ω 1 , ω 2 ) = f 1,2 (-ω 1 , -ω 2 ). Hence ℜ f 1,2 (ω 1 , ω 2 ) = 1 2 f 1,2 (ω 1 , ω 2 ) + f 1,2 (-ω 1 , -ω 2 ) and ℑ f 1,2 (ω 1 , ω 2 ) = 1 2i f 1,2 (ω 1 , ω 2 ) -f 1,2 (-ω 1 , -ω 2 ) .
Then we can compute the components of the covariance matrix of f 1,2 (ω 1 , ω 2 ):

Var ℜ f 1,2 (ω 1 , ω 2 ) = 1 2 Var f 1,2 (ω 1 , ω 2 ) + ℜCov f 1,2 (ω 1 , ω 2 ) , f 1,2 (-ω 1 , -ω 2 ) , Cov ℜ f 1,2 (ω 1 , ω 2 ) , ℑ f 1,2 (ω 1 , ω 2 ) = 1 2 ℑCov f 1,2 (ω 1 , ω 2 ) , f 1,2 (-ω 1 , -ω 2 ) , Var ℑ f 1,2 (ω 1 , ω 2 ) = 1 2 Var f 1,2 (ω 1 , ω 2 ) -ℜCov f 1,2 (ω 1 , ω 2 ) , f 1,2 (-ω 1 , -ω 2 ) .
These formulas are used to derive the form of the covariance matrix of the limit distribution of √ R f 1,2 (ω 1 , ω 2 ) as R → ∞.

lim R→∞ L √ R f 1,2 (ω 1 , ω 2 ) -E f 1,2 (ω 1 , ω 2 ) = N 2 (0, Σ 2 ) ,
for any Fourier frequencies

ω i = ω M li = liπ M , i = 1, 2, provided that M is fixed and (i) either n is constant, T, S 1 , S 2 , R → ∞, h, ℏ → 0 with h 2 S 1 , h 2 S 2 , ℏ 2 T → ∞; (ii) or n, T, S 1 , S 2 , R → ∞, h, ℏ → 0 with h 2 S 1 , h 2 S 2 , ℏ 2 T → ∞ and n 4 (h + ℏ) → 0 and sup τiT /2≤t≤2τiT t+N -1 k=t-N X S,1 s i ,k ≤ C (6.7)
almost surely, or

sup τiT /2≤t≤2τiT t+N -1 k=t-N E X S,1 s i ,k 4+δ 1/(4+δ) ≤ C (6.8)
for some δ > 0 and for some finite constant C > 0. Here Σ 2 = Σ (ω1,ω2) τ1,τ2,u 1 ,u 2 is the limit covariance matrix of size (2 × 2) (see Proposition 3.1.) Proof. To prove the convergence in distribution of the two-dimensional random vector

√ R f 1,2 (ω 1 , ω 2 ) -E f 1,2 (ω 1 , ω 2 )
we use the Cramér-Rao device. We show the asymptotic normality of

√ R (ζ R -E (ζ R )), where ζ R = a 1 ℜ f 1,2 (ω 1 , ω 2 ) + a 2 ℑ f 1,2 (ω 1 , ω 2 )
for all real numbers a 1 and a 2 . For the sake of simplicity, we take a 1 = 1 and a 2 = 0. Thus,

ζ R = ℜ f 1,2 (ω 1 , ω 2 ) = 1 R R r=1 W R,r ,
where W R,r is the triangular random array of the form W R,r := (2π) -2 ℜ d r 1 (ω 1 ) d r 2 (ω 2 ) , for r = 1, . . . , R, R = 1, 2, . . . . Recall that T , S, ℏ and h depend on R. Moreover, since the replicates {X S,r t }, r = 1, . . . , R, have the same distribution, the random variables W R,r , r = 1, . . . , R have the same distribution. To get the asymptotic normality of

√ R (ζ R -E (ζ R ))
under assumption (MR)(ii), we apply Theorem 3.3.1 from [START_REF] Guyon | Random Fields on a Network[END_REF] (see also [START_REF] Bolthausen | On the central limit theorem for stationary mixing random fields[END_REF]) and hence we verify the following conditions:

(i) sup R E |W R,1 -E (W R,1 )| 2+δ ′ < ∞, for some δ ′ > 0; (ii) lim sup R→∞ R-1 k=-R+1 α W R (k) δ ′ /(2+δ ′ ) < ∞; (iii) R -1 Var R r=1 W R,r -→ σ 2 as R → ∞, where σ 2 is the (1,1)-component of the variance 2 × 2-matrix Σ 2 .
Recall that

d r i (ω) = t s W τi (t)w u i (s) t+N -1 k=t-N X S,r s,k e -iωk .
The relation (i) is a direct consequence of the Hölder inequality, the triangular inequality for metric and the assumption on the moment of order 4 + δ of X S,r s i ,ki . Furthermore, for each R we have that α W R (k) ≤ α X (k) and hence the relation (ii) is a consequence of the assumption k α X (k) δ/(4+δ) < ∞. Finally, thanks to Proposition 3.1, the condition (iii) is fulfilled. The proof under assumption (MR)(i) follows the same reasoning.

Remarks 1) In conditions (6.7) and (6.8), when N = nM → ∞, we essentially need that t i ≥ τ i T /2. The inequality t i ≤ 2τ i T is added to avoid considering any t i between τ i T /2 and ∞.

2) To state Theorem (3.2), it suffices to follow the proof of Proposition (6.1), noticing that the support of the function

t i → W τi (t i ) is contained in [(τ i -ℏ)T, (τ i + ℏ)T ]. Moreover, for N/T + ℏ ≤ τ i ≤ 1 -N/T -ℏ, t i ∈ [(τ i -ℏ)T, (τ i + ℏ)T ] and k i ∈ [t i -N, t i + N -1], we have 0 ≤ k i ≤ T -1. 3) As an example, consider that |X S,1 s,t | ≤ ln(t) -1 for t > 1 a.s. Let 0 < τ < 1 and t ≥ τ T /2 ≥ N + 1 > 1. Then if T > 4/τ 2 and 1 < N < min{ln T, τ T /2} -1, t+N -1 k=t-N X S,1 s,k ≤ 2N ln τ T 2 -N -1 ≤ 4,
and condition (6.7) is satisfied.

6.2.6. Proof of Corollary 3.1

Proof. We apply the delta method (see e.g. [START_REF] Van Der | Asymptotic Statistics[END_REF]) with the function ϕ : ((0, ∞] × R) 2 × R 2 → R defined by ϕ(x 1 , y 1 , x 2 , y 3 , x 3 , y 3 ) :=

x 2 3 +y 2 3 x1x2 . Indeed, ρu 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) = ϕ f , where f = f u 1 ,u 1 ,τ1,τ1 (ω 1 , ω 1 ) ′ , f u 2 ,u 2 ,τ2,τ2 (ω 2 , ω 2 ) ′ , f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) ′ ′
.

Notice that f u j ,u j ,τj ,τj (ω j , ω j ) , j = 1, 2 are real and non-negative. Since the gradient of the function ϕ is equal to for Fourier frequencies ω i,j = ω M li,j with l i,j ∈ {M, . . . , M -1}, i = 1, 2, j = 1, . . . , r. Note that the bootstrap versions f * u 11 ,u 21 ,τ11,τ21 (ω 11 , ω 21 ) , . . . , f * u 1r ,u 2r ,τ1r,τ2r (ω 1r , ω 2r ) are constructed using the same bootstrap blocks (see step 2 of the CBB algorithm). Moreover, let f * ,R 1,2 (ω 1 , ω 2 ) be the bootstrap counterpart of (6.4). Proposition 6.2. Under assumptions of Theorem 3.2 and assuming that b

∇ϕ(x 1 , y 1 , x 2 , y 2 , x 3 , y 3 ) = -2(x 2 3 + y 2 3 ) x 2 1 x 2 , 0, -2(x 2 3 + y 2 3 ) x 1 x 2 2 , 0, 2x 3 x 1 x 2 , 2y 3 x 1 x 2 ′ , Corollary
-1 + R -1 b = o(1) the CBB is consistent i.e., sup x∈R 2 P * √ R f * 1,2 (ω 1 , ω 2 ) -E * f * 1,2 (ω 1 , ω 2 ) ≤ x -P √ R f 1,2 (ω 1 , ω 2 ) -E f 1,2 (ω 1 , ω 2 ) ≤ x p -→ 0 as R -→ ∞,
for ω i = ω M li with l i ∈ {M, . . . , M -1}, i = 1, 2. Proposition 6.3. Under conditions of Proposition 6.2

sup x∈R 2r P * √ R f * -E * f * ≤ x -P √ R f -f ≤ x p -→ 0 as R -→ ∞.
Then Theorem 3.3 is almost a direct application of Theorem 4.1 in [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF] for the smooth function ϕ(x 1 , y 1 , x 2 , y 3 , x 3 , y 3 ) =

x 2 3 +y 2 3 x1x2 and the sequence 1/(2π) 2 dr u 1 ,τ1 (ω 1 ) dr u 2 ,τ2 (ω 2 ), r = 1, . . . , R. We should just be aware that the mentioned theorem cannot be applied directly and requires a small adjustment. Indeed, the considered estimator is assumed to be unbiased in [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF]. This condition does not hold in our context (see equation (6.6) given above in the proof of Lemma 6.1), but one can easily show that the conclusion of Theorem 4.1 in [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF] are yet valid under the conditions of Theorem 3.3.

Proof of Theorem 6.1

Recall that Lemma 6.1 gives us the convergence to 0 of the bias of the estimator. Below we study the behavior of the covariance and the asymptotic normality of the estimator in the Gaussian framework. This restrictive framework allows us to avoid the mixing assumption.

In the following for u i ∈ (0, 1) 2 , τ i ∈ (0, 1) and ω i ∈ [-π, π), i = 1, . . . , 4, denote as

f 1,2 := f u 1 ,u 2 ,τ1,τ2 (ω 1 , ω 2 ) and f 3,4 := f u 3 ,u 4 ,τ3,τ4 (ω 3 , ω 4 ).
where

f 2i-1,2i (ω 2i-1 , ω 2i ) = 1 4π 2 R R r=1 d r 2i-1 d r 2i . and d r j = t s W τj (t)w u j (s) t+N -1 k=t-N X S,r s,k e -iωj k .
From the multilinearity of the cumulants, we have

(4π 2 R) p Cum f 1,2 (ω 1 , ω 2 ), . . . , f 2p-1,2p (ω 2p-1 , ω 2p ) = R r1=1 • • • R rp=1
= ν Cum d i,l : (i, l) ∈ ν 1 × • • • × Cum d i,l : (i, l) ∈ ν q ,
where the summation is over all the indecomposable partitions ν = ν 1 ∪• in [START_REF] Brillinger | Time Series[END_REF] for the definition of indecomposable partitions of a table. Thus, we have

(4π 2 R) p Cum f 1,2 (ω 1 , ω 2 ), . . . , f 2p-1,2p (ω 2p-1 , ω 2p ) = ν R r1=1 • • • R rp=1 Cum d i,l : (i, l) ∈ ν 1 × • • • × Cum d i,l : (i, l) ∈ ν qν . (6.9)
Since the random array d i,l : i = 1, . . . , p; l = 1, 2 is centered Gaussian, the cumulant Cum d i,l : (i, l) ∈ ν k is null except when the set ν k has only two elements. Hence the terms of the sum (6.9) that are not necessarily null correspond to the partitions ν for which all their components ν k , k = 1, . . . , q, have exactly two elements. In this case q = p. The number of such partitions is bounded by (2p)!. Besides we know that

Cum d i1,l1 , d i2,l2 = Cov d i1,l1 , d i2,l2 = E d i1,l1 d i2,l2 = t1 t2 s1 s2 W τ 2i 1 -2+l 1 (t 1 )W τ 2i 2 -2+l 2 (t 2 )w u 2i 1 -2+l 1 (s 1 )w u 2i 2 -2+l 2 (s 2 ) × t1+N -1 k1=t1-N t2+N -1 k2=t2-N e i -(-1) l 1 +1 ωi 1 k1+(-1) l 2 +1 ωi 2 k2 E X ri 1 s 1 ,k1 X ri 2 s 2 ,k2 for 2N T + ℏ ≤ τ 2i1-2+l1 , τ 2i2-2+l2 ≤ 1 -ℏ -2N -1 T . However, E X ri 1 s 1 ,k1 X ri 2 s 2 ,k2 = π 2 N 2 N -1 g1=-N N -1 g2=-N f S,ri 1 -ri 2 s 1 ,s 2 ,k1,k2 ω N g1 , ω N g2 e i ω N g 1 k1-ω N g 2 k2 .
Moreover, under assumption (LGR) we get that

f S,ri 1 -ri 2 s 1 ,s 2 ,k1,k2 ω N g1 , ω N g2 = f ri 1 -ri 2 s1 ,s 2 , k1, k2 ω N g1 , ω N g2 + O(S -1 1 + S -1 2 + T -1 = f ri 1 -ri 2 u 1 ,u 2 ,τ1,τ2 ω N g1 , ω N g2 + O S -1 1 + S -1 2 + N T -1 + h + h for |τ i -ki | ≤ ℏ -N T and ∥u i -si ∥ ≤ h, i = 1, 2.
When ω i1 and ω i2 are the Fourier frequencies, say ω i1 = ω M pi 1

and

ω i2 = ω M pi 2 , N = nM , we obtain that Cum d i1,l1 , d i2,l2 = Cov d i1,l1 , d i2,l2 = E d i1,l1 d i2,l2 = t1 t2 s1 s2 W τ 2i 1 -2+l 1 (t 1 )W τ 2i 2 -2+l 2 (t 2 )w u 2i 1 -2+l 1 (s 1 )w u 2i 2 -2+l 2 (s 2 ) × π 2 N 2 N -1 g1=-N N -1 g2=-N t1+N -1 k1=t1-N t2+N -1 k2=t2-N f ri 1 -ri 2 u 2i 1 -2+l 1 ,u 2i 2 -2+l 2 ,τ 2i 1 -2+l 1 ,τ 2i 2 -2+l 2 ω N g1 , ω N g2 × e i ω N j 1 +(-1) l 1 ωi 1 k1 e -i ω N j 2 +(-1) l 2 ωi 2 k2 + N 2 O(S -1 1 + S -1 2 + N T -1 + h + h = 1 + O S -1 1 h -2 1 + O S -1 2 h -2 1 + O T -1 ℏ -2 × 4π 2 f ri 1 -ri 2 i1,l1,i2,l2 + N 2 O S -1 1 + S -1 2 + N T -1 + h + h ,
where for simplicity we denote

f ri 1 -ri 2 i1,l1,i2,l2 := f ri 1 -ri 2 u 2i 1 -2+l 1 ,u 2i 2 -2+l 2 ,τ 2i 1 -2+l 1 ,τ 2i 2 -2+l 2 (-1) l1-1 ω M pi 1 , (-1) l2-1 ω M pi 2 .
Now fix a partition ν which is significant for the sum (6.9), that is, ν = ν 1 ∪ • • • ∪ ν p is an indecomposable partition of Table 1 and each of its components has two elements.

First, we can state that there exists a path starting at (i 1 , l 1 ) ∈ ν 1 visiting only once every (i j , l j ) of the table and such that (i j , l j ) and (i j , 3 -l j ) do not belong to the same component ν k of ν. More precisely, we can build a sequence (i 1 , l 1 ), (i 1 , 3 -l 1 ), (i 2 , l 2 ), (i 2 , 3 -l 2 ), (i 3 , l 3 ), . . . , (i j , l j ), (i j , 3 -l j ), . . . , (i p , l p ), (i p , 3 -l p ) of all the elements of Table 1 with ν 1 = {(i p , 3 -l p ), (i 1 , l 1 )}, ν kj = {(i j-1 , 3 -l j-1 ), (i j , l j )}, for j = 2, . . . , p and where {k 2 , . . . , k p } is some permutation of {2, . . . , p}. See the proof of Lemma 6.3 below. Then, we deduce that

Complement : Indecomposable partition

Let ν = ν 1 ∪ • • • ∪ ν p be an indecomposable partition of Table 1, such #ν 1 = • • • = #ν p = 2. We can state the following elementary properties 1. If ν k ↔ ν k ′ (ν k and ν k ′ hook, and k ̸ = k ′ ) then ∃(i 1 , l 1 ) ∈ ν k , (i 2 , l 2 ) ∈ ν k ′ such that i 1 = i 2 . Of course in this case, l 1 ̸ = l 2 and, since l 1 , l 2 ∈ {1, 2}, we have l 2 = 3 -l 1 . 2. If p ≥ 2, (i, 1) ∈ ν k and (i, 2) and(i, l) ∈ q j=1 ν kj ⇒ (i, 3 -l) ∈ q j=1 ν kj then q = p and ν = p j=1 ν kj . Lemma 6.3. There exists at least one path passing through each (i, l) of Table 1 and only once.

∈ ν k ′ then k ̸ = k ′ . ν k ∩ ν ′ k = ∅. 3. If ν k1 ↔ ν k2 ↔ • • • ↔ ν kq ,
Proof. We build such a path by "recurrence", applying the previous properties.

(1) Let k 1 = 1 and (i 1 , l 1 ) ∈ ν k1 .

(2) Since 2 ≤ p, let k 2 ≤ p for which (i 1 , 3 -l 1 ) ∈ ν k2 . Then k 2 ̸ = k 1 and there exists a unique (i 2 , l 2 ) ∈ ν k2 such that i 2 ̸ = i 1 . (3) If 3 ≤ p, let k 3 ≤ p for which (i 2 , 3 -l 2 ) ∈ ν k3 . Then k 3 / ∈ {k 1 , k 2 } and there exists a unique (i 3 , l 3 ) ∈ ν k3 such that i 3 ̸ = i 2 . We have also i 3 / ∈ {i 1 , i 2 }. . . . (j) If 3 ≤ j ≤ p, let k j ≤ p for which (i j-1 , 3 -l j-1 ) ∈ ν kj . Then k j / ∈ {k 1 , . . . , k j-1 }, and let (i j , l j ) ∈ ν kj such that i j / ∈ {i 1 , . . . , i j-1 }. Notice that k j and (i j , l j ) are unique. . . . (p) When j = p, we see that k p is the only value of {1, . . . , p}, which have not yet been considered: {k p } = {1, . . . , p} \ {k 1 , . . . , k p-1 }. Moreover, ν = j=1 ν kj and (i p , 3 -l p ) ∈ ν k1 .

Thus we have built the sequence:

(i 1 , l 1 ), (i 1 , 3 -l 1 ), (i 2 , l 2 ), (i 2 , 3 -l 2 ), (i 3 , l 3 ), . . . , (i j , l j ), (i j , 3 -l j ), . . . , (i p , l p ), (i p , 3 -l p )

where ν k1 = {(i p , 3 -l p ), (i 1 , l 1 )} and ν kj = {(i j-1 , 3 -l j-1 ), (i j , l j )}, j = 2, . . . , p. Hence Lemma 6.3 is proved.

Visual Working Memory performance experiment

  for some δ > 0. The mixing coefficients are defined as follows α X (κ) := sup r sup S sup A∈Fr (S) B∈F r+κ (S) P(A ∩ B) -P(A)P(B) ,

Fig 1 .

 1 Fig 1. Angular errors associated with each replicate; x-axis: orientation error; y-axis: color error. Red and blue colors identify the subgroup resulting from unsupervised clustering.

Figure ( 2 )

 2 we present the estimated spatial dual-frequency coherence for frequency bands (Ω θ , Ω α ) in the group of poor VWM scores ρ τ,τ,u 1 ,u 2 (Ω θ , Ω α ) at time τ = 1.2s. The graph on the left shows the location of u 1 ∈ M. It is specified by the user. The graph on the right contains the output of our software, i.e. a topographic map of the spatial coherence ρ τ,τ,u 1 ,u 2 (Ω θ , Ω α ) , u 2 ∈ [0, 1] 2 .

Fig 2 .

 2 Fig 2. Estimated spatial (Ω θ , Ωα) dual-frequency coherence function associated with poor VWM scores. On the left-hand side: location of u 1 (user-specified); on the right-hand side: the estimated spatial coherence for all spatial locations u 2 .

Fig 3 .

 3 Fig 3. Statistically and neurophysiologically significant functional connections at (Ωα, Ω β ) frequency bands associated with poor (left graph) and with good memorization scores (right graph).

∈ N * 2

 2 and any R ∈ N * . Under the stationarity condition (SR), we have Cov X ′ t2for any κ ∈ Z and for any positive integers r and r ′ > -κ. Denote C S,κ (t 1 , t 2 ) := Cov X

. 5 ) 6 . 2 . 3 .

 5623 Proof of Theorem 3.1

6. 2 . 5 .

 25 Proof of Theorem 3.2 Theorem 3.2 is a direct application of the Cramér-Wold device for the univariate asymptotic normality result, which we derive from Corollary 6.1 and from Proposition 6.1 shown below. Proposition 6.1. Let assumptions (SR), (MR) and (LR) hold. Then

T -1 ℏ - 2 × 4π 2 f 1 ip, 3 - 2 .R p/ 2

 221322 Cum d i,l : (i, l) ∈ ν 1 × • • • × Cum d i,l : (i, l) ∈ ν p Cum d i,l : (i, l) ∈ ν k1 × • • • × Cum d i,l : (i, l) ∈ ν kp ≤ ri p -ri 1 ip,3-lp,i1,l1 + N 2 O S -1 1 + S -1 2 + N T -1 + h + h × 4π 2 f ri 1 -ri 2 i1,3-l1,i2,l2 + N 2 O S -1 1 + S -1 2 + N T -1 + h + h × • • • • • • × 4π 2 f ri p-1 -ri p ip-1,3-lp-1,ip,lp + N 2 O S -1 1 + S -1 2 + N T -1 + h + h . Setting κ j = r ij-1 -r ij , j = 2, . . . , p, p ≥ 3, we can write lp,i1,l1 × f ri 1 -ri 2 i1,3-l1,i2,l2 × • • • × f ri p-1 -ri p ip-1,3-lp-1,ip,lp lp,i1,l1 × f κ2 i1,3-l1,i2,l2 × • • • × f κp ip-1,3-lp-1,ip,lpUnder the condition (LGR) (convergence assumptions) we deduce that limR→∞ Cum d i,l : (i, l) ∈ ν 1 × • • • × Cum d i,l : (i, l) ∈ ν p = 0 for any indecomposable partition ν = ν 1 ∪ • • • ∪ ν p , with #ν 1 = • • • = #ν p = 2, p ≥ 3.Finally, using the relation (6.9) and the Gaussianity of the observations we get that lim R→∞ Cum f 1,2 (ω 1 , ω 2 ), . . . , f 2p-1,2p (ω 2p-1 , ω 2p ) = 0 for any integer p ≥ 3 and any Fourier frequencies ω 1 = ω M l1 , ω 2 = ω M l2 , . . . , ω 2p = ω M l2p , l 1 , l 2 , . . . , l 2p = -M, . . . , M -1. This completes the proof of Theorem 6.1.

Figure 4

 4 Figure 4 illustrates the experiment performed. It shows two rectangles with arrows. On the left panel, you can see the examples of arrows appearing on the screen, which can have different orientations and colors. The test subject has to memorize them. The right panel shows a possible answer of the test subject. The graph at the bottom presents the timeline of the experiment described in Section 4. The experiment involved 6 participants, each of whom performed 2400 repetitions. For simplicity, in our illustrative data analysis, we considered only one subject. The aim of this experiment is to possibly identify in the EEG traces, recorded during the memory set step, specific brain mechanisms that could be related to the errors committed.

Fig 4 .

 4 Fig 4. Schematic representation of the experiment.

Figure 5

 5 Figure5shows the software interface through which the subject answers questions about the color he has memorized. He indicates the color from among the continuous color scale wrapped on a circle. As a result, the error made by the test subject is measured as the angle between the truth and his answer.

Fig 5 .

 5 Fig 5. Software interface for color specification.

  3.1 is a direct consequence of Theorem 3.1 in[START_REF] Van Der | Asymptotic Statistics[END_REF] and Theorem 3.2.Below we state bootstrap consistency and its multivariate version, which are direct application of Theorem 3.2 from[START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF]. Recall that by P * and E * we denote the conditional probability and conditional expectation given the sample. Let f * = f * u 11 ,u 21 ,τ11,τ21 (ω 11 , ω 21 ) , . . . , f * u 1r ,u 2r ,τ1r,τ2r (ω 1r , ω 2r ) f * u 11 ,u 21 ,τ11,τ21 (ω 11 , ω 21 ) , . . . , E * f * u 1r ,u 2r ,τ1r,τ2r (ω 1r , ω 2r )

	6.2.7. Proof of Theorem 3.3
	′

′

and

E * f * = E *

Table 1

 1 • •∪ν qν of Table 1. See also Theorem 2.3.2

	(1,1)	(1,2)
	. . .	. . .
	(p,1) (p,2)

Acknowledgements

We thank the editor and the anonymous reviewers for their comments, which helped us to improve the manuscript.

* Anna Dudek acknowledges support from the King Abdullah University of Science and Technology (KAUST) Research Grant OSRare funded by James S. McDonnell Foundation 21st Century Science Initiative in Understanding Human Cognition (grant number 220020370; received by Denes Szucs)

Data Availability statement

The data that support the findings of this study are available on request from Prof. Denes Szucs. The data are not publicly available due to privacy and ethical restrictions.

Asymptotic behavior of the covariance and consistency

First notice that

Using the fact that the observations are Gaussian we know that

Additionally, from the stationarity condition (SR) with respect to the replications, we have

Consequently,

T and ω i = ω M li = liπ M , i = 1, 2, 3, 4. Then using Lemma 6.1 we obtain the consistency of the estimator f 1,2 (ω 1 , ω 2 ). Moreover, we get the convergence of R Cov f 1,2 (ω 1 , ω 2 ) , f 3,4 (ω 3 , ω 4 ) . This concludes the proof.

Proof of the asymptotic normality

By Lemma P4.5 in [START_REF] Brillinger | Time Series[END_REF], it remains to prove that every cumulant of order p ≥ 3 of √ R f converges to 0: lim R→∞ R p/2 Cum f 1,2 (ω 1 , ω 2 ), . . . , f 2p-1,2p (ω 2p-1 , ω 2p ) = 0,