# Adaptive Frequency Band Analysis for Functional Time Series

Pramita<EMAIL_ADDRESS>[ Scott<EMAIL_ADDRESS>[
Department of Statistics, George Mason University

###### Abstract

The frequency-domain properties of nonstationary functional time series often
contain valuable information. These properties are characterized through its
time-varying power spectrum. Practitioners seeking low-dimensional summary
measures of the power spectrum often partition frequencies into bands and
create collapsed measures of power within bands. However, standard frequency
bands have largely been developed through manual inspection of time series
data and may not adequately summarize power spectra. In this article, we
propose a framework for adaptive frequency band estimation of nonstationary
functional time series that optimally summarizes the time-varying dynamics of
the series. We develop a scan statistic and search algorithm to detect changes
in the frequency domain. We establish theoretical properties of this framework
and develop a computationally-efficient implementation. The validity of our
method is also justified through numerous simulation studies and an
application to analyzing electroencephalogram data in participants alternating
between eyes open and eyes closed conditions.

62M15,

62G20,

Frequency band estimation,

functional time series,

locally stationary,

multitaper estimation,

spectrum analysis.,

###### keywords:

[class=MSC2020]

###### keywords:

and

##  1 Introduction

Functional data has emerged as an important object of interest in statistics
in recent years as advances in technology have led to an abundance of high-
dimensional and high resolution data. While classical statistical methods
often fail in this setting, functional data analysis techniques use the smooth
structure of the observed process in order to model non-sparse, high-
dimensional and high-resolution data. The term “functional data analysis” was
coined by [26] and [27], but the history of this area is much longer, dating
back to [14] and [29]. The intrinsic high, or rather infinite, dimensionality
of such data poses interesting challenges in both theory and computation and
has garnered a considerable amount of attention within the statistics
community. Functional time series data often arise in many important problems,
such as analyzing forward curves derived from commodity futures [16], daily
patterns of geophysical and environmental data [31], demographic quantities,
such as age-specific fertility or mortality rates studied over time [32], and
neurophysiological data, such as electroencephalography (EEG) and functional
magnetic resonance imaging (fMRI), recorded at various locations in the brain
[13, 33]. For example, NASA records surface temperatures for more than 5000
locations, and these readings are used to identify yearly temperature
anomalies, which are crucial in studying global warming patterns [19]. In
practice, such data are typically analyzed separately for each location, which
is computationally expensive and fails to account for the spatial structure of
the data. However, it is reasonable to assume temperatures vary smoothly
across locations, and thus ideal to analyze the collection of readings across
locations as a functional data object.

The frequency-domain properties of time series data, including the
aforementioned functional time series data examples, often contain valuable
information. These properties are characterized through its power spectrum,
which describes the contribution to the variability of a time series from
waveforms oscillating at different frequencies. Practitioners seeking
practical, low-dimensional summary measures of the power spectrum often
partition frequencies into bands and create collapsed measures of power within
these bands. In practice, frequency band summary measures are used in a wide
variety of contexts, such as to summarize seasonal patterns in environmental
data, and to measure association between frequency-domain characteristics of
EEG data and cognitive processes [17]. In the scientific literature, standard
frequency bands used for analysis have largely been developed through manual
inspection of time series data. This is accomplished by noting prominent
oscillatory patterns in the data and forming frequency bands that largely
account for these dominant waveforms. For example, frequency-domain analysis
of heart rate variability (HRV) data began in the late 1960s [5] and led to
the development of three primary frequency bands used to summarize power
spectra: very low frequency (VLF) ($\leq$ 0.04 Hz), low frequency (LF)
(0.04-0.15 Hz), and high frequency (HF) (0.15-0.4 Hz) [20]. Within these
frequency bands, collapsed measures of power are used to summarize the
frequency-domain properties of the data and provide a basis for comparison.
However, fixed bands that are not allowed to vary across subjects, covariates,
experimental variables, or other settings may not adequately summarize the
power spectrum. For example, [18] shows that in the study of EEG data,
different frequency bands within the fixed alpha band (8-13 Hz) reflect quite
different cognitive processes. This indicates that the _number_ of frequency
bands needed to adequately summarize the power spectrum may depend on
experimental factors. [11] also notes the _endpoints_ for the alpha frequency
band may vary across subjects and proposes a data-adaptive definition for
alpha band power. Both of these examples illustrate the need for a
standardized, quantitative approach to frequency band estimation that provides
a data-driven determination of the number of frequency bands and their
respective endpoints.

For univariate nonstationary time series, a data-adaptive framework for
identifying frequency bands that best preserve time-varying dynamics has been
introduced in [7]. They develop a frequency-domain scan and hypothesis testing
procedure and search algorithm to detect and validate changes in nonstationary
behavior across frequencies. A rigorous framework for studying the frequency
domain properties of functional time series was established in [23] and [24].
The spectral characteristics of locally stationary functional time series have
been developed by [39] and have been subsequently studied in several papers
including [1], [38] and [37]. In this paper, we develop a scan statistic for
identifying the optimal frequency band structures for characterizing
functional time series data. Developing appropriate scan statistics for the
functional domain poses two important challenges. As the periodogram itself is
an inconsistent estimator of the power spectrum, the univariate scan statistic
was developed using a local multitapered periodogram estimator [7]. The
corresponding test for detecting changes in the time-varying dynamics of the
power spectrum across frequencies uses as asymptotic $\chi^{2}$ distribution
approximation of the local multitapered periodograms. However, the asymptotic
behavior of the multitaper periodogram estimator for functional time series is
not well studied. In fact, the approximation of such a quantity by a
$\chi^{2}$ distribution is not possible in the functional case. Hence, a
functional generalization of this result requires a completely different
approach. In this paper, we derive a central limit theorem type result for the
functional multitaper periodogram and establish a uniform distributional
approximation of the functional scan statistics to a quadratic functional of a
Gaussian process. The second challenge is implementation of the developed
method, which is challenging largely due to the non-standard limit
distribution and high-dimensionality of the data. In this work, we propose an
efficient algorithm to detect the frequency band structure based on a finite
dimensional projection of the functional scan statistic. Moreover, we propose
a computationally-efficient and memory-smart modification of the algorithm
based on the intrinsic structure of the asymptotic distribution. These
modifications enable the proposed methodology to be applied to longer and more
densely-observed functional time series typically encountered in practice.

The rest of the paper is organized as follows. Section 2 introduces the
locally stationary functional time series model and offers a definition of the
frequency-banded power spectrum. Section 3 first provides an overview of the
local multitaper periodogram estimator for the power spectrum, which is used
as a first step in the proposed procedure. It should be noted that the
proposed procedure can be carried out with any consistent time-varying
spectral estimator, but this work focuses on the local multitaper periodogram
estimator due to favorable empirical and theoretical properties established
herein. Section 3 then details the components of the proposed analytical
framework, including the scan statistics and their theoretical properties, an
iterative algorithm for identifying multiple frequency bands, and another test
statistic to determine if the power spectrum within a frequency band is
stationary with respect to time. Section 4 contains simulation results to
evaluate empirical performance in estimating the number of frequency bands and
their endpoints and provides an application to frequency band analysis of EEG
data for an individual alternating between eyes open (EO) and eyes closed (EC)
resting states. Section 5 offers a discussion of the results and concluding
remarks. Proofs for all theoretical results are provided in the Appendix of
this paper.

##  2 Nonstationary Functional Time Series in Frequency Domain

###  2.1 Notation and Functional Set-up

We begin by introducing some notation and results from functional analysis
used in the paper. We typically assume observations are elements of a
separable Hilbert space, $H$. That is, $H$ is equipped with an inner product
$\langle\cdot,\cdot\rangle$ and an associated norm $\|\cdot\|$ defined as
$\|x\|=\langle x,x\rangle$ for $x\in H$. Let $\mathcal{L}$ be the space of
bounded linear operators on $H$ with the norm

$\displaystyle\|\Phi\|_{\mathcal{L}}=\sup\\{\|\Phi(x)\|:\|x\|\leq 1\\}.$

An operator $\Phi\in\mathcal{L}$ is compact if there exists two orthonormal
bases, $\\{v_{j}\\}$ and $\\{f_{j}\\}$, and a real sequence converging to
zero, $\\{\lambda_{j}\\}$, such that for $x\in H$,

$\displaystyle\Phi(x)=\sum_{j=1}^{\infty}\lambda_{j}\langle x,v_{j}\rangle
f_{j}.$

This representation is called the singular value decomposition. A compact
operator with such a representation is said to be a Hilbert-Schmidt operator
if $\sum_{j=1}^{\infty}\lambda_{j}^{2}<\infty$. The space $\mathcal{S}$ of
Hilbert-Schmidt operators is itself a separable Hilbert space with inner
product

$\displaystyle\langle\Phi_{1},\Phi_{2}\rangle_{\mathcal{S}}=\sum_{i=1}^{\infty}\langle\Phi_{1}(e_{i}),\Phi_{2}(e_{i})\rangle$

where $\\{e_{i}\\}$ is an arbitrary orthonormal basis whose selection does not
alter the value of $\langle\Phi_{1},\Phi_{2}\rangle_{\mathcal{S}}$. One can
show that $\|\Phi\|_{\mathcal{S}}^{2}=\sum_{j=1}^{\infty}\lambda_{j}^{2}$ and
$\|\Phi\|_{\mathcal{L}}\leq\|\Phi\|_{\mathcal{S}}$. An operator
$\Phi\in\mathcal{L}$ is symmetric if $\langle\Phi(x),y\rangle=\langle
x,\Phi(y)\rangle$ for all $x,y\in H$ and positive semi-definite if
$\langle\Phi(x),x\rangle\geq 0$ for all $x\in H$.

In particular, we are interested in the Hilbert space
$L^{2}([0,1]^{k},\mathbb{C})$, the set of complex-valued measurable functions
$x$ defined on $[0,1]^{k}$ satisfying $\int_{[0,1]^{k}}|x^{2}(u)|du<\infty$,
equipped with the inner product

$\displaystyle\langle x,y\rangle=\int_{[0,1]^{k}}x(t)\overline{y(t)}dt,$

where $\overline{y}$ denotes the complex conjugate of $y\in\mathbb{C}$. Note
for $x,y\in L^{2}([0,1]^{k},\mathbb{C})$, we denote by $x=y$ the fact that
$\|x-y\|=\int_{[0,1]^{k}}|x(t)-y(t)|^{2}dt=0$. The space of square integrable
real valued functions $L^{2}([0,1]^{k},\mathbb{R})$ is defined similarly.

An important class of integral operators on $L^{2}([0,1]^{k},\mathbb{C})$ are
those defined by

$\displaystyle\Phi[x](t)=\int_{[0,1]^{k}}\phi(t,s)x(s)ds$

for $x\in L^{2}([0,1]^{k},\mathbb{C})$ with kernel function $\phi\in
L^{2}\left([0,1]^{k}\times[0,1]^{k},\mathbb{C}\right)$. These operators are
Hilbert-Schmidt if and only if $\int\int|\phi^{2}(s,t)|dsdt<\infty,$ in which
case,

$\displaystyle\|\Phi\|_{\mathcal{S}}^{2}=\int_{[0,1]^{k}}\int_{[0,1]^{k}}|\phi^{2}(s,t)|dsdt.$

If $\phi(s,t)=\phi(t,s)$ and $\int_{0}^{1}\int_{0}^{1}\phi(s,t)^{2}dsdt\geq
0$, the integral operator $\Phi$ is symmetric and positive semi-definite.

###  2.2 Power Spectrum of Locally Stationary Functional Time Series

Let $\\{X_{h}\\}_{h\in\mathbb{Z}}$ be a weakly stationary functional time
series such that $X_{t}$ is a random element of $L^{2}([0,1],\mathbb{R})$ for
each $t\in\mathbb{Z}$ with expectation $\mu$ and auto-covariance kernel
$r_{h}\in L([0,1]^{2},\mathbb{R})$ at lag $h$ defined as

$r_{h}(s,t)=\mathbb{E}\left((X_{h}(t)-\mu(t))(X_{0}(s)-\mu(s))\right),~{}~{}~{}~{}~{}~{}~{}~{}~{}t,s\in[0,1].$

Let $\mathcal{R}_{h}$ be the corresponding autocovariance operator, i.e., the
integral operator induced by the autocovariance kernel $r_{h}$. The second
order dynamics of this time series can be completely described by the spectral
density operator, defined as the Fourier transform of
$\\{\mathcal{R}_{h}\\}_{h\in\mathbb{Z}}$, acting on $L^{2}([0,1],\mathbb{C})$,
i.e.,

$\mathcal{F}_{\omega}=\sum_{h\in\mathbb{Z}}\mathcal{R}_{h}e^{-\mathrm{i}2\pi\omega
h}\quad\omega\in(0,0.5).$

We will assume our data are centered and hence $\mu=0$. When stationarity is
violated, we can no longer define the spectral density in this manner for all
time points. To meaningfully model nonstationarity, we consider a triangular
array $\\{X_{t,T}:1\leq t\leq T\\}_{T\in\mathbb{N}}$ as a doubly-indexed
functional time series, where $X_{t,T}$ is a random element with values in
$L^{2}([0,1],\mathbb{R})$ for each $1\leq t\leq T$ and $T\in\mathbb{N}$. The
processes $\\{X_{t,T}:1\leq t\leq T\\}$ are extended on $t\in\mathbb{Z}$ by
setting $X_{t,T}=X_{1,T}$ for $t<1$ and $X_{t,T}=X_{T,T}$ for $t>T$. The
sequence of stochastic processes $\\{X_{t,T}:t\in\mathbb{Z}\\}$ indexed by
$T\in\mathbb{N}$ is called locally stationary if for all rescaled times
$u\in[0,1]$, there exists an $L^{2}([0,1],\mathbb{R})$-valued strictly
stationary process $\\{X^{(u)}_{t}:t\in\mathbb{Z}\\}$ such that

$\Bigl{\|}X_{t,T}-X^{(u)}_{t}\Bigr{\|}_{2}\leq\big{(}\big{|}\tfrac{t}{T}-u\big{|}+\tfrac{1}{T}\big{)}\,P_{t,T}^{(u)}\qquad
a.s.$ (2.2)

for all $1\leq t\leq T$, where $P_{t,T}^{(u)}$ is a positive real-valued
process such that for some $\rho>0$ and $C<\infty$ the process satisfies
$\mathbb{E}\big{(}\big{|}P_{t,T}^{(u)}\big{|}^{\rho}\big{)}<C$ for all $t$ and
$T$ and uniformly in $u\in[0,1]$. If the second-order dynamics change
gradually over time, the second order dynamics of the stochastic process
$\\{X_{t,T}:t\in\mathbb{Z}\\}_{T\in\mathbb{N}}$ are completely described by
the time-varying spectral density kernel given by

$\displaystyle
f_{u,\omega}(\tau,\sigma)=\sum_{h\in\mathbb{Z}}\mathbb{E}\left(X_{t+h}^{(u)}(\tau)X_{t}^{(u)}(\sigma)\right)e^{-\mathrm{i}2\pi\omega
h},~{}~{}~{}~{}~{}~{}~{}\omega\in(0,1/2).$ (2.3)

The time varying spectral density operator $F_{u,\omega}$ is the operator
induced by $f_{u,\omega}$ by right integration. Following the set-up of [1],
in order to establish our asymptotic results, we impose the following
technical assumptions on the time series under consideration.

######  Assumption 2.1.

Assume $\\{\,X_{t,T}\colon t\in\mathbb{Z}\\}_{T\in\mathbb{N}}$ is a locally
stationary zero-mean stochastic process and let
$~{}\kappa_{{k};t_{1},\ldots,t_{{k}-1}}$ be a positive sequence in
$L^{2}([0,1]^{k},\mathbb{R})$ independent of $T$ such that, for all
$j=1,\ldots,k-1$ and some $\ell\in\mathbb{N}$,

$\displaystyle\sum_{t_{1},\ldots,t_{k-1}\in\mathbb{Z}}(1+|t_{j}|^{\ell})\|\kappa_{{k};t_{1},\ldots,t_{{k}-1}}\|_{2}<\infty.$
(2.4)

Let us denote

$Y^{(T)}_{t}=X_{t,T}-X_{t}^{(t/T)}\qquad\text{and}\qquad
Y_{t}^{(u,v)}=\frac{X^{(u)}_{t}-X^{(v)}_{t}}{(u-v)}$ (2.5)

for $T\geq 1$, $1\leq t\leq T$ and $u,v\in[0,1]$ such that $u\neq v$. Suppose
furthermore that $k$-th order joint cumulants satisfy

  1. (i)

$\|\text{Cum}(X_{t_{1},T},\ldots,X_{t_{k-1},T},Y^{(T)}_{t_{k}})\|_{2}\leq\frac{1}{T}\|\kappa_{k;t_{1}-t_{k},\ldots,t_{k-1}-t_{k}}\|_{2}$,

  2. (ii)

$\|\text{Cum}(X_{t_{1}}^{(u_{1})},\ldots,X_{t_{k-1}}^{(u_{k-1})},Y_{t_{k}}^{(u_{k},v)})\|_{2}\leq\|\kappa_{k;t_{1}-t_{k},\ldots,t_{k-1}-t_{k}}\|_{2}$,

  3. (iii)

$\sup_{u}\|\text{Cum}(X_{t_{1}}^{(u)},\ldots,X_{t_{k-1}}^{(u)},X_{t_{k}}^{(u)})\|_{2}\leq\|\kappa_{k;t_{1}-t_{k},\ldots,t_{k-1}-t_{k}}\|_{2}$,

  4. (iv)

$\sup_{u}\|\frac{\partial}{\partial
u}\text{Cum}(X_{i,t_{1}}^{(u)},\ldots,X_{i,t_{k-1}}^{(u)},X_{i,t_{k}}^{(u)})\|_{2}\leq\|\kappa_{k;t_{1}-t_{k},\ldots,t_{k-1}-t_{k}}\|_{2}$.

Under these assumptions, the spectral density operator
$\mathcal{F}_{u,\omega}$ is a Hilbert-Schmidt operator and the spectral
density kernel $f_{u,\omega}\in L^{2}\left([0,1]^{2},\mathbb{C}\right)$ is
twice-differentiable with respect to $u$ and $\omega$.

We consider the demeaned power spectrum,

$g_{u,\omega}(\tau,\sigma):=f_{u,\omega}(\tau,\sigma)-\int_{0}^{1}f_{u,\omega}(\tau,\sigma)du,$
(2.6)

and equivalently
$\mathcal{G}_{u,\omega}=\mathcal{F}_{u,\omega}-\int_{0}^{1}\mathcal{F}_{u,\omega}du.$
We assume $g$ admits a partition of the frequency space such that

$g_{u,\omega}(\tau,\sigma)=\left\\{\begin{array}[]{ll}g^{(1)}_{u}(\tau,\sigma)&~{}~{}\text{for}~{}\omega\in(0,\omega_{1}),\\\
g^{(2)}_{u}(\tau,\sigma)&~{}~{}\text{for}~{}\omega\in[\omega_{1},\omega_{2}),\\\
~{}~{}~{}~{}\vdots\\\
g^{(p)}_{u}(\tau,\sigma)&~{}~{}\text{for}~{}\omega\in[\omega_{p-1},0.5).\end{array}\right.$
(2.7)

Equivalently, the spectral density operator $\mathcal{G}$ admits a partition

$\mathcal{G}_{u,\omega}=\left\\{\begin{array}[]{ll}\mathcal{G}^{(1)}_{u}&~{}~{}\text{for}~{}\omega\in(0,\omega_{1}),\\\
\mathcal{G}^{(2)}_{u}&~{}~{}\text{for}~{}\omega\in[\omega_{1},\omega_{2}),\\\
~{}~{}~{}~{}\vdots\\\
\mathcal{G}^{(p)}_{u}&~{}~{}\text{for}~{}\omega\in[\omega_{p-1},0.5).\end{array}\right.$
(2.8)

We want to estimate the number of segments in the partition $p$ and the
associated partition points $(\omega_{1},\omega_{2},\dots,\omega_{p})$ of the
frequency space.

##  3 Empirical Band Analysis

###  3.1 Estimation of the Power Spectrum and Proposed Statistic

We start by approximating the time-varying power spectrum by multitaper local
periodograms based on the data. We consider $B$ equally-sized non-overlapping
temporal blocks. For the sake of notational convenience, suppose that $T$ is a
multiple of $B$, and let $T_{B}=T/B$ be the number of observations in each
temporal block.

For $b=1,2,\dots,B$, $\omega\in(0,0.5)$, $k=1,2,\dots,K$ and $T\geq 1$, the
functional discrete Fourier transform (fDFT) is defined as a random function
with values in $L^{2}([0,1],\mathbb{C})$ given by

$\widetilde{X}_{T}^{(k),b,\omega}:=\sum_{t=1}^{T}v_{b}^{k}(t){I}_{b}(t/T)X_{t,T}e^{-\mathrm{i}2\pi\omega
t},$ (3.1)

where $v_{b}^{k}:\\{1,2,\dots,T\\}\mapsto\mathbb{R}$ is the $k$-th data taper
for the $b$-th time segment. We propose the use of sinusoidal tapers of the
form

$v_{b}^{k}(t)=\sqrt{\frac{2}{T_{B}+1}}\sin\frac{\pi
k[t-(b-1)T_{B}]}{T_{B}+1},\;\mathrm{for}\;k=1,\ldots,K,$ (3.2)

which are orthogonal for $t\in[(b-1)T_{B}+1,\ldots,bT_{B}]$. Sinusoidal tapers
are more computationally efficient than the Slepian tapers proposed in [35],
which require numerical eigenvalue decomposition to construct the tapers, and
can achieve similar spectral concentration with significantly less local bias
[30]. Another advantage is that the bandwidth, $bw$, which is the minimum
separation in frequency between approximately uncorrelated spectral estimates,
can be fully determined by setting the number of tapers, $K$, appropriately
[40]. More specifically,

$bw=\frac{K+1}{T_{B}+1}.$ (3.3)

The k-th single local periodogram kernel is then defined by

$\widehat{f}_{b,\omega}^{(k)}(\tau,\sigma)=\widetilde{X}_{T}^{(k),b,\omega}(\tau)\widetilde{X}_{T}^{(k),b,-\omega}(\sigma)=\widetilde{X}_{T}^{(k),b,\omega}(\tau)\overline{\widetilde{X}_{T}^{(k),b,\omega}(\sigma)}.$

The multitaper estimator of the local power spectrum within the $b$-th block
is then defined to be the average of all $K$ single taper estimators

$\displaystyle\widehat{f}_{b,\omega}^{(mt)}(\tau,\sigma)=\frac{1}{K}\sum_{k=1}^{K}\widehat{f}_{b,\omega}^{(k)}(\tau,\sigma).$
(3.4)

The final estimator of the time-varying power spectrum is then given by

$\displaystyle\widehat{f}_{u,\omega}^{(mt)}(\tau,\sigma)=\sum_{b=1}^{B}I_{b}(u)\widehat{f}_{b,\omega}^{(mt)}(\tau,\sigma).$
(3.5)

In practice, this estimator requires appropriate selection of two tuning
parameters: $B$ and $K$. The number of time blocks, $B$, balances frequency
and temporal properties. It should be selected small enough to ensure
sufficient frequency resolution and so the central limit theorem holds for the
local tapered periodgrams. It should be selected large enough so that data
within each block are approximately stationary, which depends on the signal
under study. Asymptotic results presented in the next section indicate the
optimal rate for $B$ is $T^{1/2}$. In the absence of scientific guidelines, we
recommend selecting $B$ as the factor of $T$ closest to $T^{1/2}$. The number
of tapers, $K$, controls the smoothness of local spectral estimates.
Asymptotic results presented in the next section indicate an appropriate
choice would be
$K=\max\left(1,\left\lfloor\min(T_{B},B)^{1/2}\right\rfloor\right)$.

We further define the estimated demeaned time-varying spectra as

$\widehat{g}_{b/B,\omega}(\tau,\sigma)=\widehat{f}_{b/B,\omega}^{(mt)}(\tau,\sigma)-\frac{1}{B}\sum_{l=1}^{B}\widehat{f}_{l/B,\omega}^{(mt)}(\tau,\sigma).$
(3.6)

Let $\widetilde{g}_{u,\omega_{0},\delta}$ be the average of demeaned time-
varying spectrum for $\omega\in[\omega_{0},\omega_{0}+\delta).$ In order to
identify frequency partition points, we consider an integrated scan statistic
defined as

$\displaystyle Q_{\omega_{0},\delta}$
$\displaystyle~{}=~{}\displaystyle\sum_{b=1}^{B}\int_{0}^{1}\int_{0}^{1}\left(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma)-\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma)\right)^{2}d\tau
d\sigma.$ (3.7)

###  3.2 Asymptotic Properties of the Scan Statistic

The basic intuition behind this scan statistic is that the estimator
$\widehat{g}$ defined in (3.6) asymptotically behaves like the theoretical
demeaned power spectrum defined in (2.7). The scan statistic compares the
value of $\widehat{g}$ at $\omega_{0}+\delta$ with the average of the same
quantity over the interval $(\omega_{0},\omega_{0}+\delta)$. Hence, if one of
the partition points $\omega_{k}$ defined in (2.7) is present in the interval
$(\omega_{0},\omega_{0}+\delta)$, the scan statistic should take a large
value, and it should be bounded and close to zero otherwise. For rest of the
paper, we will consider the following asymptotic scheme.

######  Assumption 3.1.

Assume $T\to\infty$, $B\to\infty$, $K\to\infty$, $T/K^{2}\to\infty$ and
$T/BK\to\infty,$

Our first result characterizes the asymptotic behavior of the estimated
demeaned power spectrum.

######  Lemma 3.1.

Under Assumption 2.1 and Assumption 3.1, for every
$b,b_{1},b_{2}\in\\{1,2,\dots,B\\},\omega\in(0,0.5)$, and
$\tau,\tau_{1},\tau_{2},\sigma,\sigma_{1},\sigma_{2}\in[0,1]$,

$\displaystyle\mathbb{E}\left(\widehat{g}_{b/B,\omega}(\tau,\sigma)\right)=g_{u_{b},\omega}(\tau,\sigma)+O\left(\log(T_{B})/T_{B}\right)+\log\left(1/T\right)$
$\displaystyle\text{Cov}\left(\widehat{g}_{b/B,\omega_{1}}(\tau_{1},\sigma_{1}),\widehat{g}_{b/B,\omega_{2}}(\tau_{2},\sigma_{2})\right)=\frac{(1-2/B)}{K}F(u_{b},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})$
$\displaystyle~{}~{}\hskip
162.6075pt+\frac{1}{B^{2}K}\sum_{l=1}^{B}F(u_{l},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})+o(1).$
$\displaystyle\text{Cov}\left(\widehat{g}_{b_{1}/B,\omega_{1}}(\tau_{1},\sigma_{1}),\widehat{g}_{b_{2}/B,\omega_{2}}(\tau_{2},\sigma_{2})\right)=-\frac{1}{BK}F(u_{b_{1}},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{BK}F(u_{b_{2}},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})+\frac{1}{B^{2}K}\sum_{l=1}^{B}F(u_{b_{l}},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})+o(1),$

where

$F(u,\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}):=f_{u,\omega_{1}}(\tau_{1},\tau_{2})f_{u,\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u,\omega_{2}}(\tau_{1},\sigma_{2})f_{u,\omega_{2}}(\tau_{2},\sigma_{1}),$

and $u_{b}$ is the mid-point of the $b$-th block.

######  Remark 3.1.

Note that the assumptions $T/BK\to\infty$ and $K\to\infty$ together guarantee
$T_{B}\to\infty$. Therefore Lemma 3.1 implies that
$\mathbb{E}\left(\widehat{g}_{b/B,\omega}(\tau,\sigma)\right)\to
g_{u_{b},\omega}(\tau,\sigma)$ under our asymptotic scheme. Moreover, the
covariance kernel of the process $\widehat{g}$ is $O(1/K)$. This guarantees
for fixed $b$ and $\omega$, $\widehat{g}_{b/B,\omega}$ consistently estimates
$g_{u_{b},\omega}$. A stronger uniform result on $b$ and $\omega$ can be
established under additional moment assumptions. Additionally, the covariance
across different time blocks is of the order $O(1/BK)$, which converges to $0$
faster than the covariance kernel within the same block.

Next, we investigate the asymptotic behavior of the scan statistic itself. The
next two theorems describe the asymptotic behavior of the scan statistic under
the absence and presence of a jump point in the interval
$(\omega_{0},\omega_{0}+\delta)$ respectively.

######  Theorem 3.1.

Assume $g_{u,\omega}(\tau,\sigma)=g^{(0)}_{u}(\tau,\sigma)$ for all
$\tau,\sigma\in[0,1]$ and $\omega\in[\omega_{0},\omega_{0}+\delta].$ Under
Assumption 2.1 and Assumption 3.1,

$Q_{\omega_{0},\delta}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{K}\sum_{b=1}^{B}\left(\|\mathcal{G}_{b}\|_{2}^{2}+o_{p}(1)\right),$

where $\\{\mathcal{G}_{b}\\}$ is a collection of zero-mean Gaussian processes
in $L^{2}\left([0,1]^{2}\right)$ with covariance structure given in (.3) and
(.3).

By Lemma .4, the quantity
$\frac{1}{B}\sum_{b}\|\mathcal{G}_{b}\|^{2}=O_{p}(1)$ as $B\to\infty$.
Therefore, Theorem 3.1 suggests that the quantity $K/BQ_{\omega_{0},\delta}$
is asymptotically bounded. Hence, as long as the last quantity diverges under
our asymptotic set-up, we can construct a consistent test for the existence of
a partition point in frequency domain. The next theorem guarantees that is
indeed the case.

######  Theorem 3.2.

Assume that

$\displaystyle
g_{u,\omega}(\tau,\sigma)=\left\\{\begin{array}[]{ll}g^{(1)}_{u}(\tau,\sigma)&\text{for
}\omega\in[\omega_{0},\omega^{*})\\\ g^{(2)}_{u}(\tau,\sigma)&\text{for
}\omega\in[\omega^{*},\omega_{0}+\delta].\end{array}\right.$

Under Assumption 2.1 and Assumption 3.1,

$Q_{\omega_{0},\delta}\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{K}\sum_{b=1}^{B}\left(\|\mathcal{G}_{b}\|_{2}^{2}+o_{p}(1)\right)+B\left(\frac{\omega^{*}-\omega_{0}}{\delta}\right)^{2}\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\left(g^{(1)}_{u}(\tau,\sigma)-g^{(2)}_{u}(\tau,\sigma)\right)^{2}d\tau
d\sigma du+O_{p}(B/K),$

where $\\{\mathcal{G}_{b}\\}$ are as defined in Theorem 3.1.

######  Remark 3.2.

If there is no significant frequency partition point in
$[\omega_{0},\omega_{0}+\delta]$, we have

$Q_{\omega_{0},\delta}=O(B/K),$

and in the presence of a change,

$Q_{\omega_{0},\delta}=O(B/K)+O(B).$

Therefore, we expect to see a large spike in the scan statistic around the
frequency where the dynamics of the power spectrum changes.

######  Remark 3.3.

Suppose we are interested in testing the null hypothesis that the kernel
$g_{u,\omega}$ does not change in $\omega$ on the interval
$(\omega_{0},\omega_{0}+\delta)$. Let $c_{\alpha}$ be the $(1-\alpha)$-th
quantile of $\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|_{2}^{2}$ for large
$B$. Then an asymptotic level $\alpha$ test for this hypothesis can be
constructed by rejecting the null when $KQ_{\omega_{0},\delta}>Bc_{\alpha}$.
In fact, the next Corollary suggests that this is a consistent test.

######  Corollary 3.1.

Under the assumptions of Theorem 3.2, we have

$P\left(Q_{\omega_{0},\delta}>\frac{1}{K}\sum_{b=1}^{B}\VERT\mathcal{G}_{b}\VERT_{2}^{2}\right)\to
1.$

###  3.3 An Iterative Algorithm

The results discussed in the previous section provide a consistent test to
find a single partition point in the frequency space. In this section, we
extend this framework to detect multiple frequency partition points using an
iterative search algorithm that uses this scan statistic to efficiently
explore the frequency space and identify all frequency partition points.

For computational efficiency, we search across frequencies in batches of size
$n_{\max}$. This limits the number of calculations required to approximate the
null distribution of the scan statistic to avoid undue computational burden.
In particular, $\omega_{0}$ is first fixed at a value near $0$, and we conduct
a test for partition points in the interval $(\omega_{0},\omega_{0}+\delta]$
comparing the statistic $Q_{\omega_{0},\delta}$ against the null distribution
described in Theorem 3.1 for different values of $\delta<n_{\max}/T_{B}$. A
Hochberg step-up procedure [15] is then used to test for the presence of a
partition in the frequency domain at each particular frequency
$\omega_{0}+\delta$. If the procedure returns a set of frequencies for which
the null hypothesis (no partition point) is rejected, we then add the smallest
frequency in this set to the estimated frequency band partition, increase
$\omega_{0}$ to be just larger than this newly found frequency partition
point, and repeat the process. If the procedure does not return any
frequencies for which the null hypothesis is rejected, we increase
$\omega_{0}$ to be just larger than the largest frequency tested in the
current batch and repeat the process. The procedure continues until all
frequencies have been evaluated as potential partition points. To better
visualize this procedure, consider that the traversal of this procedure across
frequencies resembles the movement of an inchworm. A complete algorithmic
representation of the search procedure in its entirety is available in
Algorithm 1.

P-values are obtained by comparing observed test statistics with a simulated
distribution of the limiting random variable described in Theorem 3.1. To
simulate from the null distribution, we generate $d_{0}$ draws from the
collection of zero-mean Gaussian processes $\\{\mathcal{G}_{b}\\}$ with
covariance structure given in (.3) and (.3). The covariance is estimated using
the estimates of the demeaned time-varying power spectrum based on the local
multitaper power spectrum estimator $\hat{f}^{(mt)}_{u,\omega}$ described in
Section 3.1.

Input : Demeaned time-varying multitaper power spectrum estimates,
$\hat{g}_{b/B,\omega_{k}}(\tau,\sigma)$, for $b=1,\ldots,B$ and
$\omega_{k}=k/T_{B}$ for $k=1,\ldots,N_{B}=\lfloor T_{B}/2\rfloor-1$
Number of tapers $K$, significance level $\alpha$, number of frequencies
tested in each pass $n_{\max}$, number of draws for approximating $p$-values,
$d_{0}$

Output : Estimated number of partition points, $\widehat{p}$
Estimated partition points,
$\widehat{\boldsymbol{\omega}}_{\widehat{p}}=\\{\widehat{\omega}_{1},\widehat{\omega}_{2},\ldots,\widehat{\omega}_{\widehat{p}-1}\\}$

$\widehat{p}\leftarrow 1$,
$\widehat{\boldsymbol{\omega}}_{\widehat{p}}\leftarrow\\{\\}$,
$\mathrm{stop}\leftarrow 0$, $\epsilon\leftarrow\frac{K+1}{T_{B}+1}$,
$k^{*}\leftarrow\left\lceil T_{B}\epsilon\right\rceil$,
$\omega_{0}\leftarrow\omega_{k^{*}}$

while  _$\mathrm{stop}\neq 1$_ do

$k_{\min}\leftarrow\left\lceil T_{B}\epsilon\right\rceil$,
$k_{\max}\leftarrow\min\left(k_{\min}+n_{\max}-1,\left\lfloor
N_{B}-k^{*}-T_{B}\epsilon\right\rfloor\right)$

Compute test statistics
$Q_{\omega_{0},\delta_{k}}\;\forall\;\delta_{k}=k/T_{B},k\in\\{k_{\min},k_{\min}+1,\ldots,k_{\max}\\}.$

Simulate $d_{0}$ draws,
$\mathbf{Q}^{H_{0}}_{\omega_{0},\delta_{k}}=\left\\{Q^{H_{0},i}_{\omega_{0},\delta_{k}}\right\\}_{i=1}^{d_{0}}$,
from the limiting null distribution of $Q_{\omega_{0},\delta_{k}}$.

Approximate $p$-values,
$\hat{p}(k)=\frac{1}{d_{0}}\sum_{i=1}^{d_{0}}I\left(Q^{H_{0},i}_{\omega_{0},\delta_{k}}>Q_{\omega_{0},\delta_{k}}\right)$
to test $H_{0}(k):g_{u,\omega}(\tau,\sigma)=g_{u}^{(0)}(\tau,\sigma)$ for
$\omega\in[\omega_{0},\omega_{0}+\delta_{k}]$.

Identify
$\mathbf{R_{\alpha}}=\\{\omega_{0}+\delta_{k}:H_{0}(k)\;\mathrm{rejected}\\}$
using level-$\alpha$ Hochberg step-up procedure.

if  _$\mathbf{R_{\alpha}}=\\{\\}$ _ then

$k^{*}\leftarrow T_{B}(\omega_{0}+\delta_{k_{\max}}),\;\omega_{0}\leftarrow
k^{*}/N_{B}$

else

$\hat{\omega}^{*}\leftarrow\min\mathbf{R_{\alpha}},\;\widehat{\boldsymbol{\omega}}_{\widehat{p}}\leftarrow\widehat{\boldsymbol{\omega}}_{\widehat{p}}\cup\\{\hat{\omega}^{*}\\},\;\widehat{p}\leftarrow\widehat{p}+1\;$

$k^{*}\leftarrow\left\lceil
T_{B}(\hat{\omega}^{*}+\epsilon)\right\rceil,\;\omega_{0}\leftarrow
k^{*}/N_{B}$

end if

if  _$\omega_{0} >\omega_{N_{B}}-2\epsilon$_ then stop $\leftarrow 1$

end while

return _$\widehat{p},\widehat{\boldsymbol{\omega}}_{\widehat{p}}$_

Algorithm 1 Inchworm Frequency Band Search Algorithm

######  Remark 3.4.

Computation cost is a critical aspect for analyzing functional time series
data. In order to make our algorithm more efficient, we use two particular
adjustments

  1. (i)

A block diagonal approximation of the asymptotic covariance structure: While
approximating the asymptotic quantiles by simulating $\\{\mathcal{G}_{b}\\}$,
we ignore the covariance across different time blocks $b$ given by (.3). This
allows for simulating from a collection of $B$ independent, lower dimensional
Gaussian processes, which can be carried out in parallel for computational
efficiency. Note that the covariance kernel for $\mathcal{G}_{b}$ given in
(.3) is $O(1)$ and the covariance kernel in (.3) is $O(1/B)$. Therefore, for
large $B$, the approximation works reasonably well.

  2. (ii)

Choice of the tuning parameter $n_{\max}$: Notice that number of calculations
required to estimate the covariance kernel given in (.3) grows with the number
of frequencies between $\omega_{0}$ and $\omega_{0}+\delta$. By traversing
frequencies by testing in smaller batches, this reduces computation times
significantly compared to letting $\delta$ vary over the whole frequency
domain in a single pass.

All computations in what follows were performed using R 4.0.3 [25]. With the
computational improvements introduced above, the run time using the proposed
search algorithm for a single realization of functional white noise of length
$T=2000$ with $R=5$ observed points in the functional domain using $B=5$ time
blocks and $K=15$ tapers with $n_{\max}=40$ and $d_{0}=100000$ is 5 to 6
minutes on a laptop with 32GB RAM and a 8-core 2.4 GHz processor on a Mac
operating system. Without these improvements, the run time for the same data
settings is over 10 hours, and many of the larger data settings considered in
Section 4.1 are not computationally feasible.

###  3.4 Test for Stationarity Within Frequency Blocks

The proposed methodology produces homogeneous regions of frequencies in which
the power spectrum varies only across time. It is natural to then seek to
identify frequency bands for which the second order structure is stationary.
Specifically, within a frequency band $[\omega_{1},\omega_{2}]$, the power
spectrum $f_{u,\omega}(\tau,\sigma)=f_{u}(\tau,\sigma)$ for all
$u,\tau,\sigma\in[0,1]$ and $\omega\in[\omega_{1},\omega_{2}].$ Furthermore,
if the process in stationary within this region, the power spectrum is
constant across time $u$, i.e., $f_{u,\omega}(\tau,\sigma)=f(\tau,\sigma)$
within than band. In this situation, the demeaned power spectrum

$g_{u,\omega}(\tau,\sigma)=f_{u,\omega}(\tau,\sigma)-\int_{0}^{1}f_{u,\omega}(\tau,\sigma)du\equiv
0,$

for for all $u,\tau,\sigma\in[0,1]$ and $\omega\in[\omega_{1},\omega_{2}].$
Therefore, to test stationarity within a frequency band, we consider null
hypothesis

$H_{0}:g_{u,\omega}(\tau,\sigma)\equiv 0,~{}~{}\text{ almost everywhere
}u,\tau,\sigma\in[0,1],~{}\omega\in[\omega_{1},\omega_{2}],$

against the alternative

$H_{a}:g_{u,\omega}(\tau,\sigma)\neq 0,~{}~{}\text{on a set of positive
Lebesgue measure.}$

To test this hypothesis, we propose the test statistic

$Q^{0}(\omega_{1},\omega_{2})=\frac{1}{B}\sum_{b=1}^{B}\int_{\omega_{1}}^{\omega_{2}}\int_{0}^{1}\int_{0}^{1}\left|\widehat{g}_{b/B,\omega}(\tau,\sigma)\right|^{2}d\tau
d\sigma d\omega.$

As $\widehat{g}$ is a consistent estimator of the true demeaned power spectrum
$g$, this statistic is close to $0$ under $H_{0}$ and takes large positive
values under the alternative. Now we formalize this idea through the
asymptotic distribution of this statistic.

######  Theorem 3.3.

Assume that $g_{u,\omega}(\tau,\sigma)=g_{u}(\tau,\sigma)=0$ for all
$u,\tau,\sigma\in[0,1]$ and $\omega\in[\omega_{1},\omega_{2}].$ Under
Assumption 2.1 and Assumption 3.1,

$Q^{0}(\omega_{1},\omega_{2})\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{K}\sum_{b=1}^{B}\left(\|\mathcal{H}_{b}\|_{2}^{2}+o(1)\right),$
(3.8)

where $\\{\mathcal{H}_{b}\\}$ is a collection of zero-mean Gaussian Process in
$L^{2}[0,1]^{2}$ with covariance structure given in (A.7).

######  Remark 3.5.

Theorem 3.3 suggests that under $H_{0}$, $g_{u}(\tau,\sigma)\equiv 0$, the
quantity $K/B\times
Q^{0}(\omega_{1},\omega_{2})=\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{H}_{b}\|^{2}+o(1)$.
Therefore, approximate $p$-values for the test of stationarity can be
constructed by comparing $KQ^{0}(\omega_{1},\omega_{2})/B$ with simulated
quantiles of $\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{H}_{b}\|^{2},$ for large
$B$. Note that $\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{H}_{b}\|^{2}$ is
$O_{p}(1)$ as $B\to\infty$ by Lemma .4.

The following Lemma guarantees consistency of the proposed test.

######  Lemma 3.2.

Assume that $g_{u,\omega}(\tau,\sigma)=g_{u}(\tau,\sigma)\neq 0$ on a set of
positive Lebesgue measure. Under Assumption 2.1 and Assumption 3.1,
$KQ^{0}(\omega_{1},\omega_{2})/B\to\infty$ in probability.

##  4 Finite Sample Properties

###  4.1 Simulation Studies

In order to evaluate the performance of the search algorithm in finite
samples, we consider three simulation settings representing appropriate
extensions of [7] for functional time series. A B-spline basis with 15 basis
functions is used to generate random realizations of the functional time
series. All settings can be represented as
$f_{u,\omega}(\tau,\sigma)=\phi_{u,\omega}f(\tau,\sigma)$ for $u\in[0,1]$ and
$\omega\in(0,0.5)$ and

$\phi_{u,\omega}=1\;\mathrm{for}\;\omega\in(0,0.5),$ (4.1)
$\phi_{u,\omega}=\begin{cases}10-9u\;\mathrm{for}\;\omega\in(0,0.15)\\\
5\;\mathrm{for}\;\omega\in[0.15,0.35)\\\
1+9u\;\mathrm{for}\;\omega\in[0.35,0.5),\\\ \end{cases}$ (4.2)

and

$\phi_{u,\omega}=\begin{cases}2+\sin(8\pi
u-\pi/2)\;\mathrm{for}\;\omega\in(0,0.15]\\\ 2+\cos(8\pi
u)\;\mathrm{for}\;\omega\in(0.15,0.35]\\\ 2+\cos(16\pi
u)\;\mathrm{for}\;\omega\in(0.35,0.5).\\\ \end{cases}$ (4.3)

See Figure 1 for an illustration of the estimated time-varying auto spectrum
for a single point in the functional domain for each setting. In the first
setting, we consider functional white noise in order to ensure the method
maintains appropriate control for false positives. In the second and third
settings, both linear and non-linear nonstationary dynamics are considered
within frequency bands. These settings are specially designed to assess
performance in detecting time-varying dynamics of different forms, as well as
subtle changes in dynamics over frequencies.

Figure 1: Local multitaper estimates of the time-varying autospectra for a
single component from each of the simulation settings for time series with
$N=$1000 observations per time segment and $B=$20 time segments. Solid green
lines represent the true partition of frequencies.

Table 1 reports the means and standard deviations over 100 replications for
the estimated number of frequency bands, $\hat{p}$, and Rand indices,
$R(\hat{\boldsymbol{\omega}},\boldsymbol{\omega})$. The Rand index [28]
summarizes the similarity between the estimated frequency partition,
$\hat{\boldsymbol{\omega}}$, and the true partition, $\boldsymbol{\omega}$.
Let $a$ be the number of pairs of Fourier frequencies in the same frequency
band in $\hat{\boldsymbol{\omega}}$ and the same frequency band in
$\boldsymbol{\omega}$ and $b$ be the number of pairs of Fourier frequencies in
different frequency bands in $\hat{\boldsymbol{\omega}}$ and different
frequency bands in $\boldsymbol{\omega}$. Then the Rand index,
$R(\hat{\boldsymbol{\omega}},\boldsymbol{\omega})=(a+b)/\binom{N_{B}}{2}$,
where $N_{B}=\lfloor T_{B}/2\rfloor-1$ is the number of Fourier frequencies.
The Rand index can take on values from 0 to 1 with values close to 1
indicating good estimation of the true frequency partition. For these results,
we consider performance under different combinations of settings for the
number of observations in each time block, $T_{B}$, the number of
approximately stationary time blocks, $B$, and the number of observations in
the functional domain, $R$. Furthermore, we fix the family-wise error rate
(FWER) control level, $\alpha=0.05$, number of frequencies tested in each
pass, $n_{\max}=30$, number of draws to approximate the null distribution of
the test statistics $d_{0}=100,000$, local multitaper estimator bandwidth,
$bw=0.05$, such that the number of tapers, $K=\left\lfloor
bw(N+1)\right\rfloor-1$, and we use a block diagonal approximation for the
covariance function of the Gaussian process used to approximate the limiting
null distribution of the scan statistics and select four approximately
equally-spaced points in the functional domain for testing. The proposed
method provides good estimation accuracy for both the number and location of
frequency band partition points, and performance generally improves as $T_{B}$
and $B$ increase.

|  | R=5 | R=10
---|---|---|---
$T_{B}$ | B | $\hat{p}$ | $R(\hat{\boldsymbol{\omega}},\boldsymbol{\omega})$ | $\hat{p}$ | $R(\hat{\boldsymbol{\omega}},\boldsymbol{\omega})$
White noise ($p=1$)
200 | 5 | 1.000(0.000) | 1.000(0.000) | 1.020(0.141) | 0.996(0.030)
10 | 1.000(0.000) | 1.000(0.000) | 1.000(0.000) | 1.000(0.000)
20 | 1.000(0.000) | 1.000(0.000) | 1.000(0.000) | 1.000(0.000)
500 | 5 | 1.970(0.810) | 0.655(0.248) | 1.850(0.796) | 0.701(0.250)
10 | 1.350(0.539) | 0.862(0.209) | 1.260(0.525) | 0.900(0.195)
20 | 1.130(0.367) | 0.944(0.155) | 1.130(0.367) | 0.948(0.146)
1000 | 5 | 1.800(0.841) | 0.735(0.255) | 1.780(0.786) | 0.736(0.251)
10 | 1.060(0.239) | 0.974(0.105) | 1.150(0.411) | 0.941(0.157)
20 | 1.020(0.141) | 0.992(0.059) | 1.010(0.100) | 0.996(0.042)
Linear ($p=3$)
200 | 5 | 2.120(0.498) | 0.738(0.129) | 2.210(0.498) | 0.762(0.115)
10 | 2.200(0.471) | 0.761(0.102) | 2.240(0.495) | 0.767(0.105)
20 | 2.460(0.501) | 0.816(0.083) | 2.530(0.502) | 0.826(0.082)
500 | 5 | 3.260(0.463) | 0.936(0.050) | 3.220(0.504) | 0.923(0.066)
10 | 3.050(0.261) | 0.955(0.041) | 3.050(0.297) | 0.953(0.049)
20 | 3.030(0.171) | 0.963(0.036) | 3.010(0.100) | 0.967(0.033)
1000 | 5 | 3.450(0.609) | 0.915(0.033) | 3.370(0.580) | 0.921(0.035)
10 | 3.080(0.273) | 0.933(0.023) | 3.030(0.171) | 0.934(0.021)
20 | 3.020(0.141) | 0.931(0.012) | 3.020(0.141) | 0.931(0.014)
Sinusoidal ($p=3$)
200 | 5 | 1.060(0.239) | 0.354(0.088) | 1.080(0.273) | 0.359(0.097)
10 | 2.000(0.201) | 0.732(0.062) | 1.980(0.200) | 0.727(0.072)
20 | 2.030(0.171) | 0.753(0.026) | 2.060(0.239) | 0.754(0.028)
500 | 5 | 2.600(0.865) | 0.726(0.143) | 2.590(0.805) | 0.718(0.153)
10 | 2.950(0.458) | 0.897(0.085) | 2.910(0.534) | 0.893(0.088)
20 | 3.040(0.243) | 0.952(0.039) | 2.990(0.266) | 0.945(0.054)
1000 | 5 | 2.300(0.916) | 0.651(0.182) | 2.140(0.954) | 0.620(0.192)
10 | 2.660(0.623) | 0.829(0.106) | 2.720(0.533) | 0.854(0.102)
20 | 3.040(0.197) | 0.932(0.016) | 3.010(0.100) | 0.934(0.012)

Table 1: Mean(standard deviation) for the estimated number of frequency bands,
$\hat{p}$, and Rand index values,
$R(\boldsymbol{\hat{\omega}},\boldsymbol{\omega})$, for $R=100$ replications.

When the number of time blocks, $B$, is smaller, the proposed method slightly
overestimates the number of partition points. This is not unexpected, since
the proposed method uses the asymptotic behavior of the scan statistics for
testing purposes, which may not hold for smaller values of $B$. For the first
setting, Table 1 indicates good performance in correctly estimating only one
frequency band. Since the search algorithm is designed to control FWER, the
false positive rate remains under control, even as the number of distinct
frequencies tested increases with larger values of $T_{B}$. In some cases,
FWER may be controlled more tightly than $\alpha=0.05$ due to dependence among
significance tests across frequencies [12, 34, 8]. For the second setting and
third settings, performance generally improves as $T_{B}$ and $B$ increase.
However, the accuracy is also impacted by the magnitude of the differences
between the underlying demeaned time-varying power spectra across frequency
bands, as well as the number of time blocks used to approximate the time-
varying behavior. The time-varying dynamics of the power spectrum for adjacent
frequency bands for the second setting are dissimilar for nearly all time
points. Accordingly, the algorithm is able to correctly estimate the frequency
band partition with smaller $T_{B}$ and $B$ for this setting. On the other
hand, the two frequency bands covering higher frequencies in the third setting
have similar time-varying dynamics in the power spectrum at particular time
points due to their periodicities (see Figure 1). Also, without a sufficient
number of time blocks, it is difficult to distinguish the different periodic
time-varying behavior across frequency bands for this setting. Taken together,
this explains the need for relatively larger values of $T_{B}$ and $B$ for
accurate estimation of the frequency band partition for the third setting
compared to the second setting.

###  4.2 Frequency Band Analysis for EEG Data

To illustrate the usefulness of the proposed methodology for functional time
series analysis, we turn to frequency band analysis of EEG signals. Analyzing
EEG signals as nonstationary functional time series is warranted by the high-
dimensionality, nonstationarity, and strong dependence typically observed for
EEG signals. Frequency bands are also commonly used in the scientific
literature to generate summary measures of EEG power spectra, so a principled
approach to frequency band estimation would be a welcomed development. In the
scientific literature, there is significant variability in frequency ranges
used to define traditional EEG frequency bands [22]. For the purposes of
comparison with the proposed method, we use the following definitions (Hz)
[21]: delta $(0,4)$, theta $[4,7)$, alpha $[7,12)$, beta $[12,30)$, and gamma
$[30,100)$. We analyze a 4-minute segment of a 72-channel 256 Hz EEG signal
from a single participant from the study described in [36]. Participants in
the study sat in a wakeful resting state and alternated between eyes open (EO)
and eyes closed (EC) conditions at 1 minute intervals. Given this particular
design, we expect to see periodic nonstationary behavior similar to that of
the third simulation setting introduced in Section 4.1.

For computational efficiency, we downsample the time series at 32 Hz to
produce a series of length $T=7680$. For the local multitaper estimator of the
time-varying power spectrum, we use 10-second segments resulting in
$T_{B}=32\times 10=320$ observations per segment and $B=7680/320=24$ segments.
Other parameters (e.g. FWER control level, number of frequencies tested in
each pass, local multitaper bandwidth, etc.) follow the same settings used for
the simulations in Section 4.1. Since it is possible that different brain
regions may be characterized by different frequency band structures, we
consider two groups of channels, one representing the parietal and occipital
lobes (PO7, PO3, POz, PO4, and PO8), which is associated with attention, and
one representing the frontal and central lobes (F1, Fz, F2, FC1, FCz, FC2),
which is associated with sensorimotor function [41]. By applying the proposed
method to each group, we can better understand if the frequency band
structures that characterize the time-varying dynamics of the power spectrum
are different within sub-regions of the functional domain.

Figure 2 presents the log autospectra for the five channels associated with
attentional system areas along with traditional EEG frequency bands and
estimated EEG frequency bands using the proposed method. Applying the proposed
methodology to this data revealed three frequency bands with different time-
varying dynamics (Hz): (0,3.8), [3.8,12.6), [12.6,16). Comparing with the
traditional frequency band partition, the proposed low frequency band
$(0,3.8)$ coincides with the traditional delta band $(0,4)$, but the next
proposed frequency band $[3.8,12.6)$ covers both the traditional theta $[4,7)$
and alpha $[7,12)$ frequency bands. This suggests that while the delta band
exhibits different time-varying characteristics from other bands, the theta
and alpha bands exhibit similar time-varying characteristics for these
channels. These results are not surprising, as it is well-known that alpha
band power increases during EC conditions and attenuates with visual
stimulation during EO conditions. As is the case with this participant,
similar behavior has also been observed for theta band power, with larger
reductions observed in the posterior regions of the brain, including the
attentional system area channels currently under study [2, 3]. Hence, it is
reasonable that the time-varying dynamics of the alpha and theta bands may be
similar for this particular set of EEG channels.

Figure 2: Log time-varying autospectra for 5 EEG channels and frequency bands
determined by the proposed methodology (solid green lines) and traditional
frequency bands (dashed blue lines) are displayed. These EEG channels measure
activity in the parietal (P) and occipital (O) lobes, which are associated
with attention and visual processing.

To better understand differences in the time-varying behavior across the
estimated bands, Figure 3 displays a smoothed estimate of the frequency-band
specific demeaned time-varying autospectra,
$\hat{g}_{u}^{(p)}(\tau_{i},\tau_{i}),p=1,2,3$ for each channel
$i=1,\ldots,5$. Smoothing was performed using cubic splines with 4 knots to
better visualize the slowly-evolving time-varying dynamics under the
assumption of local stationarity. It can be seen that the time-varying
dynamics for the estimated frequency band corresponding to the traditional
delta band, (0, 3.8), coincides with the time-varying dynamics of the
estimated frequency band covering the theta and alpha bands, [3.8, 12.6), for
some channels (PO7, PO4, PO8), but not others (PO3, POz). Since the search
algorithm relies on an integrated scan statistic that integrates over the
functional domain, it is sensitive to differences in the time-varying dynamics
for proper subsets of the functional domain. Accordingly, the proposed method
correctly distinguishes between lower frequencies (0,3.8) and higher
frequencies in the estimated frequency band structure.

Figure 3: Smoothed estimator of the frequency band specific demeaned time-
varying power spectra $g^{(p)}_{u}$ $p=1,2,3$ for 5 EEG channels using cubic
splines with 4 knots. Frequency bands are estimated used the proposed method
and these channels measure activity in the parietal (P) and occipital (O)
lobes.

Turning to the group of six EEG channels covering the frontal and central
lobes, Figure 4 presents the log autospectra along with traditional EEG
frequency bands and estimated EEG frequency bands using the proposed method.
Applying the proposed methodology to this data revealed five frequency bands
with different time-varying dynamics (Hz): (0,2.2), [2.2,4.9), [4.9, 8.1),
[8.1,11.8), [11.8,16.0). These estimated bands align reasonably well with the
traditional EEG frequency bands. However, the estimated bands suggest that the
traditional delta band, (0,4), should be characterized by two sub-bands,
(0,2.2) and [2.2, 4.9), which exhibit significantly different time-varying
behavior of the power spectrum. Such findings have been noted in the
scientific literature, in which the so-called “slow delta” (0.7-2 Hz) and
“fast delta” (2-4 Hz) bands exhibit different behavior during the wake-sleep
transition [4]. [2] also observed that the magnitude of the difference in
theta power between the EO and EC conditions is less for frontal and central
brain regions compared to posterior regions, while the magnitude of the
difference for the alpha band is similar across regions. This is supported by
the current analysis and can help explain why the theta and alpha bands are
estimated to have similar time-varying behavior for the group of posterior
region EEG channels (Figure 2), and different time-varying dynamics for the
group of central and frontal region EEG channels (Figure 4).

Figure 4: Log time-varying autospectra for 6 EEG channels and frequency bands
determined by the proposed methodology (solid green lines) and traditional
frequency bands (dashed blue lines) are displayed. These EEG channels measure
activity in the frontal (F) and central (C) lobes, which are associated with
sensorimotor function.

The smoothed estimates of the frequency-band specific demeaned time-varying
autospectra for these channels (see Figure 5) illustrate the different time-
varying behavior of the estimated frequency bands captured by the search
algorithm. The two sub-bands covering the traditional delta band indicate very
different time-varying behavior. Also, the estimated band that roughly
corresponds to the traditional alpha band, [8.1, 11.8), has a more regular and
pronounced time-varying behavior, corresponding to the alternating EO and EC
conditions, compared to the other estimated frequency bands. In summary, the
proposed search algorithm estimates frequency bands that can better
characterize the power spectrum for the particular functional EEG time series
under study, compared to traditional EEG frequency band analysis, and can be
used to construct customized frequency band summary measures for
characterizing different brain regions.

Figure 5: Smoothed estimator of the frequency band specific demeaned time-
varying power spectra $g^{(p)}_{u}$ $p=1,2,3$ for 5 EEG channels using cubic
splines with 4 knots. Frequency bands are estimated used the proposed method
and these channels measure activity in the frontal (F) and central (C) lobes.

##  5 Discussion

The frequency band analysis framework for nonstationary functional time series
introduced in this article offers a quantitative approach to identifying
frequency bands that best preserve the nonstationary dynamics of the
underlying functional time series. This framework allows for estimation of
both the number of frequency bands and their corresponding endpoints through
the use of a sensible integrated scan statistic within an iterative search
algorithm. Another test statistic is also offered to determine which bands, if
any, are stationary with respect to time. Motivated by the application to EEG
frequency band analysis, it would be interesting to consider extensions of
this framework that enable localization of the frequency band estimation
framework in the time and functional domains. Such extensions would allow for
the frequency band estimation framework to automatically adapt to local
spectral characteristics without needing to pre-specify particular time
segments or subsets of the functional domain for analysis. However, these
extensions present significant computational challenges associated with
searching over multiple spaces simultaneously.

We have focused on estimation of frequency bands for a single nonstationary
functional time series, but this framework can also be extended for the
analysis of multiple functional time series. For example, extending this
framework for estimating frequency bands for classification and clustering of
functional time series would provide researchers with optimal frequency band
features for supervised and unsupervised learning tasks. This extension could
be very useful in the study of EEG and fMRI signals to construct frequency
band features that are associated with clinical and behavioral outcomes or
that can be used to identify groups of time series with similar spectral
characteristics.

## Technical Details

###  .1 Properties of Multitaper Periodogram of Functional Time Series

We first establish some asymptotic properties of the multitaper periodogram
estimator defined in (3.5). Following the notations in [10] we define the
quantities

$H_{k}(\omega)=\sum_{t=1}^{T_{B}}e^{-2\pi\omega t}v_{b}^{k}(t),$
$H_{k,l}(\omega)=\sum_{t=1}^{T_{B}}e^{-2\pi\omega t}v_{b}^{k}(t)v_{b}^{l}(t),$

Note that,

$\int_{-1}^{1}|H_{k}(\omega)|^{2}d\omega=\sum_{t=1}^{T_{B}}(v_{b}^{k}(t))^{2}=1.$

Moreover by [10] we can write an upper bound $L(\omega)$ for both the
functions $\sqrt{T_{B}}H_{k}(\omega)$ and ${T_{B}}H_{k,k}(\omega)$ where

$L(\omega)=\left\\{\begin{array}[]{ll}{T_{B}},&\text{if }|\omega|\leq
1/T_{B}\\\ \frac{1}{|\omega|},&\text{otherwise}.\end{array}\right.$

######  Lemma .1.

Let $\tilde{X}_{T}^{(k),b,\omega}$ be the $k$-th functional discrete Fourier
transformation of the observed time series at $b$-th block and frequency
$\omega$, as defined in (3.1). Under Assumption 2.1 and Assumption 3.1, the
$l$-th order cumulant kernel of the fDFT is given by

$\displaystyle
cum\left(\tilde{X}_{T}^{(k_{1}),b_{1},\omega_{1}}(\tau_{1}),\dots,\tilde{X}_{T}^{(k_{l}),b_{l},\omega_{l}}(\tau_{l})\right)$
$\displaystyle=\left\\{\begin{array}[]{ll}H_{k_{1},\dots,k_{l}}\left(\sum_{j=1}^{l}\omega_{j}\right)f_{u_{b},\omega_{1},\dots,\omega_{l-1}}(\tau_{1},\dots,\tau_{l})+o(1),&\text{if
}b_{1}=b_{2}=\dots=b_{l}=b\\\ 0,&\text{otherwise. }\end{array}\right.$

where
$H_{k_{1},\dots,k_{l}}(\omega)=\displaystyle\sum_{t=1}^{T_{B}}v_{b}^{k_{1}}(t)v_{b}^{k_{2}}(t)\dots
v_{b}^{k_{l}}(t)e^{-i\omega t}.$

###### Proof.

Let $B_{b}$ be the $b$-th time block and write

$\displaystyle
cum\left(\tilde{X}_{T}^{(k_{1}),b_{1},\omega_{1}}(\tau_{1}),\dots,\tilde{X}_{T}^{(k_{l}),b_{l},\omega_{l}}(\tau_{l})\right)$
$\displaystyle=$ $\displaystyle cum\left(\sum_{t\in
B_{b_{1}}}v_{b_{1}}^{k_{1}}(t)X_{t,T}(\tau_{1})e^{-\mathrm{i}2\pi\omega_{1}t},\dots,\sum_{t\in
B_{b_{l}}}v_{b_{l}}^{k_{l}}(t)X_{t,T}(\tau_{l})e^{-\mathrm{i}2\pi\omega_{l}t}\right)$
$\displaystyle=$ $\displaystyle\sum_{t_{1}\in B_{b_{1}}}\dots\sum_{t_{l}\in
B_{b_{l}}}v_{b_{1}}^{k_{1}}(t_{1})\dots
v_{b_{l}}^{k_{l}}(t_{l})e^{-\mathrm{i}2\pi\sum_{k=1}^{l}\omega_{k}t_{k}}cum\left(X_{t_{1},T}(\tau_{1}),\dots,X_{t_{l},T}(\tau_{l})\right)$
$\displaystyle=$ $\displaystyle\sum_{t_{1}\in B_{b_{1}}}\dots\sum_{t_{l}\in
B_{b_{l}}}v_{b_{1}}^{k_{1}}(t_{1})\dots
v_{b_{l}}^{k_{l}}(t_{l})e^{-\mathrm{i}2\pi\sum_{k=1}^{l}\omega_{k}t_{k}}cum\left(X_{t_{1}}^{(u_{b_{1}})}(\tau_{1}),\dots,X_{t_{l}}^{(u_{b_{l}})}(\tau_{l})\right)+o(1).$

If we have $b_{i}\neq b_{j}$ for any pair $i,j\in\\{1,2,\dots,l\\}$, the last
expression converges to 0 by Assumption 2.1.

For the case $b_{1}=b_{2}=\dots=b_{l}=b$, we write $t_{j}=t_{1}+v_{j}$ for
$j=2,\dots,l$. The expression is then simplified to

$\displaystyle\sum_{v_{2}}\dots\sum_{v_{l}}\exp\left(-\mathrm{i}\sum_{j=2}^{l}v_{j}\omega_{j}\right)cum\left(X_{t_{1}}^{(u_{b})}(\tau_{1}),\dots,X_{t_{1}+v_{l}}^{(u_{b})}(\tau_{l})\right)\sum_{t_{1}}v_{b}(t_{1})v_{b}(t_{1}+v_{2})\times$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\dots\times
v_{b}(t_{1}+v_{l})\exp\left(-\mathrm{i}\sum_{j=1}^{l}\omega_{j}t_{1}\right)$

The result then follows by Lemma P.4.1 an Lemma P.4.2 from [6]. ∎

######  Lemma .2.

Let $\widehat{f}^{(mt)}_{b,\omega}$ be the multitaper periodogram estimator,
as defined in (3.5). Under Assumption 3.1, we have

$\displaystyle\mathbb{E}\hat{f}^{(mt)}_{b,\omega}(\tau,\sigma)=f_{u_{b},\omega}(\tau,\sigma)+O\left(\log(T_{B})/T_{B}\right)+O(1/T).$
(A.1)
$\displaystyle\text{Cov}\left(\hat{f}^{(mt)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(mt)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)=\displaystyle{\left\\{\begin{array}[]{ll}\frac{f_{u_{b},\omega}(\tau_{1},\tau_{2})f_{u_{b},\omega}(\sigma_{1},\sigma_{2})}{K}+O(1/T_{B})+O(1/T),&\text{
if }\omega_{1}=\omega_{2}=\omega\\\ O(1/T_{B})+O(1/T),&\text{ if
}\omega_{1}\neq\omega_{2}.\end{array}\right.}$ (A.3)

###### Proof.

We start by noting that,

$\mathbb{E}\hat{f}^{(k)}_{b,\omega}(\tau,\sigma)=\int_{-1}^{1}|H_{k}(\alpha)|^{2}f_{u_{b},\omega-\alpha}d\alpha+O(1/T),$

where $u_{b}$ is the midpoint of the $b$-th segment. This follows from Theorem
5.2.3 from [6] applied to the approximating series $X^{(u_{b})}_{t}$ within
the $b$-th block. Now write
$f_{u_{b},\omega-\alpha}=f_{u_{b},\omega}+O(|\alpha|)$ and using the form of
$L$, we have

$\mathbb{E}\hat{f}^{(k)}_{b,\omega}(\tau,\sigma)=f_{u_{b},\omega}(\tau,\sigma)+O\left(\log(T_{B})/T_{B}\right)+O(1/T).$

Taking an average of the last expression over different tapers, we get

$\mathbb{E}\hat{f}^{(mt)}_{b,\omega}(\tau,\sigma)=f_{u_{b},\omega}(\tau,\sigma)+O\left(\log(T_{B})/T_{B}\right)+O(1/T).$

To calculate the covariance, first note that

$\text{Cov}\left(\hat{f}^{(k)}_{b_{1},\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(k)}_{b_{2},\omega_{2}}(\tau_{2},\sigma_{2})\right)=0,~{}~{}~{}~{}~{}\text{if
}b_{1}\neq b_{2}.$

Similar calculations as in the proof of Theorem 5.2.8 in [6] yields

$\displaystyle\text{Cov}\left(\hat{f}^{(k)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(k)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)$
$\displaystyle=$
$\displaystyle\left[|H_{k,k}(\omega_{1}+\omega_{2})|^{2}f_{u_{b},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b},\omega_{2}}(\tau_{2},\sigma_{1})+|H_{k,k}(\omega_{1}-\omega_{2})|^{2}f_{u_{b},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{b},\omega_{2}}(\sigma_{1},\sigma_{2})\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}+O(1/T_{B})+O(1/T).$

Therefore if $\omega_{1}\neq\omega_{2}$ we have

$\text{Cov}\left(\hat{f}^{(k)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(k)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)=O(1/T_{B})+O(1/T),$

and

$\text{Cov}\left(\hat{f}^{(k)}_{b,\omega}(\tau_{1},\sigma_{1}),\hat{f}^{(k)}_{b,\omega}(\tau_{2},\sigma_{2})\right)=f_{u_{b},\omega}(\tau_{1},\tau_{2})f_{u_{b},\omega}(\sigma_{1},\sigma_{2})+O(1/T_{B})+O(1/T).$

For local periodograms calculated for different tapers,

$\displaystyle\text{Cov}\left(\hat{f}^{(k)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(l)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)$
$\displaystyle=$
$\displaystyle\left[|H_{k,l}(\omega_{1}+\omega_{2})|^{2}f_{u_{b},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b},\omega_{2}}(\tau_{2},\sigma_{1})+|H_{k,l}(\omega_{1}-\omega_{2})|^{2}f_{u_{b},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{b},\omega_{2}}(\sigma_{1},\sigma_{2})\right]$
$\displaystyle~{}~{}~{}~{}~{}~{}+O(1/T_{B})+O(1/T).$

By the orthogonality of the tapers, $H_{k,l}(0)=0$ for $k\neq l.$ And by
Cauchy-Schwartz, we have
$H_{k,l}(\omega)\leq\sqrt{H_{k,k}(\omega)}\sqrt{H_{l,l}(\omega)}$, and hence
$L$ is indeed an upper bound for $T_{B}H_{k,l}(\omega).$ Therefore, in general
for $k\neq l$,

$\displaystyle\text{Cov}\left(\hat{f}^{(k)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(l)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)=O(1/T_{B})+O(1/T).$

Therefore we have,

$\displaystyle\text{Cov}\left(\hat{f}^{(mt)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(mt)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)=\frac{1}{K^{2}}\sum_{k=1}^{K}\sum_{l=1}^{K}\text{Cov}\left(\hat{f}^{(k)}_{b,\omega_{1}}(\tau_{1},\sigma_{1}),\hat{f}^{(l)}_{b,\omega_{2}}(\tau_{2},\sigma_{2})\right)$
$\displaystyle=$
$\displaystyle\frac{1}{K}\left[f_{u_{b},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{b},\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u_{b},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b},\omega_{2}}(\tau_{2},\sigma_{1})\right]+O(1/T_{B})+O(1/T).$

∎

The next theorem establishes a Central Limit Theorem type result for the
multitaper periodogram estimator and is essential for proving Theorem 3.1 and
Theorem 3.2.

######  Theorem .1.

Consider the processes $E_{b,j}\in L^{2}[0,1]^{2}$ defined as

$E_{b,j}(\tau,\sigma)=\sqrt{K}\left(\hat{f}^{(mt)}_{b,\omega_{j}}(\tau,\sigma)-f_{u_{b},\omega_{j}}(\tau,\sigma)\right)$

for $b=1,2,\dots,B$ and $j=1,2,\dots,J$. For fixed $B$, as $T\to\infty$,
$K\to\infty$ and $T/BK\to\infty$, the finite dimensional distributions of
$\\{E_{b,j}\\}_{j,b}$ converge to a multivariate normal distribution. More
precisely, for all
$(\tau_{1},\sigma_{1}),\dots,(\tau_{d},\sigma_{d})\in[0,1]^{2}$ and for all
$d\in\mathbb{N}$,

$\\{E_{b,j}(\tau_{1},\sigma_{1}),\dots,E_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}\stackrel{{\scriptstyle
d}}{{\to}}\\{Z_{b,j}(\tau_{1},\sigma_{1}),\dots,Z_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$

where
$\\{Z_{b,j}(\tau_{1},\sigma_{1}),\dots,Z_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$
is a multivariate normal random vector with zero mean and covariance structure

$\displaystyle\text{Cov}\left(Z_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),Z_{b_{2},j_{2}}(\tau_{2},\sigma_{2})\right)$
$\displaystyle=\left\\{\begin{array}[]{ll}f_{u_{b},\omega_{j_{1}}}(\tau_{1},\tau_{2})f_{u_{b},\omega_{j_{2}}}(\sigma_{1},\sigma_{2})+f_{u_{b},\omega_{j_{1}}}(\tau_{1},\sigma_{2})f_{u_{b},\omega_{j_{2}}}(\tau_{2},\sigma_{1}),&\text{if
}~{}b_{1}=b_{2}=b,\\\ 0,&\text{otherwise.}\end{array}\right.$

###### Proof.

We will show the cumulants of the vector
$\\{E_{b,j}(\tau_{1},\sigma_{1}),\dots,E_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$
converge to the cumulants of the vector
$\\{Z_{b,j}(\tau_{1},\sigma_{1}),\dots,Z_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$.
As the cumulants of order $l$ of the Gaussian distribution are zero for $l>2$,
we will show that as $T\to\infty$ and $K\to\infty$,

$\displaystyle
cum\left(E_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),\dots,\dots,E_{b_{l},j_{l}}(\tau_{l},\sigma_{l})\right)$
$\displaystyle=\left\\{\begin{array}[]{ll}o(1)&\text{if }l\neq 2\\\
\text{Cov}\left(Z_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),Z_{b_{2},j_{2}}(\tau_{2},\sigma_{2})\right)+o(1)&\text{if
}l=2.\end{array}\right.$

Note that the equality for $l=1$ and $2$ follow from the earlier expectation
and variance calculation. Hence we have to show the result for $l\geq 3.$

Note that

$\displaystyle
cum\left(E_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),\dots,\dots,E_{b_{l},j_{l}}(\tau_{l},\sigma_{l})\right)$
$\displaystyle=K^{l/2}cum\left(\hat{f}^{(mt)}_{b_{1},\omega_{j_{1}}}(\tau_{1},\sigma_{1}),\dots,\dots,\hat{f}^{(mt)}_{b_{l},\omega_{j_{l}}}(\tau_{l},\sigma_{l})\right)$
$\displaystyle=\frac{1}{K^{l/2}}\sum_{k_{1}}\dots\sum_{k_{l}}cum\left(\hat{f}^{(k_{1})}_{b_{1},\omega_{j_{1}}}(\tau_{1},\sigma_{1}),\dots,\dots,\hat{f}^{(k_{l})}_{b_{l},\omega_{j_{l}}}(\tau_{l},\sigma_{l})\right)$
$\displaystyle=\frac{1}{K^{l/2}}\sum_{k_{1}}\dots\sum_{k_{l}}cum\left(Y_{11}Y_{12},\dots,Y_{l1}Y_{l2}\right)$

where $Y_{i1}=\widetilde{X}_{T}^{(k_{i}),b_{i},\omega_{i}}(\tau_{i})$ and
$Y_{i2}=\widetilde{X}_{T}^{(k_{i}),b_{i},\omega_{i}}(\sigma_{i})$.
Using Theorem 2.3.2 from [6], the last quantity is equal to

$\frac{1}{K^{l/2}}\sum_{k_{1}}\dots\sum_{k_{l}}\sum_{\nu}cum(Y_{ij}:ij\in\nu_{1})\dots
cum(Y_{ij}:ij\in\nu_{p})=:\sum_{\nu}C(\nu)$

where the sum is over all indecomposable partitions
$\nu=\nu_{1}\cup\nu_{2}\cup\dots\cup\nu_{p}$ of

$\begin{array}[]{cc}(1,1)&(1,2)\\\ (2,1)&(2,2)\\\ \vdots&\vdots\\\
(l,1)&(l,2).\end{array}$

As there are a finite number of partitions, it is enough to show $C(\nu)=o(1)$
for all indecomposible partitions $\nu$ for $l>2$.

To this end, we note that the function
$H_{k_{1},k_{2},\dots,k_{m}}(\omega)=O(T_{B}^{-m/2})$ if $\omega\neq 0.$ By
the orthonormality and symmetry of the wave function,
$H_{k_{1},k_{2},\dots,k_{m}}(\omega)=O(T_{B}^{1-m/2})$ if $\omega=0$ and all
distinct $k_{i}$’s appear an even number of times in the index
$k_{1},\dots,k_{m}$, and $H_{k_{1},k_{2},\dots,k_{m}}(\omega)=0$, if
$\omega=0$ and any $k_{i}$ appear an odd number of times.

Let $|\nu_{i}|$ denote the number of elements of the set $\nu_{i}$. Note that
$\nu$ is a partition of a set of $2l$ elements, and therefore
$\sum_{i=1}^{p}|\nu_{i}|=2l.$ By Lemma .1 and the property of the
$H_{k_{1},\dots,k_{l}}$ functions, we note that $C(\nu)=0$ if any $\nu_{i}$ in
$\nu$ has at least one $k_{i}$ an odd number of times and otherwise

$C(\nu)=\frac{1}{K^{l/2}}\sum_{k_{1}}\dots\sum_{k_{l}}\prod_{i=1}^{p}O\left(T_{B}^{-|\nu_{i}|/2+1\left(\sum_{j\in\nu_{i}}\omega_{j}=0\mod
2\pi\right)}\right)=\sum_{r=1}^{l}O(K^{r-l/2})O\left(T_{B}^{-l+s(\nu)}\right),$

where $r$ is the distinct number of $k_{i}$’s in a collection
$k_{1},k_{2},\dots,k_{l}$ and $s(\nu)$ = the number of $\nu_{i}$ in $\nu$,
such that $\sum_{j\in\nu_{i}}\omega_{j}=0.$

If $r>l/2$, then at least one of the $k_{i}$’s appear just once, and therefore
one of the sets in the partition must have a single occurrence of that index.
So it is enough to consider the case where $r\leq l/2$. Now consider the
possibilities for the $O\left(T_{B}^{-l+s(\nu)}\right)$ term.

Case 1: If $p<l$ then $s(\nu)\leq p<l$, and therefore
$O\left(T_{B}^{-l+s(\nu)}\right)=o(1).$

Case 2: If $p>l$, at least $2(p-l)$ sets of the partitions have just one
element. (To see this, suppose $l_{1}$ is the number of sets with one element.
Then we have $2l\geq l_{1}+2(p-l_{1})$.) For all those one element sets
$\sum_{j}\omega_{j}\neq 0\mod 2\pi.$ Therefore $s(\nu)\leq p-2(p-l)=2l-p<l$
and consequently $O\left(T_{B}^{-l+s(\nu)}\right)=o(1).$

Case 3: If $p=l$ and at least one set in the partition has a single element,
for that set $\sum_{j}\omega_{j}\neq 0\mod 2\pi$, therefore $s(\nu)\leq l-1$
and $O\left(T_{B}^{-l+s(\nu)}\right)=o(1).$

Case 4: Consider the case where $p=l$ and all the partitions have 2 elements.
Note that as $s(\nu)\leq p$, if $r<l/2$ the product
$O(K^{r-l/2})O\left(T_{B}^{-l+s(\nu)}\right)=o(1).$ Therefore it is enough to
consider the case where $l$ is even and $r=l/2.$ If $l>2$, this means there
are at least two distinct $k_{i}$ in the collection and by indecomposibility,
one of the sets $\nu_{i}$ in the partition must have two distinct $k_{i}$,
making $C(\nu)=0$. ∎

###  .2 Proof of Results in Section 3

#### Proof of Lemma 3.1

Noting the definition of $\widehat{g}_{u,\omega}$ in (3.6), the result follows
from Lemma .2 and some simple algebra.

######  Lemma .3.

Consider the processes $H_{b,j}\in L^{2}[0,1]^{2}$ defined as

$H_{b,j}(\tau,\sigma)=\sqrt{K}\left(\widehat{g}_{b/B,\omega_{j}}(\tau,\sigma)-g_{u_{b},\omega_{j}}(\tau,\sigma)\right)$

for $b=1,2,\dots,B$ and $j=1,2,\dots,J$. For fixed $B$, as $T\to\infty$,
$K\to\infty$ and $T/BK\to\infty$, the finite dimensional distributions of
$\\{H_{b,j}\\}_{j,b}$ converges to a multivariate normal distribution. More
precisely, for all
$(\tau_{1},\sigma_{1}),\dots,(\tau_{d},\sigma_{d})\in[0,1]^{2}$ and for all
$d\in\mathbb{N}$

$\\{H_{b,j}(\tau_{1},\sigma_{1}),\dots,H_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}\stackrel{{\scriptstyle
d}}{{\to}}\\{Z^{\prime}_{b,j}(\tau_{1},\sigma_{1}),\dots,Z^{\prime}_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$

where
$\\{Z^{\prime}_{b,j}(\tau_{1},\sigma_{1}),\dots,Z^{\prime}_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$
is a multivariate normal random vector with zero mean and covariance structure

$\displaystyle\text{Cov}\left(Z^{\prime}_{b,j_{1}}(\tau_{1},\sigma_{1}),Z^{\prime}_{b,j_{2}}(\tau_{2},\sigma_{2})\right)=\left(1-\frac{2}{B}\right)F(u_{b},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})$
$\displaystyle~{}~{}\hskip
162.6075pt+\frac{1}{B^{2}}\sum_{l=1}^{B}F(u_{l},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}).$
$\displaystyle\text{Cov}\left(Z^{\prime}_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),Z^{\prime}_{b_{2},j_{2}}(\tau_{2},\sigma_{2})\right)=-\frac{1}{B}F(u_{b_{1}},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{B}F(u_{b_{2}},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})+\frac{1}{B^{2}}\sum_{l=1}^{B}F(u_{b_{l}},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}),$
(A.4)

where $b_{1}\neq b_{2}$ and

$F(u,\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}):=f_{u,\omega_{1}}(\tau_{1},\tau_{2})f_{u,\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u,\omega_{2}}(\tau_{1},\sigma_{2})f_{u,\omega_{2}}(\tau_{2},\sigma_{1}).$
(A.5)

###### Proof.

In view of Lemma 3.1 it is enough to show that the joint cumulants of
$\\{H_{b,j}(\tau_{1},\sigma_{1}),\dots,H_{b,j}(\tau_{d},\sigma_{d})\\}_{b,j}$
for $b=1,2,\dots,B$ and $j=1,2,\dots,J$ of order $>2$ converges to $0$, as
$T\to\infty,$ $K\to\infty$ and $T/K\to\infty$. Note that by definition of
$\widehat{g}$, we have

$H_{b,j}(\tau,\sigma)=E_{b,j}(\tau,\sigma)-\frac{1}{B}\sum_{b=1}^{B}E_{b,j}(\tau,\sigma).$

Therefore using the linearity of cumulants (Theorem 2.3.1 (i) & (iii) from
[6]) we can write
$cum\left(H_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),\dots,\dots,H_{b_{l},j_{l}}(\tau_{l},\sigma_{l})\right)$
as sum of $2^{l}$ terms, where each term is of the form

$cum\left(\widetilde{E}_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),\dots,\dots,\widetilde{E}_{b_{l},j_{l}}(\tau_{l},\sigma_{l})\right)$

where $\widetilde{E}_{b,j}(\tau,\sigma)$ is either ${E}_{b,j}(\tau,\sigma)$ or
the average $\frac{1}{B}\sum_{b=1}^{B}{E}_{b,j}(\tau,\sigma)$. Without loss of
generality, consider a term where the first $k$
$\widetilde{E}_{b,j}(\tau,\sigma)$’s are the averages and the last $(l-k)$ are
${E}_{b,j}$. That particular term can then be simplified to

$\frac{1}{B^{k}}\sum_{b_{1}=1}^{B}\dots\sum_{b_{k}=1}^{B}cum\left(E_{b_{1},j_{1}}(\tau_{1},\sigma_{1}),\dots,\dots,E_{b_{l},j_{l}}(\tau_{l},\sigma_{l})\right).$

As $B$ is finite, the last sum converges to 0 by Theorem .1. As all the
$2^{l}$ terms converge to 0 and $l$ is finite, this imples the convergence of
joint cumulants of $H_{b,j}(\tau,\sigma)$ for order $l>2$ to zero.

∎

#### Proof of Theorem 3.1

We write
$Q_{\omega_{0},\delta}=\int_{0}^{1}\int_{0}^{1}Q_{\omega_{0},\delta}(\tau,\sigma)d\tau
d\sigma,$ where

$Q_{\omega_{0},\delta}(\tau,\sigma)=\sum_{b=1}^{B}\left(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma)-\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma)\right)^{2}.$

Noting than the functional $:L^{2}\left([0,1]^{2}\right)\mapsto\mathbb{R}$ is
continuous it is enough to show that

$Q_{\omega_{0},\delta}(\tau,\sigma)\stackrel{{\scriptstyle
d}}{{=}}\frac{1}{K}\sum_{b=1}^{B}\left(\mathcal{G}_{b}^{2}(\tau,\sigma)+o_{p}(1)\right),$

This will be proved in two steps. Specifically, we will show that:

  * (i)

The result holds in finite dimensional distribution, i.e., for any $k$ and any
fixed $(\tau_{1},\sigma_{1}),\dots,(\tau_{k},\sigma_{k})\in[0,1]^{2}$ the
distribution of the random vector

$(Q_{\omega_{0},\delta}(\tau_{1},\sigma_{1}),\dots,Q_{\omega_{0},\delta}(\tau_{k},\sigma_{k}))$

is asymptotically equal to the distribution of

$\frac{1}{K}\left(\sum_{b=1}^{B}\left(\mathcal{G}_{b}^{2}(\tau_{1},\sigma_{1})+o_{p}(1)\right),\dots,\sum_{b=1}^{B}\left(\mathcal{G}_{b}^{2}(\tau_{1},\sigma_{1})+o_{p}(1)\right)\right).$

  * (ii)

The process
$\\{Q_{\omega_{0},\delta}(\tau,\sigma)\\}_{(\tau,\sigma)\in[0,1]^{2}}$ is
asymptotically tight as a process in $L^{2}([0,1]^{2}).$

Without the loss of generality, we will prove (i) for $k=1$, the result for
general $k$ can be proved similarly with some more notations. Note that the
scan statistics can be written as

$Q_{\omega_{0},\delta}(\tau,\sigma)=\frac{1}{K}\sum_{b=1}^{B}A_{b,K,T}$

where

$A_{b,K,T}(\tau,\sigma)=K\left(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma)-\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma)\right)^{2}.$

Therefore it is enough to show the process $A_{b,k,T}\in L^{2}([0,1]^{2})$
converges in distribution to $\mathcal{G}_{b}^{2}$, where $\mathcal{G}_{b}$ is
the Gaussian process defined in the statement of the Theorem, uniformly over
$b\in\\{1,2,\dots,B\\}.$ In order to establish that write

$\displaystyle A_{b,K,T}(\tau,\sigma)=$ $\displaystyle
K\left(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma)-\mathbb{E}(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma))-\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma)+\mathbb{E}(\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma))\right)^{2}$
$\displaystyle+2K\left(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma)-\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma)\right)\left(\mathbb{E}(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma))-\mathbb{E}(\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma))\right)$
$\displaystyle-K\left(\mathbb{E}(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma))-\mathbb{E}(\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma))\right)^{2}.$
(A.6)

Suppose $L_{\delta}$ is the number of frequencies $\omega_{j}$ in the interval
$[\omega_{0},\omega_{0}+\delta)$. Using Lemma 3.1, under $H_{0}$ we have,

$\mathbb{E}(\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma))=\frac{1}{L_{\delta}}\sum_{\omega_{j}\in[\omega_{0},\omega_{0}+\delta)}\mathbb{E}(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma))=g_{u_{b},\omega_{0}+\delta}(\tau,\sigma)+O(\log(T_{B})/T_{B})+O(1/T).$

Therefore by Lemma 3.1, the last two terms of (.2) is of the order
$O\left(K\times\frac{\log^{2}{T_{B}}}{T_{B}^{2}}\right)+O\left(\frac{K}{T^{2}}\right)$,
which converges to zero under Assumption 3.1. Note that the order of these
residuals are independent of the choice of block $b$.

The first term of (.2) can be written as
$T^{2}(E_{b,j_{1}}(\tau,\sigma),\dots,E_{b,j_{L_{\delta}+1}}(\tau,\sigma)),$
where $E_{b,j}$ is the process defined in Theorem .1, the set of frequencies
$\\{\omega_{j_{1}},\dots,\omega_{j_{L_{\delta}}}\\}=\\{\omega:\omega\in[\omega_{0},\omega_{0}+\delta)\\}$,
$\omega_{j_{L_{\delta}+1}}=\omega_{0}+\delta$ and the function
$T:\mathbb{R}^{L_{\delta}+1}\mapsto\mathbb{R}$ is defined as

$T(x_{1},\dots,x_{L_{\delta}+1})=x_{L_{\delta}+1}-\frac{1}{L_{\delta}}\sum_{i=1}^{L_{\delta}}x_{i}.$

Therefore an application of the Delta method along with Theorem .1 guarantee

$T(E_{b,j_{1}}(\tau,\sigma),\dots,E_{b,j_{L_{\delta}+1}}(\tau,\sigma))\stackrel{{\scriptstyle
d}}{{\to}}\mathcal{G}_{b}(\tau,\sigma)$

where $\mathcal{G}_{b}$ is a zero mean Gaussian process with covariance kernel
given in Theorem 3.1. The weak convergence of $A_{b,K,T}(\tau,\sigma)$ to
$\mathcal{G}_{b}^{2}(\tau,\sigma)$ then follows by the Continuous Mapping
Theorem and Slutzky’s Theorem. This in turn proves the asymptotic finite
dimensional distributional equivalence in (i).

Part (ii) follows from Lemma .1 by (4.3) of Theorem 2 from [9].

#### Proof of Theorem 3.2

The proof is similar to the proof of Theorem 3.1. The only difference is in
the treatment of the residual (second and third) terms of
$A_{b,K,T}(\tau,\sigma)$ defined in (.2). Note that under the alternative
specified in the statement of the Theorem,

$\displaystyle\mathbb{E}(\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma))=$
$\displaystyle\frac{1}{L_{\delta}}\sum_{\omega_{j}\in[\omega_{0},\omega_{0}+\delta)}\mathbb{E}(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma))$
$\displaystyle=$
$\displaystyle\frac{\\#\\{\omega:\omega\in(\omega_{0},\omega^{*})\\}}{L_{\delta}}g^{(1)}_{u_{b}}(\tau,\sigma)+\frac{\\#\\{\omega:\omega\in[\omega^{*},\omega_{0}+\delta)\\}}{L_{\delta}}g^{(2)}_{u_{b}}(\tau,\sigma)$
$\displaystyle~{}~{}+O(\log(T_{B})/T_{B})+O(1/T).$ $\displaystyle=$
$\displaystyle\frac{\omega^{*}-\omega_{0}}{\delta}g^{(1)}_{u_{b}}(\tau,\sigma)+\frac{\omega_{0}+\delta-\omega^{*}}{\delta}g^{(2)}_{u_{b}}(\tau,\sigma)+O(\log(T_{B})/T_{B})+O(1/T).$
$\displaystyle=$ $\displaystyle
g^{(2)}_{u_{b}}(\tau,\sigma)+\frac{\omega^{*}-\omega_{0}}{\delta}\left(g^{(1)}_{u_{b}}(\tau,\sigma)-g^{(2)}_{u_{b}}(\tau,\sigma)\right)+O(\log(T_{B})/T_{B})+O(1/T).$

Note that the second equality follows from the fact that $\omega_{j}$’s are
chosen equally spaced. Therefore the residual is

$K\left(\mathbb{E}(\widehat{g}_{b/B,\omega_{0}+\delta}(\tau,\sigma))-\mathbb{E}(\widetilde{g}_{b/B,\omega_{0},\delta}(\tau,\sigma))\right)^{2}+o_{p}(1)=\frac{\omega^{*}-\omega_{0}}{\delta}\left(g^{(1)}_{u_{b}}(\tau,\sigma)-g^{(2)}_{u_{b}}(\tau,\sigma)\right)+o_{p}(1).$

The the rest of the proof is similar the proof of Theorem 3.1.

#### Proof of Theorem 3.3

We write

$Q^{0}(\omega_{1},\omega_{2})=\frac{1}{K}\sum_{b=1}^{B}\int_{\omega_{1}}^{\omega_{2}}\int_{0}^{1}\int_{0}^{1}\left(\sqrt{K}\widehat{g}_{b/B,\omega}(\tau,\sigma)\right)^{2}d\tau
d\sigma
d\omega=\int_{\omega_{1}}^{\omega_{2}}\int_{0}^{1}\int_{0}^{1}Q_{0}(\omega,\tau,\sigma),$

where

$Q_{0}(\omega,\tau,\sigma)=\frac{1}{K}\sum_{b=1}^{B}\left(\sqrt{K}\widehat{g}_{b/B,\omega}(\tau,\sigma)\right)^{2}.$

An application of Lemma 3.1 along with continuous mapping theorem guarantees
under $H_{0}$ that the finite dimensional distributions of
$Q_{0}(\omega,\tau,\sigma)$ are asymptotically equivalent to
$\frac{1}{K}\sum_{b=1}^{B}\mathcal{H}_{b}^{2}(\tau,\sigma).$ The rest of the
proof follows is similar as to the proof of Theorem 3.1.

#### Proof of Lemma 3.2

Using similar expansion as in the proof of Theorem 3.3, under the alternative
we write

$\displaystyle Q_{0}(\omega,\tau,\sigma)=$
$\displaystyle\frac{1}{K}\sum_{b=1}^{B}\left(\sqrt{K}\widehat{g}_{b/B,\omega}(\tau,\sigma)\right)^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{b=1}^{B}\left(\sqrt{K}\left[\widehat{g}_{b/B,\omega}(\tau,\sigma)-g_{u}(\tau,\sigma)+g_{u}(\tau,\sigma)\right]\right)^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{b=1}^{B}\left(\sqrt{K}\left[\widehat{g}_{b/B,\omega}(\tau,\sigma)-g_{u}(\tau,\sigma)\right]\right)^{2}+Bg_{u}^{2}(\tau,\sigma)$
$\displaystyle+\frac{2g_{u}(\tau,\sigma)}{\sqrt{K}}\sum_{b=1}^{B}\sqrt{K}\left[\widehat{g}_{b/B,\omega}(\tau,\sigma)-g_{u}(\tau,\sigma)\right].$

As the second term dominates under the asymptotic scheme in Assumption 3.1,
the quantity $Q_{0}(\omega,\tau,\sigma)=O_{p}(B)$. Taking integral over
$\omega,\tau,\sigma$ we have $Q^{0}(\omega_{1},\omega_{2})=O_{p}(B)$ and the
result follows.

###  .3 Asymptotic Distribution Properties

By Lemma 3.1 the covariance of the process $\mathcal{H}_{b}$ defined in
Theorem 3.3 is given by

$\displaystyle\text{Cov}(\mathcal{H}_{b}(\tau_{1},\sigma_{1}),\mathcal{H}_{b}(\tau_{2},\sigma_{2}))=C(b_{1},b_{2},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}),$
(A.7)

where

$\displaystyle
C(b,b,\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})=$
$\displaystyle(1-2/B)[f_{u_{b},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{b},\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u_{b},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b},\omega_{2}}(\tau_{2},\sigma_{1})]$
(A.8)
$\displaystyle+\frac{1}{B^{2}}\sum_{l=1}^{B}[f_{u_{l},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{l},\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u_{l},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{l},\omega_{2}}(\tau_{2},\sigma_{1})]$
$\displaystyle{C}\left(b_{1},b_{2},\omega_{1},\omega_{2},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}\right)=$
$\displaystyle-\frac{1}{B}\left[f_{u_{b_{1}},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{b_{1}},\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u_{b_{1}},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b_{1}},\omega_{2}}(\tau_{2},\sigma_{1})\right]$
$\displaystyle-\frac{1}{B}\left[f_{u_{b_{2}},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{b_{2}},\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u_{b_{2}},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b_{2}},\omega_{2}}(\tau_{2},\sigma_{1})\right]$
$\displaystyle~{}+\frac{1}{B^{2}}\sum_{l=1}^{B}\left[f_{u_{b_{l}},\omega_{1}}(\tau_{1},\tau_{2})f_{u_{l},\omega_{2}}(\sigma_{1},\sigma_{2})+f_{u_{b_{l}},\omega_{1}}(\tau_{1},\sigma_{2})f_{u_{b_{l}},\omega_{2}}(\tau_{2},\sigma_{1})\right],$

where $b_{1}\neq b_{2}$.
The covariance structure of the process $\mathcal{G}_{b}$ is given by

$\displaystyle
Cov(\mathcal{G}_{b}(\tau_{1},\sigma_{1}),\mathcal{G}_{b}(\tau_{2},\sigma_{2}))=$
$\displaystyle
C(b,b,\omega_{0}+\delta,\omega_{0}+\delta,\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}))$
$\displaystyle+\frac{1}{L_{\delta}^{2}}\sum_{j=1}^{L_{\delta}}\sum_{k=1}^{L_{\delta}}C(b,b,\omega_{j},\omega_{k},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}))$
$\displaystyle-\frac{1}{L_{\delta}}\sum_{j=1}^{L}C(b,b,\omega_{0}+\delta,\omega_{J},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}))$
$\displaystyle-\frac{1}{L_{\delta}}\sum_{j=1}^{L}C(b,b,\omega_{j},\omega_{0}+\delta,\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})),$
(A.9)

and for $b_{1}\neq b_{2}$,

$\displaystyle
Cov(\mathcal{G}_{b_{1}}(\tau_{1},\sigma_{1}),\mathcal{G}_{b_{2}}(\tau_{2},\sigma_{2}))=$
$\displaystyle
C(b_{1},b_{2},\omega_{0}+\delta,\omega_{0}+\delta,\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}))$
$\displaystyle+\frac{1}{L_{\delta}^{2}}\sum_{j=1}^{L_{\delta}}\sum_{k=1}^{L_{\delta}}C(b_{1},b_{2},\omega_{j},\omega_{k},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}))$
$\displaystyle-\frac{1}{L_{\delta}}\sum_{j=1}^{L}C(b_{1},b_{2},\omega_{0}+\delta,\omega_{j},\tau_{1},\sigma_{1},\tau_{2},\sigma_{2})))$
$\displaystyle-\frac{1}{L_{\delta}}\sum_{j=1}^{L}C(b_{1},b_{2},\omega_{j},\omega_{0}+\delta,\tau_{1},\sigma_{1},\tau_{2},\sigma_{2}))$
(A.10)

where, $C$ is as defined in (A.8).

######  Lemma .4.

As $B\to\infty$, the quantities
$\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|^{2}$ and
$\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{H}_{b}\|^{2}$ are $O_{p}(1).$

###### Proof.

We will show this for the case of $\mathcal{G}_{b}$. The proof for
$\mathcal{H}_{b}$ is similar.

Note that

$\displaystyle\mathbb{E}\left(\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|^{2}\right)=\frac{1}{B}\sum_{b=1}^{B}\mathbb{E}\left(\|\mathcal{G}_{b}\|^{2}\right)$
$\displaystyle=\frac{1}{B}\sum_{b=1}^{B}\mathbb{E}\left(\int_{0}^{1}\int_{0}^{1}\mathcal{G}_{b}^{2}(\tau,\sigma)d\tau
d\sigma\right)$
$\displaystyle=\frac{1}{B}\sum_{b=1}^{B}\int_{0}^{1}\int_{0}^{1}\mathbb{E}\mathcal{G}_{b}^{2}(\tau,\sigma)d\tau
d\sigma$
$\displaystyle\to\int_{0}^{1}\int_{0}^{1}\int_{0}^{1}\left[f_{u,\omega_{0}+\delta}(\tau,\tau)f_{u,\omega_{0}+\delta}(\sigma,\sigma)+f_{u,\omega_{0}+\delta}^{2}(\tau,\sigma)\right]d\tau
d\sigma du$
$\displaystyle~{}~{}+\frac{1}{L_{\delta}^{2}}\sum_{j=1}^{L_{\delta}}\sum_{k=1}^{L_{\delta}}\iiint\left[f_{u,\omega_{j}}(\tau,\tau)f_{u,\omega_{k}}(\sigma,\sigma)+f_{u,\omega_{j}}(\tau,\sigma)f_{u,\omega_{k}}(\tau,\sigma)\right]d\tau
d\sigma du$
$\displaystyle~{}~{}-\frac{2}{L_{\delta}}\sum_{j=1}^{L_{\delta}}\iiint\left[f_{u,\omega_{j}}(\tau,\tau)f_{u,\omega_{0}+\delta}(\sigma,\sigma)+f_{u,\omega_{j}}(\tau,\sigma)f_{u,\omega_{0}+\delta}(\tau,\sigma)\right]d\tau
d\sigma du,$

as $B\to\infty$. As $f_{u,\omega}(.)$ is continuous in $u$ and $f_{u,\omega}$
is square integrable for all $u\in[0,1]$ and $\omega\in(0,0.5]$, the integrals
in the limit are finite. Note that we can exchange the integral and
expectation in the second line by Fubini’s Theorem as the double integral is
finite.

$\displaystyle
Var\left(\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|^{2}\right)=\frac{1}{B^{2}}\sum_{b=1}^{B}Var\left(\|\mathcal{G}_{b}\|^{2}\right)+\frac{1}{B^{2}}\sum_{b_{1}=1}^{B}\sum_{b_{2}=1}^{B}Cov\left(\|\mathcal{G}_{b_{1}}\|^{2},\|\mathcal{G}_{b_{2}}\|^{2}\right)$

The first term can be simplified as

$\displaystyle\frac{1}{B^{2}}\sum_{b=1}^{B}Var\left(\|\mathcal{G}_{b}\|^{2}\right)$
$\displaystyle\leq\frac{1}{B^{2}}\sum_{b=1}^{B}\mathbb{E}\left(\|\mathcal{G}_{b}\|^{4}\right)=\frac{1}{B^{2}}\sum_{b=1}^{B}\mathbb{E}\left(\int_{0}^{1}\int_{0}^{1}\mathcal{G}_{b}^{2}(\tau,\sigma)d\tau
d\sigma\right)^{2}$
$\displaystyle\leq\frac{1}{B^{2}}\sum_{b=1}^{B}\left(\int_{0}^{1}\int_{0}^{1}\mathbb{E}\mathcal{G}_{b}^{2}(\tau,\sigma)d\tau
d\sigma\right)^{2}=O(1/B).$

Similarly with some standard algebra we can show that

$\displaystyle\frac{1}{B^{2}}\sum_{b_{1}=1}^{B}\sum_{b_{2}=1}^{B}Cov\left(\|\mathcal{G}_{b_{1}}\|^{2},\|\mathcal{G}_{b_{2}}\|^{2}\right)=\frac{1}{B^{2}}\sum_{b_{1}=1}^{B}\sum_{b_{2}=1}^{B}\mathbb{E}\left(\|\mathcal{G}_{b_{1}}\|^{2}\|\mathcal{G}_{b_{2}}\|^{2}\right)=O(1/B).$

Therefore as $B\to\infty$,
$\mathbb{E}\left(\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|^{2}\right)=O(1)$
and $Var\left(\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|^{2}\right)\to 0$
and hence $\frac{1}{B}\sum_{b=1}^{B}\|\mathcal{G}_{b}\|^{2}$ is $O_{p}(1)$. ∎

## Acknowledgements

Research reported in this publication was supported by the National Institute
Of General Medical Sciences of the National Institutes of Health under Award
Number R01GM140476. The content is solely the responsibility of the authors
and does not necessarily represent the official views of the National
Institutes of Health.

[id=suppA] R code for "Adaptive Frequency Band Analysis for Functional Time
Series" [url]https://github.com/sbruce23/fEBA R code, a quick start demo, and
descriptions of all functions and parameters needed to generate simulated data
introduced in Section 4.1 and to implement the proposed method on data for use
in practice can be downloaded from GitHub at this link.

## References

  * Aue and van Delft, [2020] Aue, A. and van Delft, A. (2020).  Testing for stationarity of functional time series in the frequency domain.  The Annals of Statistics, 48(5):2505 – 2547.
  * Barry et al., [2007] Barry, R. J., Clarke, A. R., Johnstone, S. J., Magee, C. A., and Rushby, J. A. (2007).  EEG differences between eyes-closed and eyes-open resting conditions.  Clinical Neurophysiology, 118(12):2765 – 2773.
  * Barry and De Blasio, [2017] Barry, R. J. and De Blasio, F. M. (2017).  EEG differences between eyes-closed and eyes-open resting remain in healthy ageing.  Biological Psychology, 129:293 – 304.
  * Benoit et al., [2000] Benoit, O., Daurat, A., and Prado, J. (2000).  Slow (0.7–2 hz) and fast (2–4 hz) delta components are differently correlated to theta, alpha and beta frequency bands during NREM sleep.  Clinical Neurophysiology, 111(12):2103 – 2106.
  * Billman, [2011] Billman, G. (2011).  Heart rate variability - A historical perspective.  Frontiers in Physiology, 2:86.
  * Brillinger, [2002] Brillinger, D. R. (2002).  Time Series: Data Analysis and Theory.  Philadelphia: SIAM.
  * Bruce et al., [2020] Bruce, S. A., Tang, C. Y., Hall, M. H., and Krafty, R. T. (2020).  Empirical frequency band analysis of nonstationary time series.  Journal of the American Statistical Association, 115:1933–1945.
  * Causeur et al., [2009] Causeur, D., Kloareg, M., and Friguet, C. (2009).  Control of the FWER in multiple testing under dependence.  Communications in Statistics - Theory and Methods, 38(16-17):2733–2747.
  * Cremers and Kadelka, [1986] Cremers, H. and Kadelka, D. (1986).  On weak convergence of integral functionals of stochastic processes with applications to processes taking paths in $L_{p}^{E}$.  Stochastic processes and their applications, 21(2):305–317.
  * Dahlhaus, [1985] Dahlhaus, R. (1985).  Asymptotic normality of spectral estimates.  Journal of Multivariate Analysis, 16(3):412–431.
  * Doppelmayr et al., [1998] Doppelmayr, M., Klimesch, W., Pachinger, T., and Ripper, B. (1998).  Individual differences in brain dynamics: Important implications for the calculation of event-related band power.  Biological Cybernetics, 79(1):49-57.
  * Efron, [2007] Efron, B. (2007).  Correlation and large-scale simultaneous significance testing.  Journal of the American Statistical Association, 102(477):93–103.
  * Glendinning and Fleet, [2007] Glendinning, R. and Fleet, S. (2007).  Classifying functional time series.  Signal Processing, 87(1):79 – 100.
  * Grenander, [1950] Grenander, U. (1950).  Stochastic processes and statistical inference.  Arkiv för matematik, 1(3):195–277.
  * Hochberg, [1988] Hochberg, Y. (1988).  A sharper Bonferroni procedure for multiple tests of significance.  Biometrika, 75(4):800–802.
  * Horváth et al., [2020] Horváth, L., Liu, Z., Rice, G., and Wang, S. (2020).  A functional time series analysis of forward curves derived from commodity futures.  International Journal of Forecasting, 36(2):646 – 665.
  * Klimesch, [1999] Klimesch, W. (1999).  EEG alpha and theta oscillations reflect cognitive and memory performance: A review and analysis.  Brain Research Reviews, 29(2):169-195.
  * Klimesch et al., [1998] Klimesch, W., Doppelmayr, M., Russegger, H., Pachinger, T., and Schwaiger, J. (1998).  Induced alpha band power changes in the human EEG and attention.  Neuroscience Letters, 244(2):73-76.
  * Lenssen et al., [2019] Lenssen, N. J. L., Schmidt, G. A., Hansen, J. E., Menne, M. J., Persin, A., Ruedy, R., and Zyss, D. (2019).  Improvements in the gistemp uncertainty model.  Journal of Geophysical Research: Atmospheres, 124(12):6307–6326.
  * Malik et al., [1996] Malik, M., Bigger, J. T., Camm, A. J., Kleiger, R. E., Malliani, A., Moss, A. J., and Schwartz, P. J. (1996).  Heart rate variability. Standards of measurement, physiological interpretation, and clinical use.  European Heart Journal, 17(3):354-381.
  * Nacy et al., [2016] Nacy, S. M., Kbah, S. N., Jafer, H. A., and Al-Shaalan, I. (2016).  Controlling a servo motor using eeg signals from the primary motor cortex.  American Journal of Biomedical Engineering, 6(5):139–146.
  * Newson and Thiagarajan, [2019] Newson, J. J. and Thiagarajan, T. C. (2019).  Eeg frequency bands in psychiatric disorders: A review of resting state studies.  Frontiers in Human Neuroscience, 12:521.
  * Panaretos and Tavakoli, [2013] Panaretos, V. M. and Tavakoli, S. (2013).  Cramér–karhunen–loève representation and harmonic principal component analysis of functional time series.  Stochastic Processes and their Applications, 123(7):2779–2807.
  * Panaretos et al., [2013] Panaretos, V. M., Tavakoli, S., et al. (2013).  Fourier analysis of stationary time series in function space.  The Annals of Statistics, 41(2):568–603.
  * R Core Team, [2020] R Core Team (2020).  R: A Language and Environment for Statistical Computing.  R Foundation for Statistical Computing, Vienna, Austria.
  * Ramsay, [1982] Ramsay, J. (1982).  When the data are functions.  Psychometrika, 47(4):379–396.
  * Ramsay and Danzell, [1991] Ramsay, J. O. and Danzell, C. J. (1991).  Some tools for functional data analysis.  Journal of the Royal Statistical Society, B, 53:539–572.
  * Rand, [1971] Rand, W. M. (1971).  Objective criteria for the evaluation of clustering methods.  Journal of the American Statistical Association, 66(336):846–850.
  * Rao, [1958] Rao, C. R. (1958).  Some statistical methods for comparison of growth curves.  Biometrics, 14(1):1–17.
  * Riedel and Sidorenko, [1995] Riedel, K. S. and Sidorenko, A. (1995).  Minimum bias multiple taper spectral estimation.  IEEE Transactions on Signal Processing, 43(1):188–195.
  * Rubín and Panaretos, [2020] Rubín, T. and Panaretos, V. M. (2020).  Sparsely observed functional time series: estimation and prediction.  Electron. J. Statist., 14(1):1137–1210.
  * Shang and Hyndman, [2017] Shang, H. L. and Hyndman, R. J. (2017).  Grouped functional time series forecasting: An application to age-specific mortality rates.  Journal of Computational and Graphical Statistics, 26(2):330–343.
  * Stoehr et al., [2020] Stoehr, C., Aston, J. A. D., and Kirch, C. (2020).  Detecting changes in the covariance structure of functional time series with application to fmri data.  Econometrics and Statistics.
  * Storey, [2007] Storey, J. D. (2007).  The optimal discovery procedure: A new approach to simultaneous significance testing.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 69(3):347–368.
  * Thomson, [1982] Thomson, D. J. (1982).  Spectrum estimation and harmonic analysis.  Proceedings of the IEEE, 70(9):1055-1096.
  * Trujillo et al., [2017] Trujillo, L. T., Stanfield, C. T., and Vela, R. D. (2017).  The effect of electroencephalogram (eeg) reference choice on information-theoretic measures of the complexity and integration of eeg signals.  Frontiers in Neuroscience, 11:425.
  * van Delft et al., [2021] van Delft, A., Characiejus, V., and Dette, H. (2021).  A nonparametric test for stationarity in functional time series.  Statistica Sinica, 31(3).
  * van Delft and Dette, [2021] van Delft, A. and Dette, H. (2021).  A similarity measure for second order properties of non-stationary functional time series with applications to clustering and testing.  Bernoulli, 27(1):469–501.
  * van Delft and Eichler, [2018] van Delft, A. and Eichler, M. (2018).  Locally stationary functional time series.  Electronic Journal of Statistics, 12(1):107–170.
  * Walden et al., [1995] Walden, A. T., McCoy, E. J., and Percival, D. B. (1995).  The effective bandwidth of a multitaper spectral estimator.  Biometrika, 82(1):201–214.
  * Wei et al., [2018] Wei, J., Chen, T., Li, C., Liu, G., Qiu, J., and Wei, D. (2018).  Eyes-open and eyes-closed resting states with opposite brain activity in sensorimotor and occipital regions: Multidimensional evidences from machine learning perspective.  Frontiers in Human Neuroscience, 12:422.