#  An algebraic estimator
for large spectral density matrices 111This work was supported by the
_Accademia Nazionale dei Lincei_ under Grant _British Academy 2018._

Matteo Barigozzi
Department of Economics, University of Bologna
and
Matteo Farnè
Department of Statistics, University of Bologna Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

We propose a new estimator of high-dimensional spectral density matrices,
called UNshrunk ALgebraic Spectral Estimator (UNALSE), under the assumption of
an underlying low rank plus sparse structure, as typically assumed in dynamic
factor models. The UNALSE is computed by minimizing a quadratic loss under a
nuclear norm plus $l_{1}$ norm constraint to control the latent rank and the
residual sparsity pattern. The loss function requires as input the classical
smoothed periodogram estimator and two threshold parameters, the choice of
which is thoroughly discussed. We prove consistency of UNALSE as both the
dimension $p$ and the sample size $T$ diverge to infinity, as well as
algebraic consistency, i.e., the recovery of latent rank and residual sparsity
pattern with probability one. The finite sample properties of UNALSE are
studied by means of an extended simulation exercise as well as an empirical
analysis of US macroeconomic data.

Keywords: Spectral density matrix, High-dimensions, Dynamic rank, Sparsity, US
macroeconomic data.

##  1 Introduction

An appealing, natural, and classical way to model time series data is through
spectral analysis (Brillinger, , 2001). Given a $p$-dimensional vector
stochastic process, its $p\times p$ spectral density matrix characterizes all
second order dependencies. Moreover, conditional second order dependencies can
also be extracted starting from the inverse of the spectral density matrix.
The spectral approach is appealing since, once we move from the time domain to
the frequency domain, data become asymptotically independent, as the sample
size $n$ grows to infinity.

Statistical methods for the study of time series based on spectral analysis
include: pseudo-maximum likelihood estimation (Dahlhaus, 2000b, ; Velasco and
Robinson, , 2000), linear regression (Harvey, , 1978), cointegration tests or
information criteria based on the zero-frequency spectral density matrix of a
vector of time series (Stock and Watson, , 1988; Barigozzi et al., , 2021),
and similarly seasonal cointegration tests based on the spectral density
matrix at selected frequencies (Joyeux, , 1992), de-trending methods (Corbae
et al., , 2002), Granger causality tests (Breitung and Candelon, , 2006; Farnè
and Montanari, , 2021), and the analysis of low frequency co-movements (Müller
and Watson, , 2018). Finally, the inverse spectral density matrix is at the
basis of graphical models and dynamic network analysis (Granger, , 1969;
Dahlhaus, 2000a, ; Eichler, , 2007; Davis et al., , 2016; Barigozzi and
Brownlees, , 2019).

The use of spectral analysis is widespread in many applied fields. Examples
are the construction of business cycle indicators in macroeconomics (Sargent
and Sims, , 1977; Altissimo et al., , 2010), portfolio optimization at
different horizons in finance (Chaudhuri and Lo, , 2015), and the study of
brain activity in biostatistics (Ombao et al., , 2001, 2005; Fiecas and Ombao,
, 2011, 2016).

All above methods and applications require as input an estimator of the
spectral density matrix or of its inverse. Just like for the covariance matrix
estimation in time domain, estimation of a spectral density matrix is a hard
problem when the dimension of the process $p$ is comparable, or even larger,
than the sample size $T$. In this case, the classical smoothed periodogram
estimator is not positive definite simply due to lack of degrees of freedom.
Given the increased availability of large datasets in the recent years, this
issue becomes of fundamental importance. Wu and Zaffaroni,  (2018) provide
consistency results for the smoothed periodogram estimator in high-dimension,
which hold uniformly over all frequencies (see also Zhang and Wu, , 2021).

To solve the problem of the curse of dimensionality, here, we start from
observing that the second moments of most high-dimensional time series tend to
have both a low rank and a sparse component. Indeed, on the one hand, most
economic datasets are known to be mainly “dense” rather than sparse (Giannone
et al., , 2017). Moreover, there exist mathematical results proving that large
dimensional panels of time series can always be represented as having a factor
structure (Forni and Lippi, , 2001; Hallin and Lippi, , 2013). On the other
hand, once the common factors are controlled for, there is evidence of
sparseness in the second order structure of the residuals (Barigozzi and
Hallin, , 2017).

In this paper, we assume that the spectral density matrix of $p$-dimensional
time series has the low rank plus sparse structure:

$\Sigma(\theta)=L(\theta)+S(\theta),\qquad\theta\in[-\pi,\pi],$

where $L(\theta)$ has rank $r$ independent of $p$ and such that $r<p$, and
$S(\theta)$ is a sparse matrix. Based on this assumption, our estimators
$\widehat{L}(\theta)$ and $\widehat{S}(\theta)$ of the two components of the
spectral density matrix are obtained by regularizing the smoothed periodogram
estimator, $\widetilde{\Sigma}(\theta)$, by means of a nuclear norm plus
$l_{1}$ norm penalization. Specifically, at each given frequency our
estimators are defined as

$\left(\widehat{L}(\theta),\widehat{S}(\theta)\right)=\arg\\!\\!\\!\min_{\underline{L}(\theta),\underline{S}(\theta)}\frac{1}{2}\|\widetilde{\Sigma}(\theta)-(\underline{L}(\theta)+\underline{S}(\theta))\|_{F}^{2}+\psi\|\underline{L}(\theta)\|_{*}+\rho\|\underline{S}(\theta)\|_{1},$

where $\underline{L}(\theta)$ and $\underline{S}(\theta)$ indicate generic
values of the matrices belonging to appropriate algebraic matrix varieties,
$\|\underline{L}(\theta)\|_{*}=\text{tr}(\underline{L}(\theta))$ and
$\|\underline{S}(\theta)\|_{1}=\sum_{i,j=1}^{p}|\underline{S}_{ij}(\theta)|$,
and $\psi$ and $\rho$ are threshold parameters.

The above optimization problem is solved by iterating between a singular value
thresholding step (Cai et al., , 2010), giving $\widehat{L}(\theta)$, and a
soft-thresholding step (Daubechies et al., , 2004), giving
$\widehat{S}(\theta)$. The algorithm we employ has also been described in Luo,
2011b for the case of covariance estimation. We also apply the un-shrinkage
step of estimated latent eigenvalues by Farnè and Montanari,  (2020), that
optimizes the finite sample Frobenius loss with respect to the smoothed
periodogram while retaining algebraic consistency. We call the resulting
estimator of the overall spectral density matrix,
$\widehat{\Sigma}(\theta)=\widehat{L}(\theta)+\widehat{S}(\theta)$, UNshrunk
ALgebraic Spectral Estimator (UNALSE).

We prove the algebraic and parametric consistency of UNALSE uniformly over
frequencies, as both the dimension $p$ and the sample size $T$ diverge. By
algebraic consistency, we mean that, with probability tending to $1$: (i) the
UNALSE low rank estimate $\widehat{L}(\theta)$ is positive semidefinite with
the true rank $r$, (ii) the UNALSE residual estimate $\widehat{S}(\theta)$ is
positive definite having the true sparsity pattern, and (iii)
$\widehat{\Sigma}(\theta)$ is positive definite. The usual parametric
consistency holds because UNALSE estimates are close to their targets in
spectral norm (rescaled by the dimension $p$) with probability approaching
$1$. Our consistency results are obtained by generalizing to our framework the
results of Wu and Zaffaroni,  (2018) for the smoothed periodogram. We also
provide a thorough discussion on the selection of the threshold parameters
$\psi$ and $\rho$.

Our approach is based on the fundamental identifiability assumptions we make
on the behavior of the eigenvalues of the spectral density matrix. We assume
the $r$ eigenvalues of the low rank component, $L(\theta)$, to be diverging at
a rate $p^{\alpha}$ with $\alpha\in[0,1]$, possibly different across
frequencies. In the language of factor models, this means we are allowing for
the presence of factors with different degrees of pervasiveness across
frequencies, i.e., both weak and strong factors. Moreover, we assume the
sparse component, $S(\theta)$, to have eigenvalues diverging at a rate
$p^{\delta}$ with $\delta\in[0,1/2]$ and $\delta<\alpha$, possibly different
across frequencies. These assumptions imply the existence of an eigen-gap in
the spectrum of the spectral density matrix, $\Sigma(\theta)$, which has
varying width across frequencies.

There exist alternative approaches to the estimation of large spectral density
matrices. Forni et al.,  (2000) propose principal component analysis in the
frequency domain to recover the low rank component, Böhm and von Sachs,
(2008, 2009) propose to shrink the smoothed periodogram towards either a
reduced rank target or the identity, respectively, Fiecas and von Sachs,
(2014) propose a penalized likelihood approach, and Fiecas et al.,  (2019)
consider constrained $l_{1}$ minimization for estimating the inverse. While
some of those works assume either a low rank or a sparsity structure, none of
them considers both assumptions jointly. Our approach is thus encompassing all
the estimators assuming one of the two settings.

Similar approaches based on a low rank plus sparse assumption exist also in
time domain, i.e., for the estimation of the covariance matrix. Fan et al.,
(2013) consider principal components to recover the low rank component and
then, in a second step, apply soft or hard thresholding to the orthogonal
complement to obtain a sparse and positive definite residual. Their resulting
estimator is called POET. Farnè and Montanari,  (2020) adopt a minimization
algorithm analogous to the one considered in this paper which recovers the
covariance matrix consistently, both algebraically and parametrically. Their
resulting estimator is called UNALCE and they show that it systematically
outperforms POET both in terms of parametric consistency, and, more
importantly, because it provides the algebraic recovery of latent rank and
sparsity pattern. A similar approach was proposed by Luo, 2011b , however it
is based on the assumption of bounded eigenvalues for the covariance matrix,
which does not allow for the joint identification of the two components.

Our assumption of a low rank plus sparse decomposition of the spectral density
matrix is strictly related to, and inspired by, the Generalized Dynamic Factor
Model (GDFM) representation of a large panel of time series, originally proved
by Forni and Lippi,  (2001). This is a very popular approach to dimension
reduction (see, e.g., the application in Altissimo et al., , 2010). In the
GDFM, $r$ latent factors are loaded by each series in a dynamic way, i.e., not
only contemporaneously but also with lags. The key assumptions are: (i)
pervasiveness of the factors modeled via $r$ leading spiking spectral
eigenvalues, and (ii) weak serial and cross-correlation in the residuals,
modeled via boundedness of the spectral eigenvalues. These in our notation
imply $\alpha=1$.

Forni et al.,  (2000, 2005, 2017) consider different estimators of the GDFM,
which are all built starting from a consistent estimator of the spectral
density matrix. In particular, in all those approaches the low rank component
of the spectral density is estimated via the $r$ leading dynamic principal
components, i.e., the principal components of the spectral density matrix
across frequencies of the smoothed periodogram (see also Brillinger, , 2001).
The consistency of this method relies on the pervasiveness of spectral
eigenvalues with respect to the dimension $p$. The spectral density of the
residual component, called idiosyncratic component in the GDFM literature, is
then estimated as the difference between the estimated spectral density of the
observed data and its estimated low rank component. Hence, by construction,
the spectral density of the idiosyncratic component has rank $p-r$, i.e., it
is not positive definite, and, therefore, not invertible. There exist also few
papers dealing with determining the dynamic rank, $r$: Hallin and Liška,
(2007) propose an information criterion, and Onatski,  (2009) proposes a test
based on the asymptotic distribution of the spectral eigenvalues.

The above approaches to the estimation of the GDFM suffer from some drawbacks.
First, any estimator of the spectral density matrix based on the principal
components of an input estimator, like the smoothed periodogram, is likely to
suffer from numerical instability, especially if $p$ is large, due to the
Marčenko and Pastur,  (1967) law. Second, the strict pervasiveness assumption
of spectral eigenvalues ($\alpha=1$) is rarely satisfied in practice, since
both the factor number and their strength might vary across frequencies, e.g.,
due to common, frequency specific, features. Third, the weak correlation
assumption increases the number of parameters when $p$ is large, which
prevents the residual component to be identified.

The estimator we propose in this paper can be used as input of all the
estimators of the GDFM considered in the literature, and, given its algebraic
consistency, it is also a consistent estimator of the latent rank $r$.
Moreover, as already noticed above, our assumptions generalize the GDFM
setting in that they are compatible with frequency specific factor numbers and
strengths, and goes beyond a factor structure in that if no factor is present,
then our method would not return any low rank plus sparse decomposition.

The paper is organized as follows. In Section 2 we present our main results
using the GDFM setting as a guiding example. In Sections 3 and 4 we present
the general framework, describe estimation, and prove consistency. In Section
5 we discuss the choice of the threshold parameters. Sections 6 and 7 present
numerical results for simulated and real datasets.

### Notation

Let us define a $p\times p$ symmetric positive-definite complex matrix ${M}$,
and denote its transposed complex conjugate as $M^{\dagger}$. We denote by
$\lambda_{i}({M})$, $i=1,\ldots,p$, the eigenvalues of ${M}$ in descending
order (note that they are all real numbers), and by $M_{ij}$ the $(i,j)$th
entry of $M$. We also define $\overline{M}_{ij}$ as the complex conjugate of
$M_{ij}$, thus the complex modulus is
$|M_{ij}|=\sqrt{M_{ij}\overline{M}_{ij}}$, while the real and imaginary parts
are indicated as $\text{Re}(M_{ij})$ and $\text{Im}(M_{ij})$, respectively. To
indicate that $M$ is positive definite or semidefinite we use the notations:
${M}\succ 0$ or ${M}\succeq 0$, respectively.

Element-wise norms: $l_{0}$ norm:
$\|M\|_{0}=\sum_{i=1}^{p}\sum_{j=1}^{p}\mathbbm{1}(M_{ij}\neq 0)$, which is
the total number of nonzeros; $l_{1}$ norm:
$\|M\|_{1}=\sum_{i=1}^{p}\sum_{j=1}^{p}|M_{ij}|$; Frobenius norm:
$\|M\|_{F}=\sqrt{\sum_{i=1}^{p}\sum_{j=1}^{p}|M_{ij}|^{2}}=\sqrt{\text{tr}(MM^{\dagger})}$;
maximum norm: $\|M\|_{\infty}=\max_{1\leq i,j\leq p}|M_{ij}|$. Vector-induced
norms: $\|M\|_{0,v}=\max_{1\leq i\leq p}\sum_{j=1}^{p}\mathbbm{1}(M_{ij}\neq
0)$, which is the maximum number of nonzeros per column, defined also as the
maximum “degree” of ${M}$; $\|M\|_{1,v}=\max_{1\leq j\leq
p}\sum_{i=1}^{p}|M_{ij}|$; $\|M\|_{\infty,v}=\max_{1\leq i\leq
p}\sum_{j=1}^{p}|M_{ij}|$; spectral norm:
$\|M\|_{2}=\sqrt{\lambda_{1}({MM^{\dagger}})}=\lambda_{1}({M})$; the nuclear
norm: $\|M\|_{*}=\mbox{tr}(M)=\sum_{i=1}^{p}\lambda_{i}({M})$. The minimum
nonzero off-diagonal element of $M$ in absolute value is denoted as
$\|M\|_{min,off}=\min_{\begin{subarray}{c}1\leq i,j\leq p\\\ i\neq
j,M_{ij}\neq 0\end{subarray}}{|M_{ij}|}$.

##  2 Main results

In this section, we present the main features of our estimator under the
assumption that the data follow a GDFM as defined by Forni and Lippi,  (2001)
and Hallin and Lippi,  (2013). The GDFM setting has to be considered just as a
motivating example, which is well suited to allow the reader to immediately
appreciate the contribution of this paper with respect to the state of art. In
the following sections, we present our theory in more detail showing that the
validity of our results is actually much broader than the case here
considered.

Let $X=\\{X_{it},i=1,\ldots,p\,,t\in\mathbb{Z}\\}$ be a $p$-dimensional panel
of time series. We assume that for any $p\in\mathbb{N}$ the process $X$ is
second-order stationary, and, without loss of generality, we also assume that
$\mathrm{E}[X_{it}]=0$ and $V(X_{it})=\mathrm{E}[{X_{it}^{2}}]>0$, for any
$i\in\mathbb{N}$. The set of all $L_{2}$-convergent linear combinations of
$X_{it}$’s and their limits, as $p\to\infty$, of $L_{2}$-convergent sequences
thereof, is a Hilbert space, denoted by $\mathcal{H}_{X}$. Hence, for all
$t\in\mathbb{Z}$ and all $p\in\mathbb{N}$, any dynamic linear combination of
$X_{it}$s, $y_{t}=\sum_{i=1}^{p}\sum_{k=-\infty}^{\infty}a_{ik}X_{i,t-k}$,
such that $\sum_{i=1}^{p}\sum_{k=-\infty}^{\infty}a_{ik}^{2}=1$, belongs to
$\mathcal{H}_{X}$. Following Definitions 2.1 and 2.2 in Hallin and Lippi,
(2013), we define as _common variable_ the $L_{2}$-limit of any standardized
dynamic linear combination of the $X$s, say $y_{t}/V(y_{t})^{1/2}$, such that
$V(y_{t})\to\infty$, as $p\to\infty$. The Hilbert space of all _common
variables_ is denoted by $\mathcal{H}_{\textit{com}}$, while its orthogonal
complement with respect to $\mathcal{H}_{X}$, denoted as
$\mathcal{H}_{\textit{idio}}$, contains all the _idiosyncratic variables_ ,
i.e., all dynamic linear combinations $y_{t}$ with bounded variance $V(y_{t})$
for all $p\in\mathbb{N}$. Hallin and Lippi,  (2013) prove that there exist two
unique stochastic processes $\chi_{it}\in\mathcal{H}_{\textit{com}}$ and
$\epsilon_{it}\in\mathcal{H}_{\textit{idio}}$, mutually orthogonal at all
leads and lags, such that

$X_{it}=\chi_{it}+\epsilon_{it}$ (1)

for all $i\in\mathbb{N}$ and $t\in\mathbb{Z}$. The process $\chi_{it}$ is
called _common component_ , the process $\epsilon_{it}$ is called
_idiosyncratic component_. Representation (1) is the GDFM. In vector terms, we
can write $X_{t}=\chi_{t}+\epsilon_{t}$, where $X_{t}$, $\chi_{t}$ and
$\epsilon_{t}$ are $p$-dimensional random vectors. The GDFM encompasses the
_approximate static_ factor models of Fan et al.,  (2013), as well as the
_exact dynamic_ factor models of Sargent and Sims,  (1977).

Let us define the spectral density matrices of $\\{\chi_{t}\\}$ and of
$\\{\epsilon_{t}\\}$ as

$L(\theta)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\Gamma_{\chi}(k)\mathrm{e}^{-{\rm
i}\theta k}\;\mbox{ and
}\;S(\theta)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\Gamma_{\epsilon}(k)\mathrm{e}^{-{\rm
i}\theta k},\qquad\theta\in[-\pi,\pi],$

respectively, where $\Gamma_{\chi}(k)=\mathrm{E}[\chi_{t+k}\chi_{t}^{\prime}]$
and $\Gamma_{\epsilon}(k)=\mathrm{E}[\epsilon_{t+k}\epsilon_{t}^{\prime}]$.
Forni and Lippi,  (2001) prove that: (i) the common component $\chi_{it}$ is
driven by a $r$-tuple of mutually orthogonal white noises loaded by a linear
time filter, and $\epsilon_{it}$ is orthogonal to those white noises at all
leads and lags, and (ii) $X_{it}$ follows the GDFM representation (1) if and
only if the $r$ eigenvalues of $L(\theta)$ diverge almost everywhere across
$[-\pi,\pi]$ as $p$ diverges, while the eigenvalues of $S(\theta)$ remain
bounded for all $p$.

Hence, the spectral density matrix of $X_{t}$ is such that
$\Sigma(\theta)=L(\theta)+S(\theta)$, for any $\theta\in[-\pi,\pi]$, and, as
usual in the GDFM literature, in this section we adopt the assumption (relaxed
later on) that the $r$ eigenvalues of $p^{-1}L(\theta)$ are bounded away from
$0$ for all $p$ almost everywhere across the frequency range $[-\pi,\pi]$.
Similarly, the definition of idiosyncratic variable leads to the condition
$\|S(\theta)\|_{2}<\infty$ almost everywhere across the frequency range
$[-\pi,\pi]$ for all $p$. These assumptions on $L(\theta)$ and $S(\theta)$
cause the fact that the gap between the $r$th and the $(r+1)$th eigenvalue of
the spectral density matrix $\Sigma(\theta)=L(\theta)+S(\theta)$ increases at
each $\theta\in[-\pi,\pi]$ as $p$ diverges, making the recovery of the low
rank component easier.

In this paper, we further control the idiosyncratic spectral density matrix
$S(\theta)$ at each $\theta$ by enforcing $\|S(\theta)\|_{0,v}$ to be bounded
and finite for all $p$. Since $\|S(\theta)\|_{2}\leq\|S(\theta)\|_{0,v}$, the
original assumption $\|S(\theta)\|_{2}<\infty$ still holds. This is done in
order to enforce element-wise sparsity on $S(\theta)$ at each $\theta$.

Suppose now that we observe a sample of $p$-dimensional data vectors with size
$T$. A classical estimator of the spectral density matrix, which is our _pre-
estimator_ , is the smoothed periodogram, defined as

$\widetilde{\Sigma}(\theta_{h})=\frac{1}{2\pi}\sum_{k=-(T-1)}^{T-1}K\left(\frac{k}{M_{T}}\right)\mathrm{e}^{-{\rm
i}\theta_{h}k}\,\widehat{\Gamma}_{X}(k),\qquad\theta_{h}=\frac{h\pi}{M_{T}},\quad|h|\leq\lfloor
M_{T}\rfloor,$ (2)

where $\widehat{\Gamma}_{X}(k)=T^{-1}\sum_{t=1}^{T-|k|}X_{t}X_{t+k}^{\prime}$,
and $K(\cdot)$ is a suitable kernel function with $M_{T}$ being the associated
smoothing parameter. According to Brillinger,  (2001), for any given
$\theta_{h}$, $\widetilde{\Sigma}(\theta_{h})$ is consistent if
$\tfrac{M_{T}}{T}\rightarrow 0$ while $M_{T}\rightarrow\infty$ and
$T\rightarrow\infty$. Wu and Zaffaroni,  (2018) prove the consistency of
$\widetilde{\Sigma}(\theta_{h})$ uniformly over the frequencies, under
appropriate assumptions to be discussed later.

Under the GDFM setting described above, augmented with the sparsity assumption
for $S(\theta)$, we define the UNshrunk ALgebraic Spectral Estimator (UNALSE)
estimator of the spectral density matrix $\Sigma(\theta)$ as
$\widehat{\Sigma}(\theta)=\widehat{L}(\theta)+\widehat{S}(\theta)$, where
$\widehat{L}(\theta)$ and $\widehat{S}(\theta)$ are such that:

$\left(\widehat{L}(\theta),\widehat{S}(\theta)\right)=\arg\\!\\!\\!\min_{\underline{L}(\theta),\underline{S}(\theta)}\frac{1}{2}\|\widetilde{\Sigma}(\theta)-(\underline{L}(\theta)+\underline{S}(\theta))\|_{F}^{2}+\psi\|\underline{L}(\theta)\|_{*}+\rho\|\underline{S}(\theta)\|_{1},$
(3)

where $\underline{L}(\theta)$ and $\underline{S}(\theta)$ indicate generic
values of the matrices belonging to appropriate algebraic matrix varieties
(see Section 3 for the details), while $\psi$ and $\rho$ are threshold
parameters. The minimization problem (3) is a non-smooth convex optimization
problem which is the tightest convex relaxation of the following NP-hard
problem:

$\min_{\underline{L}(\theta),\underline{S}(\theta)}\frac{1}{2}\|\widetilde{\Sigma}(\theta)-(\underline{L}(\theta)+\underline{S}(\theta))\|_{F}^{2}+\psi\,\text{{rk}}(\underline{L}(\theta))+\rho\|\underline{S}(\theta)\|_{0},$
(4)

which would be the natural target under the low rank plus sparse assumption.
Indeed, we know that:  (i) $\|\underline{S}(\theta)\|_{1}$is the tightest
convex relaxation of $\|\underline{S}(\theta)\|_{0}$ (Donoho, , 2006);  (ii)
$\|\underline{L}(\theta)\|_{*}$is the tightest convex relaxation of
$\text{rk}(\underline{L}(\theta))$ (Fazel et al., , 2001).

In practice, the solution of (3) is computed as follows. For any given
frequency $\theta_{h}=\frac{\pi h}{M_{T}}$, with $|h|\leq\lfloor
M_{T}\rfloor$, we apply the following iterative procedure:

  1. 1.

set $\left(L_{0}(\theta_{h}),S_{0}(\theta_{h})\right)=\left(\frac{\text{\it
diag}(\widetilde{\Sigma}(\theta))}{2},\frac{\text{\it
diag}(\widetilde{\Sigma}(\theta))}{2}\right)$, $\eta_{0}=1$, and initialize
$Y_{0}(\theta_{h})=L_{0}(\theta_{h})$ and
$Z_{0}(\theta_{h})=S_{0}(\theta_{h})$;

  2. 2.

for $k\geq 0$, repeat:

    1. (a)

compute
$\frac{\partial\frac{1}{2}\|Y_{k-1}(\theta_{h})+Z_{k-1}(\theta_{h})-\widetilde{\Sigma}(\theta_{h})\|_{F}^{2}}{\partial
Y_{k-1}(\theta_{h})}=\frac{\partial\frac{1}{2}\|Y_{k-1}(\theta_{h})+Z_{k-1}(\theta_{h})-\widetilde{\Sigma}(\theta_{h})\|_{F}^{2}}{\partial
Z_{k-1}(\theta_{h})}=Y_{k-1}(\theta_{h})+Z_{k-1}(\theta_{h})-\widetilde{\Sigma}(\theta_{h})$;

    2. (b)

apply the _singular value thresholding_ operator of Cai et al.,  (2010)
$T_{\psi}(\cdot)$ to
$\mathcal{E}_{Y,k}(\theta_{h})=Y_{k-1}(\theta_{h})-\frac{1}{2}(Y_{k-1}(\theta_{h})+Z_{k-1}(\theta_{h})-\widetilde{\Sigma}(\theta_{h}))$
and set $L_{k}(\theta_{h})=T_{\psi}(\mathcal{E}_{Y,k}(\theta_{h}))$;

    3. (c)

apply the _soft-thresholding_ operator of Daubechies et al.,  (2004)
$T_{\rho}(\cdot)$ to
$\mathcal{E}_{Z,k}(\theta_{h})=Z_{k-1}(\theta_{h})-\frac{1}{2}(Y_{k-1}(\theta_{h})+Z_{k-1}(\theta_{h})-\widetilde{\Sigma}(\theta_{h}))$
and set $S_{k}(\theta_{h})=T_{\rho}({\mathcal{E}_{Z,k}(\theta_{h})})$;

    4. (d)

set
$(Y_{k}(\theta_{h}),Z_{k}(\theta_{h}))=(L_{k}(\theta_{h}),S_{k}(\theta_{h}))+\frac{\eta_{k-1}-1}{\eta_{k}}[(L_{k}(\theta_{h}),S_{k}(\theta_{h}))-(L_{k-1}(\theta_{h}),S_{k-1}(\theta_{h}))]$
where $\eta_{k}=\frac{1+\sqrt{1+4\eta_{k-1}^{2}}}{2}$;

    5. (e)

stop if
$\frac{\|L_{k}-L_{k-1}\|_{F}}{\|1+L_{k-1}\|_{F}}+\frac{\|S_{k}-S_{k-1}\|_{F}}{\|1+S_{k-1}\|_{F}}\leq\varsigma$,
where $\varsigma$ is a prescribed precision level (we set $\varsigma=0.01$);

  3. 3.

set $\widehat{L}(\theta_{h})=Y_{k}(\theta_{h})$ and
$\widehat{S}(\theta_{h})=Z_{k}(\theta_{h})$.

The two thresholding operators introduced in the above algorithm are defined
as follows.

  1. (I)

Singular value thresholding operator: let the singular value decomposition of
a positive semi-definite complex symmetric matrix $M$ be
$M=U_{M}\Lambda_{M}U_{M}^{\dagger}$, then, define
$T_{\psi}(M)$=$U_{M}\Lambda_{\psi,M}U_{M}^{\dagger}$, where $\Lambda_{\psi,M}$
is a diagonal matrix with $i$th diagonal element
$\Lambda_{\psi,ii,M}=\max{(\Lambda_{ii,M}-\psi,0)}$.

  2. (II)

Soft-thresholding operator: for a positive definite complex symmetric $M$
define$T_{\rho}(M_{ij})=\frac{M_{ij}}{\sqrt{M_{ij}\overline{M}_{ij}}}\max\left(\sqrt{M_{ij}\overline{M}_{ij}}-\rho,0\right)$.

In this paper we prove the following results for the UNALSE estimator.

######  Theorem 2.1.

For all $p\in\mathbb{N}$, assume that: (i) the $r$ nonzero eigenvalues of
$L(\theta)$ are such that $\frac{\lambda_{j}(L(\theta))}{p}$ is finite and
bounded away from zero for all $j=1,\ldots,r$, and (ii) $\|S(\theta)\|_{0,v}$
is bounded. Then, under the regularity conditions in Section 4, there exist
finite positive reals $G_{1}$, $G_{2}$, $G_{3}$, $G_{4}$, and $G_{5}$,
independent of $p$ and $T$, such that, as $p,T\to\infty$, with probability
approaching $1$, for $\theta_{h}=\tfrac{h\pi}{M_{T}}$:

  1. 1.

$\text{{rk}}(\widehat{L}(\theta_{h}))=\text{{rk}}({L}(\theta_{h}))=r$ and
$\max_{|h|\leq\lfloor
M_{T}\rfloor}\frac{1}{p}\|\widehat{L}(\theta_{h})-{L}(\theta_{h})\|_{2}\leq
G_{1}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$;

  2. 2.

$\max_{|h|\leq\lfloor
M_{T}\rfloor}\frac{1}{p}\|\widehat{S}(\theta_{h})-{S}(\theta_{h})\|_{\infty}\leq
G_{2}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$, and, consequently,
$\max_{|h|\leq\lfloor
M_{T}\rfloor}\frac{1}{p}\|\widehat{S}(\theta_{h})-{S}(\theta_{h})\|_{2}\leq
G_{2}\|{S}(\theta_{h})\|_{0,v}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$;

  3. 3.

$\max_{|h|\leq\lfloor
M_{T}\rfloor}\frac{1}{p}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}\leq
G_{3}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$.

Furthermore, if $\frac{\lambda_{p}(S(\theta))}{p}\geq
2G_{4}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$ and
$\frac{\lambda_{p}(\Sigma(\theta))}{p}\geq
2G_{5}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$, then, as $p,T\to\infty$, with
probability approaching $1$, for $\theta_{h}=\tfrac{h\pi}{M_{T}}$:

  1. 4.

$\widehat{S}(\theta_{h})$ is positive definite and $\max_{|h|\leq\lfloor
M_{T}\rfloor}\frac{1}{p}\|\widehat{S}(\theta_{h})^{-1}\leavevmode\nobreak\
-\leavevmode\nobreak\ S(\theta_{h})^{-1}\|_{2}\leq
G_{4}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$;

  2. 5.

$\widehat{\Sigma}(\theta_{h})$ is positive definite and $\max_{|h|\leq\lfloor
M_{T}\rfloor}\frac{1}{p}\|\widehat{\Sigma}(\theta_{h})^{-1}-\Sigma(\theta_{h})^{-1}\|_{2}\leq
G_{5}\sqrt{\frac{M_{T}\log(M_{T})}{T}}$.

The convergence speed depends on the rapidity of decay of the physical
dependence, on the finite moment of highest order, and on the smoothing
parameter $M_{T}$ (see Section 4 for further details and a more detailed
exposition). We stress that the standard condition $\frac{M_{T}}{T}\rightarrow
0$ (as $M_{T},T\rightarrow\infty$) ensures that the relative errors in the
above Theorem vanish asymptotically uniformly over the frequency grid.
Moreover, if $p\sqrt{\frac{M_{T}}{T}}\to 0$, the conditions for parts 4 and 5
to hold are certainly satisfied as $p$ diverges.

The results of Theorem 2.1 contribute to the literature in three ways. First,
the exact dynamic rank recovery in part 1 allows to bypass the use of existing
criteria for determining the number of factors, like those by Hallin and
Liška,  (2007) and Onatski,  (2009). Second, assuming that $S(\theta)$ is full
rank, we derive a consistency result also for the estimator of the
idiosyncratic spectral density $\widehat{S}(\theta)$, which allows to obtain
the same error bound also for the overall estimator
$\widehat{\Sigma}(\theta)$. Third, we obtain results also for the estimators
of the inverse spectral densities.

The validity of our estimation framework lies well beyond the standard GDFM
assumptions. First, we can relax the strict pervasiveness assumption on latent
dynamic factors, by allowing the $r$ eigenvalues of the matrix
$p^{-\alpha}L(\theta)$, with $\alpha\leq 1$, to be bounded away from $0$
almost everywhere across the frequency range $[-\pi,\pi]$. Second, we allow
for the maximum number of nonzeros per row in $S(\theta)$,
$\|S(\theta)\|_{0,v}$, to be at most proportional to $p^{\delta}$, with
$\delta\in[0,1/2]$ and $\delta<\alpha$. This means that we allow the
idiosyncratic spectrum to be quite far from the diagonal matrix. Our setting
reduces to the GDFM one when $\alpha=1$. Third, the rank $r$ can be in fact
frequency dependent.

##  3 Model setup

The aim of this paper is estimating the spectral density matrix of a
$p$-dimensional process $X=\\{X_{it},i=1,\ldots,p,\,t\in\mathbb{Z}\\}$. We
consider the following data generating process for $X$:

$\displaystyle X_{t}$
$\displaystyle=\chi_{t}+\epsilon_{t},\quad\quad\quad\;\;\,t\in\mathbb{Z},$ (5)
$\displaystyle\chi_{t}$ $\displaystyle=\sum_{s=0}^{\infty}B_{s}u_{t-s},\qquad
t\in\mathbb{Z},$ (6) $\displaystyle\epsilon_{t}$
$\displaystyle=\sum_{s=0}^{\infty}C_{s}e_{t-s},\qquad t\in\mathbb{Z},$ (7)

where $X_{t}$, $\chi_{t}$, $\epsilon_{t}$, and $e_{t}$ are $p$-dimensional,
$u_{t}$ is $r$-dimensional, the $B_{s}$ are $p\times r$, and the $C_{s}$ are
$p\times p$. Note that, differently from the original works on GDFM by Forni
et al.,  (2000) and Forni and Lippi,  (2001) who derive (6) using two-sided
filters, here we follow more the recent works by Hallin and Lippi,  (2013) and
Forni et al.,  (2017) who show that also a one-sided representation is
possible. We make the following assumptions on the processes $\\{u_{t}\\}$ and
$\\{e_{t}\\}$ in (6) and (7).

######  Assumption 3.1.

  1. (i)

$\\{u_{t},t\in\mathbb{Z}\\}$ is a $r$-dimensional independent process with
$r<p$ and independent of $p$, $\mathrm{E}(u_{t})={0}_{r}$ and
$\mathrm{E}[u_{t}u_{t}^{\prime}]={I_{r}}$;

  2. (ii)

there exists a $K_{u}$ independent of $j$ such that
$\mathrm{E}[|u_{jt}|^{4+\epsilon}]\leq K_{u}$ for some $\epsilon>0$ and
$j=1,\ldots,r$;

  3. (iii)

$\\{e_{t},t\in\mathbb{Z}\\}$ is a $p$-dimensional independent process with
$\mathrm{E}(e_{t})={0}_{p}$ and $\mathrm{E}[e_{t}e_{t}^{\prime}]={I_{p}}$;

  4. (iv)

there exists a $K_{e}$ independent of $j$ such that
$\mathrm{E}[|e_{jt}|^{4+\epsilon}]\leq K_{e}$ for some $\epsilon>0$ and
$j=1,\ldots,p$;

  5. (v)

$\\{u_{t}\\}$ and $\\{e_{t}\\}$ are two mutually independent processes.

Under Assumption 3.1, processes (6) and (7) are zero-mean linear and weakly
stationary, and consequently process (5) also is.

We define for $z\in\mathbb{C}$ the filter of the common component as
$B(z)=\sum_{s=0}^{\infty}B_{s}z^{s}$, with $B_{s}$ being $p\times r$, and the
filter of the idiosyncratic component as $C(z)=\sum_{s=0}^{\infty}C_{s}z^{s}$,
with $C_{s}$ being $p\times p$. We develop our theory under two different
settings, imposing a basic linear shape for the filters, where all matrices
$B_{s}^{\prime}B_{s}$ have the same condition number and all matrices
$C_{s}C_{s}^{\prime}$ have the same sparsity pattern, and a general linear
shape, where both features are allowed to vary across frequencies.

We start from the first setting.

######  Assumption 3.2 (Basic linear filters).

(i) $B_{s}=U_{L}\Lambda_{L,s}$, where $U_{L}$ is a $p\times r$ matrix such
that $U_{L}^{\prime}U_{L}=I_{r}$, and
$\Lambda_{L,s}=\lambda_{s}\sqrt{\Lambda_{u}}$ with $\lambda_{s}\in\mathbb{R}$
and such that $\sum_{s=0}^{\infty}\lambda_{s}^{2}=1$ and with $\Lambda_{u}$ a
diagonal $r\times r$ positive definite matrix;  (ii)
$C_{s}=U_{S}\Lambda_{S,s}$, where $U_{S}$ is a $p\times p$ matrix such that
$U_{S}^{\prime}U_{S}=I_{p}$ and $\|U_{S}U_{S}^{\prime}\|_{0}=q$ with $q\ll
p^{2}$, and $\Lambda_{S,s}=\lambda_{s}\sqrt{\Lambda_{e}}$ with
$\lambda_{s}\in\mathbb{R}$ and such that
$\sum_{s=0}^{\infty}\lambda_{s}^{2}=1$ and with $\Lambda_{e}$ a diagonal
$p\times p$ positive definite matrix.

It immediately follows that
$\mathrm{E}[\chi_{t}\chi_{t}^{\prime}]=\Gamma_{\chi}(0)=U_{L}\Lambda_{u}U_{L}^{\prime}$
has rank $r$. At the same time,
$\mathrm{E}[\epsilon_{t}\epsilon_{t}^{\prime}]=\Gamma_{\epsilon}(0)=U_{S}\Lambda_{e}U_{S}^{\prime}$
is sparse with $q$ nonzero elements. Notice that in this simple setting the
scalar $\lambda_{s}$ is the same for both filters. Moreover, since the filters
are linear, for all $k\in\mathbb{Z}$, we have:

$\displaystyle\mathrm{E}[\chi_{t}\chi_{t+k}^{\prime}]$
$\displaystyle=\Gamma_{\chi}(k)=\sum_{s=0}^{\infty}B_{s}B_{s+k}^{\prime}=U_{L}\left(\sum_{s=0}^{\infty}\lambda_{s}\lambda_{s+k}\right)\Lambda_{u}U_{L}^{\prime}=a_{k}\Gamma_{\chi}(0),$
$\displaystyle\mathrm{E}[\epsilon_{t}\epsilon_{t+k}^{\prime}]$
$\displaystyle=\Gamma_{\epsilon}(k)=\sum_{s=0}^{\infty}C_{s}C_{s+k}^{\prime}=U_{S}\left(\sum_{s=0}^{\infty}\lambda_{s}\lambda_{s+k}\right)\Lambda_{e}U_{S}^{\prime}=a_{k}\Gamma_{\epsilon}(0),$

where $a_{k}=\sum_{s=0}^{\infty}\lambda_{s}\lambda_{s+k}$, with $a_{-k}=a_{k}$
and $a_{0}=1$. Obviously, $|a_{k}|\leq a_{0}$, because
$\|\Gamma_{\chi}(k)\|_{2}\leq\|\Gamma_{\chi}(0)\|_{2}$, for all $k$. Since, by
Assumption 3.1, $\\{u_{t}\\}$ and $\\{e_{t}\\}$ are uncorrelated processes, we
obtain
$\mathrm{E}[X_{t}X_{t+k}^{\prime}]=\Gamma_{X}(k)=\Gamma_{\chi}(k)+\Gamma_{\epsilon}(k)=a_{k}(\Gamma_{\chi}(0)+\Gamma_{\epsilon}(0))$
for all $k$, which has a low rank plus sparse structure.

If we define
$a(\theta)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}a_{k}\mathrm{e}^{-\rm{i}\theta
k}$, for $\theta\in[-\pi,\pi]$, then the spectral density matrices of
$\\{\chi_{t}\\}$ and $\\{\epsilon_{t}\\}$ are
$L(\theta)=a(\theta)\Gamma_{\chi}(0)$ and
$S(\theta)=a(\theta)\Gamma_{\epsilon}(0)$ respectively. Therefore, the
spectral density matrix of $\\{X_{t}\\}$, which is
$\Sigma(\theta)=L(\theta)+S(\theta)$, has a low rank plus sparse structure at
all frequencies $\theta\in[-\pi,\pi]$. Note that Assumption 3.2 describes the
simplified case where the spectrum has no phase component, i.e., it is a real
matrix, since $a(\theta)$ is real because $\mathrm{e}^{-\rm{i}\theta
k}+\mathrm{e}^{i\theta k}$ is the sum of two waves of opposite argument.

Alternatively, we assume a more general structure for the filters.

######  Assumption 3.3 (Generalized linear filters).

(i) $B_{s}=U_{L}\Lambda_{L,s}$, where $U_{L}$ is a $p\times r$ matrix such
that $U_{L}^{\prime}U_{L}=I_{r}$, and $\sum_{s=0}^{\infty}{\Lambda_{L,s}^{2}}$
is a diagonal positive definite $r\times r$ matrix,
$\|\sum_{s=0}^{\infty}B_{s}\|_{2}\leq M_{B}(p)$;  (ii)
$C_{s}=U_{S,s}\Lambda_{S,s}$, where $U_{S,s}$ is a $p\times p$ matrix such
that $U_{S,s}^{\prime}U_{S,s}=I_{r}$ and
$\|U_{S,s}U_{S,s}^{\prime}\|_{0}=q_{s}$ with $q_{s}\ll p^{2}$, and
$\sum_{s=0}^{\infty}{\Lambda_{S,s}^{2}}$ is a diagonal positive definite
$p\times p$ matrix, $\|\sum_{s=0}^{\infty}C_{s}\|_{2}\leq M_{C}(p)$.

We refer to Assumption 3.8 below for the definition of the terms $M_{B}(p)$
and $M_{C}(p)$, which for a fixed $p$ are positive constants, but might
diverge as $p$ diverges.

Under Assumption 3.3, for all $k\in\mathbb{Z}$, we have
$\mathrm{E}[\chi_{t}\chi_{t+k}^{\prime}]=\Gamma_{\chi}(k)=\sum_{s=0}^{\infty}B_{s}B_{s+k}^{\prime}=\sum_{s=0}^{\infty}U_{L}\Lambda_{L,s}\Lambda_{L,s+k}U_{L}^{\prime}$,
and
$\mathrm{E}[\epsilon_{t}\epsilon_{t+k}^{\prime}]=\Gamma_{\epsilon}(k)=\sum_{s=0}^{\infty}C_{s}C_{s+k}^{\prime}=\sum_{s=0}^{\infty}U_{S,s}\Lambda_{S,s}\Lambda_{S,s+k}U_{S,s+k}^{\prime}$,
thus leading to
$\Gamma_{\chi}(0)=U_{L}(\sum_{s=0}^{\infty}\Lambda_{L,s}^{2})U_{L}^{\prime}$,
which has rank $r$ by assumption, and
$\Gamma_{\epsilon}(0)=\sum_{s=0}^{\infty}U_{S,s}\Lambda_{S,s}^{2}U_{S,s}^{\prime}$,
which is sparse with $q_{\Gamma_{\epsilon}(0)}=\sum_{s=0}^{\infty}q_{s}$
nonzero elements, while $\Gamma_{\epsilon}(k)$ has $q_{\Gamma_{\epsilon}(k)}$
nonzero elements.

Finally, the spectral density matrix of $\\{\chi_{t}\\}$ is:
$L(\theta)=\frac{1}{2\pi}U_{L}\left(\sum_{k=-\infty}^{\infty}\left(\sum_{s=0}^{\infty}\Lambda_{L,s}\Lambda_{L,s+k}\right)\mathrm{e}^{-\rm{i}\theta
k}\right)U_{L}^{\prime}=U_{L}A(\theta)U_{L}^{\prime}$, where
$A(\theta)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}(\sum_{s=0}^{\infty}\Lambda_{L,s}\Lambda_{L,s+k})\mathrm{e}^{-\rm{i}\theta
k}$. Therefore, $L(\theta)$ has rank $r$ and has the same orthogonal base as
$\Gamma_{\chi}(0)$ at all frequencies $\theta\in[-\pi,\pi]$. Although at first
sight this might be a restrictive design we notice that all results in the
next section hold locally in the algebraic variety of rank $r$ matrices (see
Chandrasekaran et al.,  (2012) and Remark 4.4, below). At the same time, we
obtain
$S(\theta)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\left(\sum_{s=0}^{\infty}U_{S,s}\Lambda_{S,s}\Lambda_{S,s+k}^{\prime}U_{S,s+k}^{\prime}\mathrm{e}^{-\rm{i}\theta
k}\right)$, which has $q$ nonzero elements, where $q\leq
q_{\Gamma_{\epsilon}(0)}+2\sum_{k=1}^{\infty}q_{\Gamma_{\epsilon}(k)}$.
Therefore, the spectral density matrix of $\\{X_{t}\\}$, which is
$\Sigma(\theta)=L(\theta)+S(\theta)$, has a low rank plus sparse structure at
all frequencies $\theta\in[-\pi,\pi]$.

As a consequence of Assumptions 3.1 and 3.2 for the basic filter
specification, or of Assumptions 3.1 and 3.3 for the general filter
specification, the spectral density matrix has a low rank plus sparse
structure. To make this formal we introduce the following algebraic matrix
varieties:

$\displaystyle\mathcal{L}(\mathsf{r})$
$\displaystyle=\\{L\in\mathbb{R}^{p\times p}\mid{L}\succeq
0,L=UDU^{\dagger},U\in\mathbb{C}^{p\times\mathsf{r}},U^{\dagger}U=I_{\mathsf{r}},D\in\mathbb{R}^{\mathsf{r}\times\mathsf{r}}\text{diagonal},\mathsf{r}<p\\},$
(8) $\displaystyle\mathcal{S}(\mathsf{q})$
$\displaystyle=\\{S\in\mathbb{C}^{p\times p}\mid{S}\succeq
0,|\text{supp}(S)|\leq\mathsf{q},\mathsf{q}<p^{2}\\}.$ (9)

In other words, $\mathcal{L}(\mathsf{r})$ is the variety of Hermitian matrices
with at most rank $\mathsf{r}$ and $\mathcal{S}(\mathsf{q})$ is the variety of
Hermitian complex sparse matrices with at most $\mathsf{q}$ nonzero elements
($\text{supp}(S)$ is the orthogonal complement of $\text{{ker}}(S)$).
Therefore, under our assumptions $L(\theta)\in\mathcal{L}(r)$ and
$S(\theta)\in\mathcal{S}(q)$, for all $\theta\in[-\pi,\pi]$, and where $r$ is
defined in Assumptions 3.2(i) or 3.3(i), and $q$ is defined in Assumption
3.2(ii) or under Assumption 3.3(ii).

######  Remark 3.1 (Rank and sparsity).

Notice that while our model assumptions, which combine Assumptions 3.1 and 3.2
for the basic filter specification, as well as Assumptions 3.1 and 3.3 for the
general filter specification, prescribe a fixed and common rank $r$ and
sparsity degree $q$ for the spectral density matrices across frequencies, our
methodology is still able to consistently recover latent structures with
frequency-varying rank and sparsity. However, we avoid this distinction here
to keep the notation simple.

For simplicity, hereafter we adopt the following notation
$L^{*}=\Gamma_{\chi}(0)$ and $S^{*}=\Gamma_{\epsilon}(0)$. Under Assumption
3.2, $L(\theta)=a(\theta)L^{*}$, which has eigenvalues
$\lambda_{r}(L(\theta))=\lambda_{r}(L^{*})a(\theta)$. Since
$L^{*}=U_{L}\Lambda_{u}U_{L}^{\prime}$, in order to characterize the behavior
of $\lambda_{r}(L(\theta))$ we just need to focus on the diagonal elements of
$\Lambda_{u}$ (notice that $a(\theta)$ is a positive scalar).

######  Assumption 3.4.

For all $p\in\mathbb{N}$, the elements of the $r\times r$ diagonal matrix
$\Lambda_{u}$, $\Lambda_{u,j}$, are such that for all $j=1,\ldots,r$,
$M^{\min}_{j}\leq\frac{\Lambda_{u,j}}{p^{\alpha}}\leq M^{\max}_{j}$, with
$0\leq\alpha\leq 1$, and $M^{\min}_{j}$ and $M^{\max}_{j}$ independent of $p$
and $T$ such that $M^{\min}_{j}\geq M^{\max}_{j+1}$, for $j=1,\ldots,(r-1)$.
In addition, there exists some $\kappa_{a}>0$ such that
$a(\theta)\leq\kappa_{a}$ for all $\theta\in[-\pi,\pi]$. Moreover, $r$ is
finite for all $p\in\mathbb{N}$.

Under Assumption 3.3, $L(\theta)=U_{L}A(\theta)U_{L}^{\prime}$, we
characterize the eigenvalues of $L(\theta)$ by means of the following
assumption (notice that $A(\theta)$ is diagonal with positive entries).

######  Assumption 3.5.

For all $p\in\mathbb{N}$ and all $\theta\in[-\pi,\pi]$, the elements,
$A_{j}(\theta)$, of the $r\times r$ diagonal matrix $A(\theta)$, are such that
for all $j=1,\ldots,r$, $M^{\min}_{j}\leq\frac{A_{j}(\theta)}{p^{\alpha}}\leq
M^{\max}_{j}$, with $0\leq\alpha\leq 1$, and $M^{\min}_{j}$ and $M^{\max}_{j}$
independent of $p$ and $T$ such that $M^{\min}_{j}\geq M^{\max}_{j+1}$, for
$j=1,\ldots,(r-1)$. Moreover, $r$ is finite for all $p\in\mathbb{N}$.

In other words, under both assumptions, we are assuming that the eigenvalues
of $L(\theta)$, which are all real, are of order $p^{\alpha}$,
$\alpha\in[0,1]$. In this way we generalize the strict pervasiveness of the
latent factors, corresponding to the case $\alpha=1$, necessary to ensure the
effectiveness of the recovery of the number of factors in the methods by
Hallin and Liška,  (2007) and Onatski,  (2009). We notice that our results
would hold even for $r=O(\log(p))$, however, as common in the literature, and
as proved in Forni and Lippi,  (2001) for the GDFM, we consider $r$ as fixed
and independent of $p$ to avoid unnecessary complications.

In order to give a rigorous definition of sparsity, we need to introduce
further notation. The tangent spaces to $\mathcal{L}(r)$ and $\mathcal{S}(q)$
in two generic matrices $L\in\mathcal{L}(r)$ and $S\in\mathcal{S}(q)$, defined
in (8) and (9), are respectively defined as:

$\displaystyle T(L)$ $\displaystyle=\\{M\in\mathbb{C}^{p\times p}\mid
M=UY_{1}^{\dagger}+Y_{2}U^{\dagger}\mid Y_{1},Y_{2}\in\mathbb{C}^{p\times
r},U^{\dagger}LU\in\mathbb{C}^{r\times
r}\mbox{diagonal},L\in\mathcal{L}(r)\\},$ (10) $\displaystyle\Omega(S)$
$\displaystyle=\\{N\in\mathbb{C}^{p\times p}\mid\text{
supp}(N)\subseteq\text{supp}(S),S\in\mathcal{S}(q)\\}.$ (11)

The following uncertainty principle holds (Chandrasekaran et al., , 2011): if
$L(\theta)$ is nearly sparse, $S(\theta)$ cannot be recovered, and if
$S(\theta)$ is nearly low rank, $L(\theta)$ cannot be recovered. Therefore, in
order to achieve consistency we need to control for the spikiness of the
eigenvalues of $S(\theta)$ and the sparsity pattern of $L(\theta)$. To this
end, we make use of the following rank-sparsity incoherence measures between
$T(L)$ and $\Omega(S)$, introduced in Chandrasekaran et al.,  (2011):
$\xi(T(L))=\max_{\begin{subarray}{c}N\in T(L)\\\ \|N\|_{2}\leq
1\end{subarray}}{\|N\|_{\infty}}$,
$\mu(\Omega(S))=\max_{\begin{subarray}{c}N\in\Omega(S)\\\ \|N\|_{\infty}\leq
1\end{subarray}}{\|N\|_{2}}$. In order to identify $L(\theta)$ and
$S(\theta)$, we need to control these rank-sparsity incoherence measures,
which in turn implies defining the admissible sparsity patterns for
$S(\theta)$. Indeed, as proved in Section 4, a necessary condition to ensure
both parametric and algebraic consistency is

$\xi(T({L}(\theta)))\
\mu(\Omega({S}(\theta))\leq\frac{1}{54},\quad\theta\in[-\pi,\pi],$ (12)

which guarantees that ${L}(\theta)$ is far from sparsity and $S(\theta)$ is
far from rank-deficiency. Indeed, the smaller is the product between the dual
norms $\xi(T({L}(\theta)))$ and $\mu(\Omega({S}(\theta))$, the closer the two
spaces $\mathcal{L}(r)$ and $\mathcal{S}(s)$ are to orthogonality, thus making
easier to perform low rank plus sparse decomposition effectively.

Under Assumption 3.2, $T(L(\theta))=T(L^{*})$ and
$\Omega(S(\theta))=\Omega(S^{*})$, for all $\theta\in[-\pi,\pi]$. We then make
the following assumption.

######  Assumption 3.6.

For all $p\in\mathbb{N}$, there exist:  (i) $\delta\in[0,\frac{1}{2}]$, with
$\delta<\alpha$, and $\delta_{2}>0$, such that $\|S^{*}\|_{0,v}=\max_{1\leq
i\leq p}\sum_{j=1}^{p}\mathbbm{1}(S^{*}_{ij}=0)\leq\delta_{2}p^{\delta}$;
(ii) $\kappa_{L},\kappa_{S}>0$with
$\frac{\kappa_{S}}{\kappa_{L}}\leq\frac{1}{54}$ and
$\kappa_{S}\leq\delta_{2}$, such that
$\xi(T({L^{*}}))=\frac{\sqrt{r}}{\kappa_{L}p^{\delta}}$ and
$\mu(\Omega({S^{*}}))=\kappa_{S}p^{\delta}$;  (iii)
$\underline{\delta},\delta_{2}^{\prime}$with
$\underline{\delta}\leq\delta+\frac{1}{2}$ such that
$\|S^{*}\|_{1,v}\leq\delta_{2}^{\prime}p^{\underline{\delta}}$.

Let us now consider the generalized linear filter setting of Assumption 3.3.
In this case, while, as before, $\xi(T({L(\theta)}))=\xi(T({L^{*}}))$, since
the low rank variety $\mathcal{L}(r)$ does not change,
$\mu(\Omega(S(\theta)))$ is different from $\mu(\Omega(S^{*}))$ because the
underlying sparsity pattern now depends on the sparsity pattern of the matrix
$S_{\infty}=\sum_{k=-\infty}^{\infty}\Gamma_{\epsilon}(k)$. We therefore make
the following assumption.

######  Assumption 3.7.

For all $p\in\mathbb{N}$, there exist:  (i) $\delta\in[0,\frac{1}{2}]$, with
$\delta<\alpha$, and $\delta_{2}>0$, such that
$\|S_{\infty}\|_{0,v}=\max_{1\leq i\leq
p}\sum_{j=1}^{p}\mathbbm{1}(S_{\infty,ij}=0)\leq\delta_{2}p^{\delta}$;  (ii)
$\kappa_{L},\kappa_{S}>0$with $\frac{\kappa_{S}}{\kappa_{L}}\leq\frac{1}{54}$
and $\kappa_{S}\leq\delta_{2}$, such that
$\xi(T({L^{*}}))=\frac{\sqrt{r}}{\kappa_{L}p^{\delta}}$ and
$\mu(\Omega({S_{\infty}}))=\kappa_{S}p^{\delta}$;  (iii)
$\underline{\delta},\delta_{2}^{\prime}$with
$\underline{\delta}\leq\delta+\frac{1}{2}$ such that
$\|S_{\infty}\|_{1,v}\leq\delta_{2}^{\prime}p^{\underline{\delta}}$.

In both Assumptions 3.6 and 3.7, part (ii) guarantees that (12) is satisfied,
and the maximum number of nonzeros per row in $S^{*}$ or in $S_{\infty}$,
which is controlled in part (i), is crucial to this end, as it is a natural
upper bound for $\mu(\Omega({S}(\theta))$ (Chandrasekaran et al., , 2011). The
condition $\delta<\alpha$ is instead needed to preserve the identification of
the two components of the spectral density matrix. Part (iii) controls the
$l_{1,v}$ norms of $S^{*}$ or $S_{\infty}$, notice that this condition is
compatible with part (i).

Hereafter, depending on which filter setting we are considering, we define
$q^{\prime}=\|S^{*}\|_{0,v}$ or $q^{\prime}=\|S_{\infty}\|_{0,v}$. Notice that
by means of Assumptions 3.4 or 3.5 we are constraining the number of nonzeros
in $S(\theta)$ to be at most $q$, while by means of Assumption 3.6(i) or
3.7(i) we are further constraining the maximum number of nonzeros in each row
to be at most $q^{\prime}$. Obviously $q^{\prime}\leq q$ and moreover, because
of our assumptions $q=O(p^{1+\delta})$, hence it can grow with $p$ at a
maximum rate $p^{3/2}$.

Finally, by means of the next assumption and similarly to Wu and Zaffaroni,
(2018), we control the amount of physical dependence of our stochastic process
(5) across time.

######  Assumption 3.8.

There exist $\delta_{3}>0$, $\rho\in[0,1)$, and
$\delta^{\prime}\in[0,\alpha)$, such that, for all $p\in\mathbb{N}$:  (i)
$\|\sum_{s=0}^{\infty}B_{s}\|_{1,v}\leq\delta_{3}p^{\alpha/2}\frac{1}{1-\rho}$;
(ii) $\|\sum_{s=0}^{\infty}B_{s}\|_{\infty,v}\leq\delta_{3}\frac{1}{1-\rho}$;
(iii)
$\|\sum_{s=0}^{\infty}C_{s}\|_{1}\leq\delta_{3}p^{\delta^{\prime}/2}\frac{1}{1-\rho}$;
(iv) $\|\sum_{s=0}^{\infty}C_{s}\|_{\infty,v}\leq\delta_{3}\frac{1}{1-\rho}$.

For any fixed $p$, Assumption 3.8 implies geometrically decaying auto-
covariances as $\rho^{|s|}$. Similar assumptions are made by Wu and Zaffaroni,
(2018), but here we generalize them to allow for a low rank plus sparse
structure. Three comments are necessary. First, notice that we need to control
the entire $\|.\|_{1}$ norm of the coefficients $C_{s}$ of the idiosyncratic
filters. Obviously, we could do the same for the coefficients $B_{s}$ of the
common filters, since if we bound their $l_{1}$ norm their $l_{1,v}$ norm is
bounded as well, because $r$ is finite. Second, part (iv) bounds the column
sums of the coefficients of the idiosyncratic filters. This, together with
part (iii), implies that the diverging behavior of those coefficients is due
to the row sums. This means that for each given $\epsilon_{it}$ the $p$
idiosyncratic innovations have a finite impact. This is just a useful way of
parametrizing the model and we notice that, equivalently, we could also assume
the viceversa or let both row and column sums diverge (compatibly with part
(iii)). It is also straightforward to see that parts (iii) and (iv) are
compatible with Assumption 3.6(i) or 3.7(i). Third, the assumption
$\delta^{\prime}<\alpha$ ensures that the low rank component dominates the
sparse component, in agreement with the idea of common factors.

##  4 Consistency

Suppose now that we observe a sample of $p$-dimensional data vectors with size
$T$, i.e., we observe $\\{X_{it},i=1,\ldots,p,\,t=1,\ldots,T\\}$, and we
compute the estimator
$\widehat{\Sigma}(\theta)=\widehat{L}(\theta)+\widehat{S}(\theta)$ such that:

$\left(\widehat{L}(\theta),\widehat{S}(\theta)\right)=\arg\\!\\!\\!\\!\\!\\!\\!\min_{\begin{subarray}{c}\underline{L}(\theta)\in\mathcal{L}(r)\\\
\underline{S}(\theta))\in\mathcal{S}(q)\end{subarray}}\frac{1}{2}\|\widetilde{\Sigma}(\theta)-(\underline{L}(\theta)+\underline{S}(\theta))\|_{F}^{2}+\psi\|\underline{L}(\theta)\|_{*}+\rho\|\underline{S}(\theta)\|_{1},$
(13)

where $\widetilde{\Sigma}(\theta)$ is the smoothed periodogram defined in (2).

In this section, we prove the algebraic and parametric consistency of the pair
of estimates $(\widehat{L}(\theta),\widehat{S}(\theta))$, and, in order to do
this, we introduce two definitions, taken from Chandrasekaran et al.,  (2012).
First, we say that $(\widehat{S}(\theta),\widehat{L}(\theta))$ is
algebraically consistent if the following conditions hold, for any given
$\theta\in[-\pi,\pi]$:  1.
$\text{{rk}}(\widehat{L}(\theta))=\text{rk}({L}(\theta))$;  2.
$\text{sgn}(\mathrm{Re}(\widehat{S}(\theta)_{ij}))=\text{sgn}(\mathrm{Re}({S}(\theta)_{ij}))$,
for all $i,j=1,\ldots,p$ (by convention we let $\text{sgn}(0)=0$);  3.
$\widehat{L}(\theta)+\widehat{S}(\theta)$and $\widehat{S}(\theta)$ are
positive definite and $\widehat{L}(\theta)$ is positive semidefinite.  Notice
that condition 2 is often referred to also as sparsistency (Chandrasekaran et
al., , 2012).

Second, we say that (parametric) consistency holds if the estimates
$(\widehat{S}(\theta),\widehat{L}(\theta))$ are close to
$({S}(\theta),{L}(\theta))$, for any given $\theta\in[-\pi,\pi]$, with high
probability, according to the following norm:

$g_{\gamma}(\widehat{L}(\theta),\widehat{S}(\theta))=\max\left(\frac{\|\widehat{S}(\theta)-{S}(\theta)\|_{\infty}}{\gamma},\|\widehat{L}(\theta)-{L}(\theta)\|_{2}\right)$
(14)

where $\gamma=\frac{\rho}{\psi}$ is the ratio of the thresholds in (13).

By properly adapting the results of Wu and Zaffaroni,  (2018) to the
intermediate spikiness-sparsity regimes described in Section 3, we prove
uniform consistency over frequencies of the smoothed periodogram pre-estimator
(2).

######  Lemma 4.1.

Suppose that the kernel function $K(.)$ is even, bounded, with support
$[-1,1]$, and bandwidth $M_{T}$, such that:  (i) for some $\kappa>0$,
$|K(s)-1|=O(s^{\kappa})$, as $s\rightarrow 0$;  (ii)
$\int_{-\infty}^{\infty}K^{2}(s)ds<\infty$;  (iii)
$\sum_{s^{\prime}\in\mathbb{Z}}\sup_{|s^{\prime}-s^{\prime\prime}|\leq
1}|K(s^{\prime}\vartheta)-K(s^{\prime\prime}\vartheta)|=O(1)$, as
$\vartheta\rightarrow 0$;  (iv) $c_{1}T^{\underline{\zeta}}\leq M_{T}\leq
c_{2}T^{\zeta}$, for some $c_{1},c_{2}>0$ and $\zeta,\underline{\zeta}>0$,
with $0<\underline{\zeta}<\zeta<1<\underline{\zeta}(2\kappa+1)$.
Then, under Assumption 3.8, there exists some positive real $G$, independent
of $p$ and $T$, such that, as $p,T\to\infty$, for
$\theta_{h}=\frac{h\pi}{M_{T}}$: $\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widetilde{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}\leq
G\sqrt{\frac{M_{T}\log(M_{T})}{T}}\right)\rightarrow 1$.

######  Remark 4.1 (Bandwidth choice).

Notice that the bias term, which is of order $\frac{1}{M_{T}^{\kappa}}$, is
not included in the above result, since for all $M_{T}$ satisfying condition
(iv) this term is always dominated by the variance term. Indeed, while the
optimal choice balancing variance and squared bias is
$M_{T}=O(T^{1/(2\kappa+1)})$, in condition (iv) we are instead assuming
$M_{T}=O(T^{\zeta})$ with $\zeta>\underline{\zeta}>\frac{1}{2\kappa+1}$. In
other words, similarly to Wu and Zaffaroni,  (2018), with this choice of
$M_{T}$ the mean squared error of the smoothed periodogram is dominated by the
variance, while the squared bias becomes negligible, as $T\to\infty$. Typical
values of $\kappa$ are 1 if we choose the Bartlett kernel, or 2 if we choose
the Parzen kernel. All following theoretical results are unaffected if we
relaxed our bandwidth choice and we picked smaller values of $\zeta$, provided
that, when needed, we also account for the bias in the bound in Lemma 4.1.

We are now ready to show parametric and latent rank consistency of
$(\widehat{L}(\theta),\widehat{S}(\theta))$ under the basic filter setting.

######  Theorem 4.1.

Let $\Omega=\Omega({S}^{*})$ and $\mathcal{T}=T({L}^{*})$. Suppose that the
assumptions of Lemma 4.1 hold, with Assumptions 3.1, 3.2, 3.4, and 3.6. Set
$\psi=\frac{p^{\alpha}}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log M_{T}}{T}}$ and
$\rho=\gamma\psi$, where
$\gamma\in[9\xi(\mathcal{T}),\frac{1}{6\mu(\Omega)}]$. In addition, suppose
that
$\underline{\delta}_{T}p^{2(\alpha-\underline{\delta})}<T<\overline{\delta}_{T}p^{6\delta}$
for some $\underline{\delta}_{T},\overline{\delta}_{T}$ such that
$0<\underline{\delta}_{T}<\overline{\delta}_{T}$, and the minimum eigenvalue
of $L^{*}$ is such that
$\lambda_{r}(L^{*})>G_{2}\frac{\psi}{\xi^{2}(\mathcal{T})}$. Then, there
exists a positive real $G$ independent of $p$ and $T$ such that, as
$p,T\to\infty$, for $\theta_{h}=\tfrac{h\pi}{M_{T}}$:

  1. 1.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{L}(\theta_{h})-{L}(\theta_{h})\|_{2}\leq
G\frac{1}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log M_{T}}{T}}\right)\rightarrow
1$;

  2. 2.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{S}(\theta_{h})-{S}(\theta_{h})\|_{\infty}\leq
G\frac{\gamma}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log
M_{T}}{T}}\right)\rightarrow 1$ and
$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{S}(\theta_{h})-{S}(\theta_{h})\|_{2}\leq
Gq^{\prime}\sqrt{\frac{M_{T}\log M_{T}}{T}}\right)\rightarrow 1$;

  3. 3.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}\leq
G\left[\frac{1}{\xi(\mathcal{T})}+q^{\prime}\right]\sqrt{\frac{M_{T}\log
M_{T}}{T}}\right)\rightarrow 1$.

Moreover, if $\psi<C$ for some positive real $C$, then, for
$\theta_{h}=\frac{h\pi}{M_{T}}$:  4.
$\mathcal{P}\left(\text{{rk}}(\widehat{L}(\theta_{h}))=r\right)\rightarrow 1$;
and, if the minimum absolute value of the nonzero off-diagonal entries of
$S^{*}$ is such that $\|S^{*}\|_{min,off}>G_{3}\frac{\psi}{\mu(\Omega)}$ then,
for $\theta_{h}=\frac{h\pi}{M_{T}}$:  5.
$\mathcal{P}\left(\text{sgn}(\mathrm{Re}(\widehat{S}(\theta_{h})_{ij})=\text{{sgn}}(\mathrm{Re}({S}(\theta_{h})_{ij}))\right)\rightarrow
1$, for all $i,j=1,\ldots,p$.

Similarly, under the generalized linear filter setting we have the following.

######  Theorem 4.2.

Let $\Omega=\Omega({S}_{\infty})$ and $\mathcal{T}=T({L}^{*})$. Suppose that
the assumptions of Lemma 4.1 hold, with Assumptions 3.1, 3.3, 3.5, and 3.7.
Set $\psi=\frac{p^{\alpha}}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log M_{T}}{T}}$
and $\rho=\gamma\psi$, where
$\gamma\in[9\xi(\mathcal{T}),\frac{1}{6\mu(\Omega)}]$. In addition, suppose
that
$\underline{\delta}_{T}p^{2(\alpha-\underline{\delta})}<T<\overline{\delta}_{T}p^{6\delta}$
for some $\underline{\delta}_{T},\overline{\delta}_{T}$ such that
$0<\underline{\delta}_{T}<\overline{\delta}_{T}$, and that, for all
$\theta\in[-\pi,\pi]$, the minimum eigenvalue of $L(\theta)$ is such that
$\lambda_{r}(L(\theta))>G_{2}\frac{\psi}{\xi^{2}(\mathcal{T})}$. Then, there
exists a positive real $G$ independent of $p$ and $T$ such that, as
$p,T\to\infty$, for $\theta_{h}=\tfrac{h\pi}{M_{T}}$:

  1. 1.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{L}(\theta_{h})-{L}(\theta_{h})\|_{2}\leq
G\frac{1}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log M_{T}}{T}}\right)\rightarrow
1$;

  2. 2.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{S}(\theta_{h})-{S}(\theta_{h})\|_{\infty}\leq
G\frac{\gamma}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log
M_{T}}{T}}\right)\rightarrow 1$ and
$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{S}(\theta_{h})-{S}(\theta_{h})\|_{2}\leq
Gq^{\prime}\sqrt{\frac{M_{T}\log M_{T}}{T}}\right)\rightarrow 1$;

  3. 3.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}\leq
G\left[\frac{1}{\xi(\mathcal{T})}+q^{\prime}\right]\sqrt{\frac{M_{T}\log
M_{T}}{T}}\right)\rightarrow 1$.

Moreover, if $\psi<C$ for some positive real $C$, then, for
$\theta_{h}=\frac{h\pi}{M_{T}}$:  4.
$\mathcal{P}\left(\text{{rk}}(\widehat{L}(\theta_{h}))=r\right)\rightarrow 1$;
and if, for all $\theta\in[-\pi,\pi]$, the minimum absolute value of the
nonzero off-diagonal entries of $S(\theta)$ is such that
$\|S(\theta)\|_{min,off}>G_{3}\frac{\psi}{\mu(\Omega)}$, then, for
$\theta_{h}=\frac{h\pi}{M_{T}}$:  5.
$\mathcal{P}\left(\text{sgn}(\mathrm{Re}(\widehat{S}(\theta_{h})_{ij})=\text{{sgn}}(\mathrm{Re}({S}(\theta_{h})_{ij}))\right)\rightarrow
1$, for all $i,j=1,\ldots,p$.

Some important remarks follow.

######  Remark 4.2.

The upper bound $T<\overline{\delta}_{T}p^{6\delta}$ is a non-asymptotic
condition necessary to to ensure that the conditions of the Theorem and
Assumptions 3.4-3.5, requiring the eigenvalue $\lambda_{r}({L}(\theta))$ to
diverge as $p^{\alpha}$, are satisfied. Indeed, from the conditions of the
Theorem and by Assumptions 3.4-3.5, we must have:

$\lambda_{r}({L}(\theta))>G_{2}\frac{\psi}{\xi^{2}(\mathcal{T})}=G_{2}\frac{p^{\alpha}}{\xi^{3}(\mathcal{T})}\sqrt{\frac{M_{T}\log
M_{T}}{T}}=Kp^{\alpha},$

for some positive real $K$. Therefore, since by Assumption 3.6-3.7
$\xi(\mathcal{T})=\frac{\sqrt{r}}{\kappa_{L}p^{\delta}}$, and from the
conditions of Lemma 4.1 $M_{T}\geq c_{1}T^{\underline{\zeta}}$ with
$\underline{\zeta}>0$, then it must hold that $p^{3\delta}>MT^{1/2}$, i.e.
$T<\overline{\delta}_{T}p^{6\delta}$ with
$\overline{\delta}_{T}=\frac{1}{K^{2}}$. Let us stress that this is a non-
asymptotic condition, that is, it must hold for all $p$ and $T$. Notice that
equivalently this requires $p>T^{1/(6\delta)}$ which shows that the less
sparse is $S(\theta)$ (higher $\delta$) the larger must be $p$ in order to
ensure the eigen-gap to be large enough for identification of the latent rank
and the sparsity pattern (parts 4 and 5 of the theorems) to be possible even
for finite $p$ and $T$. Notice also that if $\delta=0$ we can still have
parametric consistency (parts 1, 2, and 3 of the theorems) as long as
$p,T\to\infty$ but $\frac{{T}}{p}\to 0$.

######  Remark 4.3.

The lower bound $T>\underline{\delta}_{T}p^{2(\alpha-\underline{\delta})}$ is
also a non-asymptotic condition necessary to to ensure that the conditions of
the Theorem and Assumptions 3.6-3.7, requiring that the maximum number of non-
zeros per row in $S(\theta)$ to diverge as $p^{\delta}$, are satisfied.
Indeed, from the conditions of the Theorem and Assumptions 3.6-3.7

$\delta_{2}^{\prime}p^{\underline{\delta}}\geq\|S(\theta)\|_{1,v}\geq
q^{\prime}\|S(\theta)\|_{min,off}>G_{3}\frac{q^{\prime}\psi}{\mu(\Omega)}\geq
G_{3}\frac{C^{\prime}p^{\alpha}}{\xi(\mathcal{T})\mu(\Omega)}\sqrt{\frac{M_{T}\log
M_{T}}{T}},$

since $q^{\prime}\geq C^{\prime}$ for some positive real $C^{\prime}$.
Therefore, since $\xi(\mathcal{T})\mu(\Omega)\leq\frac{1}{54}$ by Assumptions
3.6-3.7 and from the conditions of Lemma 4.1 $M_{T}\geq
c_{1}T^{\underline{\zeta}}$ with $0<\underline{\zeta}<1$, we must have
$p^{\underline{\delta}}>Cp^{\alpha}T^{-1/2}$ for some positive real $C$ or
equivalently $T>\underline{\delta}_{T}p^{2(\alpha-\underline{\delta})}$ with
$\underline{\delta}_{T}=C^{2}$. Notice that
$\underline{\delta}<\delta^{\prime}<\alpha$ by Assumption 3.8 and since
$\delta^{\prime}$ is a bound on the $l_{1}$ norm while $\underline{\delta}$ is
a bound on the $l_{1,v}$ norm. Moreover,
$\underline{\delta}\leq\delta+\frac{1}{2}$ so when $S(\theta)$ is the least
sparse possible ($\delta=\frac{1}{2}$) we do not need a very large $T$, the
lower bound being $T>\underline{\delta}_{T}p^{2(\alpha-1)}$, and if the latent
eigenvalues are very spiked ($\alpha=1$) the sparsity pattern can be
identified without imposing constraints between $T$ and $p$. In the most
sparse case ($\delta=0$) we need a
$T>\underline{\delta}_{T}p^{2(\alpha-1/2)}$, which in the spiked case
($\alpha=1$) implies $T$ at least comparable to $p$, while, for the upper
bound in Remark 4.2 to still hold, we must have both $p$ and $T$ bounded by a
constant (recall that these are non-asymptotic conditions). If this lower
bound for $T$ is not satisfied we cannot identify the sparsity pattern, but we
can still recover the latent rank, and the parametric consistency of both
$\widehat{L}(\theta_{h})$ and $\widehat{S}(\theta_{h})$ still holds.

######  Remark 4.4.

Parts 1, 2, and 3 provide bounds for the estimation error of the spectral
density matrices, which are uniform over all frequencies. Parts 4 and 5
guarantee rank consistency for $\widehat{L}(\theta_{h})$ and sparsistency for
$\widehat{S}(\theta_{h})$. In order for these conditions to be verified we
need $\psi=\frac{p^{\alpha}}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log
M_{T}}{T}}$ to be finite for all $p$ and $T$, however it is not required for
$\psi$ to decrease as $p$ and $T$ increase. Furthermore, when $\psi$ is finite
for all $p$ and $T$, then there exists a constant $\varphi$, depending on
$\mu(\Omega)$ and $\xi(\mathcal{T})$, such that if
$g_{\gamma}(\widehat{L}(\theta),\widehat{S}(\theta))\leq\varphi$ then parts 4
and 5 hold with probability 1 for all $p$ and $T$ (see Chandrasekaran et al.,
, 2012, Propositions 5.2 and 5.3, for details).

######  Remark 4.5.

The error bound in spectral norm for the sparse component in part 2 is larger
than the bound $O\left(\sqrt{\frac{\log{p}}{T}}\right)$ reported in Bickel and
Levina,  (2008), as we are allowing $q^{\prime}$ to grow as $p^{\delta}$.

######  Remark 4.6.

From parts 1 and 2, it immediately follows that
$\mathcal{P}\left(g_{\gamma}(\widehat{L}(\theta),\widehat{S}(\theta))\leq
G\frac{p^{\alpha}}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log M_{T}}{T}}\right)\to
1$. Therefore, since $\xi(\mathcal{T})=O(p^{\delta})$ by Assumption 3.7, then,
if $p^{\alpha+\delta}\sqrt{\frac{M_{T}\log M_{T}}{T}}\to 0$, as
$p,T\to\infty$, then $g_{\gamma}(\widehat{L}(\theta),\widehat{S}(\theta))\to
0$, with probability tending to 1. Notice that, to achieve parametric
consistency, we must have $T$ growing faster than $p$, and in the worst case,
i.e., $\alpha=1$ and $\delta=\frac{1}{2}$, this means that we need
$p^{3}{\frac{M_{T}\log M_{T}}{T}}\to 0$. However, notice also that if, as
common in the literature on the estimation of large matrices, we were to
consider the norm relative to the dimension $p$, then the error bound would be
$p^{\alpha+\delta-1}\sqrt{\frac{M_{T}\log M_{T}}{T}}$, which in the worst case
requires $p{\frac{M_{T}\log M_{T}}{T}}\to 0$. On the other hand in the GDFM
case, i.e., when $\alpha=1$ and $\delta=0$, we would have the same bound,
$\sqrt{\frac{M_{T}\log M_{T}}{T}}$, which was stated in Theorem 2.1. This
bound is also the same one derived for the classical smoothed periodogram
estimator in Lemma 4.1.

Finally, we have a useful Corollary about the inverses of
$\widehat{S}(\theta_{h})$ and $\widehat{\Sigma}(\theta_{h})$.

######  Corollary 4.1.

Under the assumptions of Theorem 4.1 or Theorem 4.2, there exists a positive
real $G$ independent of $p$ and $T$ such that, as $p,T\to\infty$, for
$\theta_{h}=\frac{h\pi}{M_{T}}$:

  1. 1.

$\widehat{\Sigma}(\theta_{h})$ is positive definite if
$\frac{\lambda_{p}(\Sigma(\theta_{h}))}{p^{\alpha}}>G\left(\frac{1}{\xi(\mathcal{T})\sqrt{T}}+\frac{q^{\prime}}{\sqrt{T}}\right)\sqrt{\frac{M_{T}\log
M_{T}}{T}}$;

  2. 2.

$\widehat{S}(\theta_{h})$ is positive definite if
$\frac{\lambda_{p}(S(\theta_{h}))}{p^{\alpha}}>G\frac{q^{\prime}}{\sqrt{T}}\sqrt{\frac{M_{T}\log
M_{T}}{T}}$;

  3. 3.

$\widehat{\Sigma}^{-1}(\theta_{h})$ is positive definite if
$\frac{\lambda_{p}(\Sigma(\theta_{h}))}{p^{\alpha}}\geq
2G\left(\frac{1}{\xi(\mathcal{T})\sqrt{T}}+\frac{q^{\prime}}{\sqrt{T}}\right)\sqrt{\frac{M_{T}\log
M_{T}}{T}}$;

  4. 4.

$\widehat{S}^{-1}(\theta_{h})$ is positive definite if
$\frac{\lambda_{p}(S(\theta_{h}))}{p^{\alpha}}\geq
2G\frac{q^{\prime}}{\sqrt{T}}\sqrt{\frac{M_{T}\log M_{T}}{T}}$.

In addition, it also holds for $\theta_{h}=\tfrac{h\pi}{M_{T}}$:

  1. 5.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{\Sigma}(\theta_{h})^{-1}-\Sigma(\theta_{h})^{-1}\|_{2}\leq
G\left(\frac{1}{\xi(\mathcal{T})\sqrt{T}}+\frac{q^{\prime}}{\sqrt{T}}\right)\sqrt{\frac{M_{T}\log
M_{T}}{T}}\right)\rightarrow 1$;

  2. 6.

$\mathcal{P}\left(\max_{|h|\leq
M_{T}}\frac{1}{p^{\alpha}}\|\widehat{S}(\theta_{h})^{-1}-S(\theta_{h})^{-1}\|_{2}\leq
G\frac{q^{\prime}}{\sqrt{T}}\sqrt{\frac{M_{T}\log M_{T}}{T}}\right)\rightarrow
1$.

######  Remark 4.7 (Unshrinking).

We stress that above defined estimates may suffer from systematic sub-
optimality for what concerns estimated eigenvalues. In particular, if $p$ is
large and the latent eigenvalues are spiked, the singular value thresholding
procedure may lead to the over-shrinkage of latent eigenvalues. For this
reason, following Farnè and Montanari,  (2020), we perform the un-shrinkage of
the estimated latent eigenvalues, i.e., we give back the threshold to
$\widehat{L}(\theta_{h})$. The new idiosyncratic estimate is then obtained
keeping fixed the off-diagonal sparsity pattern recovered, and deriving its
diagonal by difference from the diagonal of $\widehat{\Sigma}(\theta_{h})$.
The resulting matrix estimators are our UNshrunk ALgebraic Spectral Estimators
(UNALSE). More specifically, for any $\theta_{h}=\frac{h\pi}{M_{T}}$, consider
the spectral decomposition
$\widehat{L}(\theta_{h})=\widehat{W}(\theta_{h})\widehat{D}(\theta_{h})\widehat{W}^{\dagger}(\theta_{h})$,
then we define:

$\displaystyle\widehat{L}_{\text{\tiny{UNALSE}}}(\theta_{h})=\widehat{W}(\theta_{h})\left(\widehat{D}(\theta_{h})+{\psi}{I}_{r}\right)\widehat{W}^{\dagger}(\theta_{h}),$
$\displaystyle\text{diag}(\widehat{S}_{\text{\tiny{UNALSE}}}(\theta_{h}))=\text{diag}(\widehat{\Sigma}(\theta_{h}))-\text{diag}(\widehat{L}(\theta_{h})),\quad\text{off-
diag}(\widehat{S}_{\text{\tiny{UNALSE}}}(\theta_{h}))=\text{off-
diag}(\widehat{S}(\theta_{h})),$
$\displaystyle\widehat{\Sigma}_{\text{\tiny{UNALSE}}}(\theta_{h})=\widehat{L}_{\text{\tiny{UNALSE}}}(\theta_{h})+\widehat{S}_{\text{\tiny{UNALSE}}}(\theta_{h}),$

where ${\psi}>0$ is the same as in Theorems 4.1 or 4.2. The above defined
UNALSE estimates have two relevant optimality properties. First, they have the
smallest possible Frobenius loss from the targets into the recovered matrix
varieties. Second, they have the maximally concentrated eigenvalues into the
class of algebraically consistent estimators, for any given sample size $T$.
Under our assumptions, the above optimality properties of the un-shrinkage
procedure of Farnè and Montanari,  (2020), to which we refer for the details,
hold straightforwardly.

##  5 Threshold selection

In solving problem (3), the choice of the eigenvalue threshold $\psi$ and the
sparsity threshold $\rho$ is a nontrivial issue. Differently from the
covariance matrix context, in fact, the magnitude of the eigenvalues can vary
a lot across frequencies, which may cause the optimization of (3) to be
strongly sensitive to the magnitude of both thresholds.

Let us suppose that $\widehat{L}_{\psi,\rho}(\theta_{h})$,
$\widehat{S}_{\psi,\rho}(\theta_{h})$,
$\widehat{\Sigma}_{\psi,\rho}(\theta_{h})=\widehat{L}_{\psi,\rho}(\theta_{h})+\widehat{S}_{\psi,\rho}(\theta_{h})$
are the solutions of (3) with thresholds $\psi$ and $\rho$, under the
assumptions of Theorem 4.1. The dual norm of the composite loss (14) is
considered, from which we define at each $\theta_{h}$ the following criterion:

$MC_{h}(\psi,\rho)=\max\left\\{\frac{\widehat{r}\|\widehat{L}_{\psi,\rho}(\theta_{h})\|_{2}}{\widehat{\beta}_{\psi,\rho}(\theta_{h})},\frac{\frac{\psi}{\rho}\|\widehat{S}_{\psi,\rho}(\theta_{h})\|_{1,v}}{(1-\widehat{\beta}_{\psi,\rho}(\theta_{h}))}\right\\},$
(15)

where
$\|\widehat{S}_{\psi,\rho}(\theta_{h})\|_{1,v}=\max_{i=1,\ldots,p}\sum_{j=1}^{p}|\widehat{S}_{ij,\psi,\rho}(\theta_{h})|$
and
$\widehat{\beta}_{\psi,\rho}(\theta_{h})=\text{tr}(\widehat{L}_{\psi,\rho}(\theta_{h}))/\text{tr}(\widehat{\Sigma}_{\psi,\rho}(\theta_{h}))$
is the estimated proportion of latent variance. The optimal threshold pair
$(\breve{\psi}_{h},\breve{\rho}_{h})$ is thus selected as the mini-max
$(\breve{\psi}_{h},\breve{\rho}_{h})=\arg\min_{\psi,\rho}{MC_{h}(\psi,\rho)}$,
where $\psi$ and $\rho$ vary across pre-specified grids. This threshold
selection method penalizes solution pairs with too dispersed latent
eigenvalues and too many residual nonzeros in single rows, by comparing two
appropriately re-scaled versions of the spectral norm of the low rank solution
and the row-wise maximum norm of the residual solution.

In order to ensure the effectiveness of the above criterion, the threshold
grids need to be properly set up at each $\theta_{h}$, according to the
unknown underlying algebraic structure. We thus recall from Theorem 4.2 that
$\psi=\sqrt{\frac{pM_{T}\log M_{T}}{T}}\frac{1}{\xi(T)}$, assuming the
intermediate value $\alpha=1/2$ and recalling from Chandrasekaran et al.,
(2011) that $\text{inc}(L(\theta_{h}))\leq\xi(T)\leq
2\text{inc}(L(\theta_{h}))$, where $\text{inc}(L(\theta_{h}))$ is the
incoherence of $L(\theta_{h})$, defined as
$\text{inc}(L(\theta_{h}))=\max_{i=1,\ldots,p}\|\mathcal{P}e_{i}\|$, with
$e_{i}$ the canonical basis vector ($i$th column of the $p$ dimensional
identity matrix), and the operator $\mathcal{P}$ projecting each $e_{i}$ onto
the row/column space of $L(\theta_{h})$. The extreme incoherence values are
$\text{inc}(L(\theta_{h}))=1$, when any vector of the standard basis belongs
to the row/column space of $L(\theta_{h})$, and
$\text{inc}(L(\theta_{h}))=\sqrt{\frac{r}{p}}$, when $L(\theta_{h})$ is a
Hadamard matrix. In light of this, and since
$\sqrt{\frac{M\log{M_{T}}}{T}}>\sqrt{\frac{1}{T}}$, we initialize the grid for
the eigenvalue threshold ${\breve{\psi}}$ as the sequence of $n_{thr}$ equi-
spaced real numbers from $\sqrt{\frac{p}{T}}\frac{1}{2\widetilde{inc}}$ to
$\sqrt{\frac{p}{T}}\frac{1}{\widetilde{inc}}$, where, for any given value of
$r_{thr}$, we set $\widetilde{inc}=\sqrt[4]{\frac{r_{thr}}{p}}$, which is the
geometric mean of the minimum and maximum incoherence values, i.e.
$\sqrt{\frac{r_{thr}}{p}}$ and $1$ respectively (see Chandrasekaran et al.,
(2011)).

We start by setting $r_{thr}=1$. Then, we apply the solution algorithm of (3)
(see Section 2) with the grid for $\breve{\psi}$ as defined above. Note that
we iteratively adapt the grid for $\breve{\psi}$ during the optimization
process by dividing, at each step $k$ of the iteration, the grid components by
the Gini index of the eigenvalues of
$\mathcal{E}_{Y,k}(\theta_{h})=Y_{k-1}(\theta_{h})-\frac{1}{2}(Y_{k-1}(\theta_{h})+Z_{k-1}(\theta_{h})-\widetilde{\Sigma}(\theta_{h}))$
(see point 2b of the solution algorithm). This is done in order to adapt the
eigenvalue thresholds to the underlying degree of spikiness of latent
eigenvalues.

Then, if the ranks of $\widehat{L}_{\breve{\psi},\breve{\rho}}(\theta_{h})$
across thresholds vary too much or the eigenvalue threshold selected by
criterion (15) lies in the grid extremes, the value of $r_{thr}$ must be
changed. In particular, it must be decreased if the recovered rank is very
large uniformly across thresholds, and increased if it is very small or zero.
When the eigenvalue threshold selected by the MC criterion is far away from
the boundaries, and the recovered rank is constant and stable across
thresholds, we stop and select $r_{thr}$.

Concerning the sparsity threshold, we decide to set $\breve{\gamma}$ as the
sequence of $n_{thr}$ equi-spaced real numbers from $s_{thr}\times p^{-1/2}$
to $s_{thr}\times p^{-1/4}$. The two functions of $p$, $p^{-1/2}$ and
$p^{-1/4}$, represent two plausible extremes for residual nonzero proportions,
while $s_{thr}$ is a magnitude parameter. We first set $s_{thr}$ to $1$, and
we run the solution algorithm of (3) setting
$\breve{\rho}=\breve{\gamma}\sqrt{\frac{p}{T}}\frac{1}{\widetilde{inc}}$.

Similarly, if the sparsity threshold selected by criterion (15) lies in the
grid extremes, we decrease or increase it, in order to obtain a non-diagonal
solution with a reasonable and stable number of nonzeros. When the value of
$\breve{\gamma}$ selected by the MC criterion is far from the grid extremes
and the number of nonzeros is approximately stable across thresholds, we stop
and select $s_{thr}$. The described process allows to map the problem of
selecting thresholds $\psi$ and $\rho$ into the more intuitive problem of
selecting $r_{thr}$ and $s_{thr}$. Apart from it, the two thresholds can also
be chosen manually, as long as the recovered rank and sparsity pattern are
constant in a neighborhood of the chosen threshold pairs.

##  6 Simulation study

In order to test the performance of UNALSE under an exhaustive range of
situations, we consider three different simulated scenarios:  A : simulated
with basic filters;  B : simulated with general filters and a very sparse,
almost negligible, residual pattern;  C : simulated with general filters and a
less sparse, more relevant, residual pattern.  For each Scenario, we setup
five Settings, with different dimensions, sample sizes, spectral magnitudes,
latent ranks, condition numbers, and sparsity degrees. In particular, the
Settings from $1$ to $3$ present $p=100$ and $T=1000$. Setting 4 presents
$p=T=150$, Setting 5 presents $p=200$ and $T=100$.

Our spectra follow a reverse S-shape: the minimum latent eigenvalue decreases
across frequencies for all settings, the minimum off-diagonal nonzero residual
entry (in absolute value) varies with a similar trend across frequencies. The
proportion of latent variance differs across settings. In absolute terms, the
largest one is for Setting 4, followed by Setting 5, and Settings 1, 2 and 3.
At the same time, considering the proportion of residual covariance, Settings
1 and 4 are the most sparse, followed by Settings 2 and 5, while Setting 3 is
the least sparse. In Scenario C we have a situation with a small minimum
latent eigenvalue and minimum residual non-zero off-diagonal entry (in
absolute terms) across frequencies, and the capability of UNALSE to recover
the sparsity pattern can be meaningfully tested.

We fix the frequency grid as $\theta_{h}=\frac{\pi h}{12}$, $h=0,\ldots,5$. We
simulate $N=100$ replications of the data $X_{t}$ having a spectral density
with a low rank plus sparse structure designed according to a given Scenario
and Setting (see Section D in the supplementary appendix for details on the
simulation mechanism), computing the pre-estimator of the spectral density
matrix (2) on each simulated dataset. To this end, we adopt the classical
choice of a Bartlett kernel (Forni et al., , 2000, 2017) and we set
$M_{T}=\lfloor\sqrt{T}\rfloor$.

For each replication, we apply the solution algorithm (see Section 2) with the
threshold selection procedure described in Section 5, thus getting $100$
optimal spectral density matrix estimates
$\widehat{L}^{(b)}_{\text{\tiny{UNALSE}}}(\theta_{h})$,
$\widehat{S}^{(b)}_{\text{\tiny{UNALSE}}}(\theta_{h})$,
$\widehat{\Sigma}^{(b)}_{\text{\tiny{UNALSE}}}(\theta_{h})$, for
$b=1,\ldots,100$, $h=0,\ldots,5$. In the literature, the only existing
competitor is for the low rank component and relies on the dynamic principal
components of Brillinger,  (2001) (see also Forni et al.,  (2000)). We call it
$\widehat{L}_{DYN}(\theta_{h})$. Note however that this estimator requires a
pre-specified rank, which, in a high-dimensional setting, can be determined
via the information criterion by Hallin and Liška,  (2007) or the test by
Onatski,  (2009).

Results are reported using the frequencies ${f}_{h}=\frac{\theta_{h}}{\pi}$,
$h=0,\ldots,5$. For each $f_{h}$ and for each replication, we calculate some
relevant statistics and some relative metrics to evaluate the quality of the
rank and sparsity pattern recovery.

  1. (i)

the latent variance proportion
$\widehat{\beta}(f_{h})=\frac{1}{100}\sum_{b=1}^{100}\frac{tr(\widehat{L}^{(b)}(f_{h}))}{tr(\widehat{\Sigma}^{(b)}(f_{h}))}$;

  2. (ii)

the binary indicator sum for the correct estimated rank
$\widehat{R}=\frac{1}{6}\sum_{b=1}^{100}\sum_{h=0}^{5}\mathbb{I}(\hat{r}^{(b)}(f_{h})=r)$.

  3. (iii)

the nonzero predictive value:
$nzpv(f_{h})=\frac{1}{100}\sum_{b=1}^{100}\frac{\sum_{i=1}^{p}\sum_{j=i+1}^{p}\mathbb{I}\\{\widehat{S}^{(b)}_{ij}(f_{h})\neq
0\bigcup{S}_{ij}(f_{h})\neq
0\\}}{\sum_{i=1}^{p}\sum_{j=i+1}^{p}{\mathbb{I}\\{\widehat{S}_{ij}^{(b)}(f_{h})\neq
0\\}}}$;

  4. (iv)

the positive predictive value:
$ppv(f_{h})=\frac{1}{100}\sum_{b=1}^{100}\frac{\sum_{i=1}^{p}\sum_{j=i+1}^{p}\mathbb{I}\\{\widehat{S}^{(b)}_{ij}(f_{h})>0\bigcup{S}_{ij}(f_{h})>0\\}}{\sum_{i=1}^{p}\sum_{j=i+1}^{p}{\mathbb{I}\\{{S}_{ij}^{(b)}(f_{h})>0\\}}}$;

  5. (v)

the negative predictive value:
$npv(f_{h})=\frac{1}{100}\sum_{b=1}^{100}\frac{\sum_{i=1}^{p}\sum_{j=i+1}^{p}\mathbb{I}\\{\widehat{S}^{(b)}_{ij}(f_{h})<0\bigcup{S}_{ij}(f_{h})<0\\}}{\sum_{i=1}^{p}\sum_{j=i+1}^{p}{\mathbb{I}\\{{S}_{ij}^{(b)}(f_{h})<0\\}}}$;

  6. (vi)

the maximum of the sum of the binary indicator of nonzero recovered residual
entries of each row: $mnz_{i}(f_{h})=\max_{j=1,\ldots,p,j\neq
i}\sum_{b=1}^{100}\mathbb{I}\\{\widehat{S}^{(b)}_{ij}(f_{h})\neq 0\\}$.

In order to evaluate the properties of our estimates, we calculate also the
Frobenius loss of each estimate from the relative target, rescaled by the
dimension:

  1. (i)

the low rank component Frobenius loss
$err_{\widehat{L}}(f_{h})=\frac{1}{100}\sum_{b=1}^{100}\|\widehat{L}^{(b)}(f_{h})-{L}(f_{h})\|_{F}/p$;

  2. (ii)

the ratio between the overall UNALSE and the input Frobenius loss with respect
to the target:
$err_{ratio}(f_{h})=\frac{1}{100}\sum_{b=1}^{100}\frac{\|\widehat{\Sigma}^{(b)}(f_{h})-{\Sigma}(f_{h})\|_{F}}{\|\widetilde{\Sigma}^{(b)}(f_{h})-{\Sigma}(f_{h})\|_{F}}$.

For all above quantities, we calculate also the standard deviation across the
$100$ trials. We present here results only for Scenarios A and C and Settings
3 and 4, while all other results are available in Section D of the
supplementary appendix.

First of all, we compare our estimates of the latent rank $r$ with those
obtained with the test by Onatski,  (2009), based on sample dynamic
eigenvalues of the smoothed periodogram estimator (see Table 1). We observe
that the test presents some empirical level issues, particularly when the
eigenvalues are not so spiked. For Setting 1, for instance, the observed
proportion of correct decisions is as low as $38\%$. For the other settings of
all scenarios, we note that the same proportion increases considerably,
consistently with the increased spikiness of latent eigenvalues. At the same
time, the observed outcome is always less than our $100\%$.

Table 1: Estimation of latent rank $r$. | Scenario
---|---
| A | B | C
Setting | Onatski | UNALSE | Onatski | UNALSE | Onatski | UNALSE
1 | 38 | 100 | 90 | 100 | 82 | 100
2 | 95 | 100 | 93 | 100 | 86 | 100
3 | 93 | 100 | 94 | 100 | 90 | 100
4 | 90 | 100 | 94 | 100 | 92 | 100
5 | 96 | 100 | 98 | 100 | 95 | 100
Figure 1: Estimated latent variance proportion $\widehat{\beta}(f_{h})$ \- Scenario A. Setting 3 | Setting 4
---|---
|

In Figure 1 we show the estimated latent variance proportion
$\widehat{\beta}(f_{h})$. We notice that UNALSE estimates systematically
better than DYN the true $\beta(f_{h})$ across frequencies (results not
reported), in particularly for Setting 3, which has $p/T=0.1$. This is due to
the bad properties of sample eigenvalues with respect to the presence of non-
spiked latent eigenvalue structures. In addition, UNALSE can recover the
residual sign pattern at each frequency (see Figure 2), and the positive and
negative predictive values are pretty similar. When the overall magnitude is
larger, the true predictive rate of nonzeros across frequencies grows
considerably: indeed, in Setting 3, it overcomes $80\%$ at all frequencies.
This happens because there are many nonzeros of sufficient magnitude, even
compared to the low rank component.

Figure 2: Positive and negative predictive values $ppv(f_{h})$ and $npv(f_{h})$ \- Scenario A. $ppv(f_{h})$ Setting 3 |  $ppv(f_{h})$ Setting 4 |  $npv(f_{h})$ Setting 3 |  $npv(f_{h})$ Setting 4
---|---|---|---
|  |  |

Concerning the error metrics (see Figure 4), $err_{\widehat{L}}(f_{h})$ is
slightly worse for UNALSE than for DYN at low frequencies, and quite better at
high frequencies. Again, this is consistent with the bad properties of the
dynamic principal components under weak factors in not so large dimensions.
The same pattern is visible for $err_{ratio}(f_{h})$.

Figure 3: $err_{\widehat{L}}(f_{h})$ and $err_{ratio}(f_{h})$ \- Scenario A. $err_{\widehat{L}}(f_{h})$ Setting 3 |  $err_{\widehat{L}}(f_{h})$ Setting 4 |  $err_{ratio}(f_{h})$ Setting 3 | $err_{ratio}(f_{h})$ Setting 4
---|---|---|---
|  |  |

Scenario C also shows that the proportion of latent variance is estimated very
well by UNALSE, which is very close, for Setting 3, even to a target lower
than $40\%$ (see Figure 10). Figure 11 shows that UNALSE is still very good
regarding $err_{\widehat{L}}(f_{h})$ and $err_{ratio}(f_{h})$, particularly
for Setting 3 at high frequencies.

Figure 4: Estimated latent variance proportion $\widehat{\beta}(f_{h})$ \- Scenario C. Setting 3 | Setting 4
---|---
|
Figure 5: $err_{\widehat{L}}(f_{h})$ and $err_{ratio}(f_{h})$\- Scenario C. $err_{\widehat{L}}(f_{h})$ Setting 3 |  $err_{\widehat{L}}(f_{h})$ Setting 4 |  $err_{ratio}(f_{h})$ Setting 3 |  $err_{ratio}(f_{h})$ Setting 4
---|---|---|---
|  |  |

Concerning the sparsity pattern, the nonzeros recovered by UNALSE are stable
across frequencies and the nonzero predictive value overcomes $70\%$ at
central frequencies for Setting 3 and it is still acceptable for Setting 4
(see Figure 6). Finally, Figure 7 shows how the indicator $mnz_{i}(f_{h})$,
$h=0,1,2$, presents a very similar pattern across variables for the first
three frequencies, proving that the sparsity pattern is consistent over
frequencies. Note that this property is common to all Scenarios and Settings,
even when the predictive value is not good.

Figure 6: Nonzero predictive values $nzpv(f_{h})$ \- Scenario C. Setting 3 | Setting 4
---|---
|
Figure 7: $mnz_{i}(f_{h})$, $h=0,1,2$ \- Scenario C - Setting 3. $f_{1}=0$ | $f_{2}=0.0833$ | $f_{3}=0.1667$
---|---|---
|  |

##  7 Real data analysis

We consider a dataset of $p=101$ quarterly macroeconomic indicators regarding
the US economy observed over $T=210$ time points spanning the period
1960:Q2-2012:Q3 (see, e.g., McCracken and Ng, , 2020). Throughout, we compute
the smoothed periodogram by setting $M_{T}=\lfloor\sqrt{T}\rfloor=14$ and
using the Bartlett kernel. The analysis that follows shows that the underlying
spectral structure seems to be quite relevant at certain frequencies and
heterogeneous.

In the top left panel of Figure 8 we show the four largest eigenvalues,
rescaled by $p$, of the smoothed periodogram estimator. The top eigenvalue
shows a decreasing shape from $f=0$ and two auxiliary peaks at $f=0.07$ and
$f=0.3$, corresponding to periods of about 3.5 years and 9 months,
respectively. Note that 3.5 years is around the typical period of a business
cycle. The estimated rank by UNALSE is $\widehat{r}=2$ at all frequencies. The
top right panel of Figure 8 shows the proportion of latent variance
$\widehat{\beta}(f_{h})=\frac{\text{rk}(\widehat{L}(f_{h}))}{\text{rk}(\widehat{\Sigma}(f_{h}))}$
so the contribution of $\widehat{L}(f_{h})$, which follows the pattern of the
leading eigenvalues of $\widetilde{\Sigma}(f_{h})$, hence it captures the
business cycle frequency. The bottom left panel of Figure 8 reports the
proportion of residual covariance
$\widehat{\zeta}(f_{h})=\frac{\sum_{i=1}^{p}\sum_{j=i+1}^{p}|\widehat{S}_{ij}(f_{h})|}{\sum_{i=1}^{p}\sum_{j=i+1}^{p}|\widehat{\Sigma}_{ij}(f_{h})|}$,
summarizing the contribution of $\widehat{S}(f_{h})$, which has the main
contribution at a higher frequency $f=0.45$, corresponding to a period of 6
months. Secondary maxima are at $f=0.03$, i.e., a period of 7 years, and
$f=0.25$, corresponding to a period of 1 year. Finally, the fraction of
nonzeros has a similar pattern (see bottom right panel of Figure 8).

Figure 8: US macroeconomic data - Co-movements and sparsity. |  |  |
---|---|---|---
$\frac{\widehat{\lambda}_{j}(\widetilde{\Sigma}(f_{h}))}{p}$, $j=1,2,3,4$. |  $\widehat{\beta}(f_{h})$. |  $\widehat{\zeta}(f_{h})$. | fraction of nonzeros in $\widehat{S}(f_{h})$.

In Figure 9 we show heat-maps of $\widehat{L}(f_{h})$ at frequencies $0$ and
$0.07$.

Figure 9: US macroeconomic data - $\widehat{L}(f_{h})$. |  |
---|---|---
$\widehat{L}(f_{h}=0)$. |  $\text{Re}(\widehat{L}(f_{h}=0.07))$. |  $\text{Im}(\widehat{L}(f_{h}=0.07))$.

The variables having the largest (in term of explained variance) common
component are all related to the labor market:  (i) Civilian Employment;  (ii)
All Employees in Service-Providing Industries, Total Private Industries, and
in Trade, Transportation & Utilities.

In Figure 10 we show heat-maps of $\widehat{S}(f_{h})$ at frequencies $0.03$
and $0.45$.

Figure 10: US macroeconomic data - $\widehat{S}(f_{h})$. |  |  |
---|---|---|---
$\text{Re}(\widehat{S}(f_{h}=0.03))$. |  $\text{Im}(\widehat{S}(f_{h}=0.03))$. |  $\text{Re}(\widehat{S}(f_{h}=0.45))$. |  $\text{Im}(\widehat{S}(f_{h}=0.45))$.

It is worth mentioning some prominent co-spectral relationship at selected
frequencies. At frequency $f=0.03$, the following pairs display strong co-
dependence:  (i) Consumer Loans at All Commercial Banks and Total Consumer
Credit (Owned and Securitized);  (ii) 3-Year Treasury Constant Maturity Rate
and 10-Year Treasury Constant Maturity Rate;  (iii) Compensation Per Hour in
the Manufacturing Sector and in the Business Sector.

At frequency $f=0.45$, we observe a strong relationship for the following
variable pairs:  (i) Consumer Price Index for All Urban Consumers of All Items
Less Energy and of All Items Less Food & Energy;  (ii) Real Imports and
Exports of Goods & Services Per-Capita;  (iii) Real Exports of Goods &
Services and Gross Private Domestic Investment Price Index;  (iv) Unit Labor
Cost in the Business Sector and Output Per Hour of All Persons in the Business
Sector.

##  8 Conclusions

In this paper, we consistently estimate the spectral density matrix under the
assumption of a dynamic factor model with a sparse idiosyncratic multivariate
spectrum for the data, where the relative pervasiveness of the common and the
idiosyncratic components are allowed to vary across frequencies. We prove that
the nuclear norm plus $l_{1}$ norm heuristics consistently recovers across
frequencies the spectral components and their sum, as well as the dynamic rank
and the residual sparsity pattern. We call the resulting estimators UNALSE
(UNshrunk ALgebraic Spectral Estimator).

The empirical implications of the UNALSE approach are discussed on a US
macroeconomic dataset, showing that UNALSE is able to catch the driving
variables of the latent dynamics as well as the particular strength of
specific relationships at all frequencies. This opens up the way to enhanced
dynamic factor scores estimation and temporal network analysis.

##  Appendix A Proofs

#### Proof of Lemma 4.1

First, due to the uncorrelation between $u_{t}$ and $e_{t}$, we can decompose
the expected value of the $(ij)-$ entry of the raw periodogram
$\widehat{\Sigma}^{X}_{Raw,ij}(\theta_{h})$ at each frequency
$\theta_{h}=\pi\frac{h}{M_{T}}$, $h\leq|M_{T}|$, as follows:

$\mathrm{E}(\widehat{\Sigma}^{X}_{Raw,ij}(\theta_{h}))=\mathrm{E}(d_{i}^{X}(\theta_{h})\overline{d_{j}^{X}(\theta_{h})})=\mathrm{E}(d_{i}^{\chi}(\theta_{h})\overline{d_{j}^{\chi}(\theta_{h})})+\mathrm{E}(d_{i}^{\epsilon}(\theta_{h})\overline{d_{j}^{\epsilon}(\theta_{h})}),$

where
$d_{i}^{X}(\theta_{h})=\frac{1}{\sqrt{T}}\sum_{t=1}^{T}X_{i,t}\mathrm{e}^{-\mathrm{i}\theta_{h}t}$,
$d_{i}^{\chi}(\theta_{h})=\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\chi_{i,t}\mathrm{e}^{-\mathrm{i}\theta_{h}t}$,
$d_{i}^{\epsilon}(\theta_{h})=\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\epsilon_{i,t}\mathrm{e}^{-\mathrm{i}\theta_{h}t}$
are the Fourier transforms at frequency $\theta_{h}$ of $X_{t}$, $\chi_{t}$
and $\epsilon_{t}$ respectively.

At each frequency $\theta_{h}=\pi\frac{h}{M_{T}}$, $h\leq|M_{T}|$, we define
the common component of $\widehat{\Sigma}^{X}_{ij,Raw}(\theta_{h})$ as
$\widehat{\Sigma}^{\chi}_{ij,Raw}(\theta_{h})=d_{i}^{\chi}(\theta_{h})\overline{d_{j}^{\chi}(\theta_{h})}$,
and the idiosyncratic component of $\widehat{\Sigma}^{X}_{ij,Raw}(\theta_{h})$
as
$\widehat{\Sigma}^{\epsilon}_{ij,Raw}(\theta_{h})=d_{i}^{\epsilon}(\theta_{h})\overline{d^{\epsilon}_{j}(\theta_{h})}$.

Then, we derive the product $T\widehat{\Sigma}^{\chi}_{ij,Raw}(\theta_{h})$ as

$\displaystyle Td_{i}^{\chi}(\theta_{h})\overline{d_{j}^{\chi}(\theta_{h})}$
$\displaystyle=\left(\sum_{t=1}^{T}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}u_{k,t-s}\mathrm{e}^{-\mathrm{i}\theta_{h}t}\right)\left(\sum_{t=1}^{T}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{jk,s}u_{k,t-s}\mathrm{e}^{\mathrm{i}\theta_{h}t}\right)$
$\displaystyle=\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}u_{k,t-s}u_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})},$

and the product $T\widehat{\Sigma}^{\epsilon}_{ij,Raw}(\theta_{h})$ as

$\displaystyle
Td_{i}^{\epsilon}(\theta_{h})\overline{d_{j}^{\epsilon}(\theta_{h})}$
$\displaystyle=\left(\sum_{t=1}^{T}\sum_{s=0}^{\infty}\sum_{k=1}^{r}c_{ik,s}e_{k,t-s}\mathrm{e}^{-\mathrm{i}\theta_{h}t}\right)\left(\sum_{t=1}^{T}\sum_{s=0}^{\infty}\sum_{k=1}^{r}c_{jk,s}e_{k,t-s}\mathrm{e}^{\mathrm{i}\theta_{h}t}\right)$
$\displaystyle=\sum_{s=0}^{\infty}\sum_{k=1}^{r}c_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}c_{jk^{\prime},s^{\prime}}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}e_{k,t-s}e_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}.$

Similarly, at each $\theta_{h}$, $h\leq|M_{T}|$, we define the common
component of the input spectral density matrix
$\widetilde{\Sigma}(\theta_{h})$ (see (2)) as

$\displaystyle\widetilde{\Sigma}_{ij}(\theta_{h})^{\chi}$
$\displaystyle=\frac{1}{2\pi}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\widehat{\Sigma}_{ij,Raw}^{\chi}(\theta_{h})$
$\displaystyle=\frac{1}{2\pi}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\frac{1}{T}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}u_{i,t-s}u_{j,t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})},$

and the idiosyncratic component of the input $\widetilde{\Sigma}(\theta_{h})$
as

$\displaystyle\widetilde{\Sigma}_{ij}(\theta_{h})^{\epsilon}$
$\displaystyle=\frac{1}{2\pi}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\widehat{\Sigma}_{ij,Raw}^{\epsilon}(\theta_{h})$
$\displaystyle=\frac{1}{2\pi}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\frac{1}{T}\sum_{s=0}^{\infty}\sum_{k=1}^{p}c_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{p}c_{jk^{\prime},s^{\prime}}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}e_{i,t-s}e_{j,t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})},$

where $K^{*}(\theta_{h})$ is the spectral window associated to the lag window
$K\left(\frac{k}{M_{T}}\right)$ satisfying the kernel assumption (see e.g.
Priestley,  (1981)).

We can now apply the framework of Wu and Zaffaroni,  (2018), disentangling the
common and the idiosyncratic components of the input estimator. Let
$P_{tt^{\prime},ij}^{\chi}=\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}u_{i,t-s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}u_{j,t^{\prime}-s^{\prime}}$
and
$P_{tt^{\prime},ij}^{\epsilon}=\sum_{s=0}^{\infty}\sum_{k=1}^{p}c_{ik,s}e_{i,t-s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{p}c_{jk^{\prime},s^{\prime}}e_{j,t^{\prime}-s^{\prime}}$.

Then

$\displaystyle\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}|\widetilde{\Sigma}_{ij}(\theta_{h})-{\Sigma}_{ij}(\theta_{h})|^{2}$
(16) $\displaystyle=\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\frac{1}{2\pi
T}\Bigg{|}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}(P_{tt^{\prime},ij}^{\chi}+P_{tt^{\prime},ij}^{\epsilon})\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-({\sigma}_{ij}(\theta_{h})^{\chi}+{\sigma}_{ij}(\theta_{h})^{\epsilon})\Bigg{|}^{2}$
$\displaystyle=\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\frac{1}{2\pi
T}\Bigg{|}\left(\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}P_{tt^{\prime},ij}^{\chi}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-{\sigma}_{ij}(\theta_{h})^{\chi}\right)$
$\displaystyle+\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\left(\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}P_{tt^{\prime},ij}^{\epsilon}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-{\sigma}_{ij}(\theta_{h})^{\epsilon}\right)\Bigg{|}^{2}$
$\displaystyle\leq\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\frac{1}{2\pi
T}\Bigg{|}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}P_{tt^{\prime},ij}^{\chi}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-{\sigma}_{ij}(\theta_{h})^{\chi}\Bigg{|}^{2}$
$\displaystyle+\Bigg{|}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}P_{tt^{\prime},ij}^{\epsilon}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-{\sigma}_{ij}(\theta_{h})^{\epsilon}\Bigg{|}^{2}$

which in turn is equal to

$\displaystyle=\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\Bigg{|}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}$
$\displaystyle\left[\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}u_{k,t-s}u_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\right]\Bigg{|}^{2}$
(17) $\displaystyle+\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\Bigg{|}\sum_{s=0}^{\infty}\sum_{k=1}^{p}c_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{p}c_{jk^{\prime},s^{\prime}}$
$\displaystyle\left[\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}e_{k,t-s}e_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\right]\Bigg{|}^{2}$
(18)

Under Assumption 3.8 and the kernel assumption, due to Theorem 1 in (Wu and
Zaffaroni, , 2018, Section 4.2), which states that

$\max_{1\leq i,j\leq p}\mathrm{E}\max_{h\leq|M_{T}|}\Bigg{|}\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}u_{k,t-s}u_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\Bigg{|}^{2}=O\left(\frac{M_{T}\log(M_{T})}{T}\right),$

for the common component of (16), which is (17), it holds

$\displaystyle\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}(\Bigg{|}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}$
$\displaystyle\left[\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}u_{k,t-s}u_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\right]\Bigg{|}^{2}$
$\displaystyle\leq\Bigg{\|}\sum_{s}B_{s}\Bigg{\|}_{\infty,v}^{2}\;\Bigg{\|}\sum_{s}B_{s}\Bigg{\|}_{\infty,v}^{2}\frac{M_{T}\log(M_{T})}{T}=O\left(\frac{M_{T}\log(M_{T})}{T}\right).$
(19)

Similarly, due to the same reasons, for the idiosyncratic component of (16),
which is (18), it holds

$\displaystyle\max_{1\leq i,j\leq
p}\mathrm{E}\max_{h\leq|M_{T}|}\Bigg{|}\sum_{s=0}^{\infty}\sum_{k=1}^{r}c_{ik,s}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}c_{jk^{\prime},s^{\prime}}$
$\displaystyle\left[\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}e_{k,t-s}e_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\right]\Bigg{|}^{2}$
$\displaystyle\leq\Bigg{\|}\sum_{s}C_{s}\Bigg{\|}_{\infty,v}^{2}\;\Bigg{\|}\sum_{s}C_{s}\Bigg{\|}_{\infty,v}^{2}\frac{M_{T}\log(M_{T})}{T}=O\left(\frac{M_{T}\log(M_{T})}{T}\right).$
(20)

Let us now consider the expected maximum overall Frobenius loss of the input
across frequencies:

$\displaystyle\mathrm{E}\max_{h\leq|M_{T}|}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{F}=\sqrt{\sum_{i=1}^{p}\sum_{j=1}^{p}(\widehat{\Sigma}_{ij}(\theta_{h})-\Sigma_{ij}(\theta_{h}))^{2}}$
$\displaystyle\leq\mathrm{E}\max_{h\leq|M_{T}|}\sqrt{\sum_{i=1}^{p}\sum_{j=1}^{p}(\widehat{\Sigma}_{ij}^{\chi}(\theta_{h})-\Sigma_{ij}^{\chi}(\theta_{h}))^{2}}+\mathrm{E}\max_{h\leq|M_{T}|}\sqrt{\sum_{i=1}^{p}\sum_{j=1}^{p}(\widehat{\Sigma}_{ij}^{\epsilon}(\theta_{h})-\Sigma_{ij}^{\epsilon})^{2}}.$
(21)

For the first term in the rhs of inequality (A) it holds

$\displaystyle\mathrm{E}\max_{h\leq|M_{T}|}\Bigg{|}\sum_{i=1}^{p}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\sum_{j=1}^{p}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}$
$\displaystyle\left[\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}u_{k,t-s}u_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\right]\Bigg{|}^{2}$
$\displaystyle\leq\Bigg{|}\sum_{i=1}^{p}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\sum_{j=1}^{p}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}\Bigg{|}^{2}G\left(\frac{M_{T}\log(M_{T})}{T}\right)$
$\displaystyle\leq\Bigg{|}\sum_{i=1}^{p}\sum_{s=0}^{\infty}\sum_{k=1}^{r}b_{ik,s}\Bigg{|}^{2}\;\Bigg{|}\sum_{j=1}^{p}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{r}b_{jk^{\prime},s^{\prime}}\Bigg{|}^{2}G\left(\frac{M_{T}\log(M_{T})}{T}\right)$
$\displaystyle\leq{r^{2}}\Bigg{|}\Bigg{|}\sum_{s}B_{s}\Bigg{|}\Bigg{|}_{1,v}^{2}\;\Bigg{|}\Bigg{|}\sum_{s}B_{s}\Bigg{|}\Bigg{|}_{1,v}^{2}G\left(\frac{M_{T}\log(M_{T})}{T}\right)$
$\displaystyle\leq
Gp^{2\alpha}\left(\frac{M_{T}\log(M_{T})}{T}\right)=O\left(\frac{p^{2\alpha}M_{T}\log(M_{T})}{T}\right)$
(22)

because Assumptions 3.8 and the kernel assumption hold.

Similarly, due to the same reasons, for the second term on the rhs of equation
(A) it holds:

$\displaystyle\mathrm{E}\max_{h\leq|M_{T}|}\Bigg{|}\sum_{i=1}^{p}\sum_{s=0}^{\infty}\sum_{k=1}^{p}c_{ik,s}\sum_{j=1}^{p}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{p}c_{jk^{\prime},s^{\prime}}$
$\displaystyle\left[\frac{1}{2\pi
T}\sum_{h=-M_{T}}^{M_{T}}K^{*}(\theta_{h})\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}e_{k,t-s}e_{k^{\prime},t^{\prime}-s^{\prime}}\mathrm{e}^{-\mathrm{i}\theta_{h}(t-t^{\prime})}-\mathrm{e}^{-\mathrm{i}\theta_{h}(s-s^{\prime})}\right]\Bigg{|}^{2}$
$\displaystyle\leq\Bigg{|}\sum_{i=1}^{p}\sum_{s=0}^{\infty}\sum_{k=1}^{p}c_{ik,s}\sum_{j=1}^{p}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{p}c_{jk^{\prime},s^{\prime}}\Bigg{|}^{2}G\left(\frac{M_{T}\log(M_{T})}{T}\right)$
$\displaystyle\leq\Bigg{|}\sum_{i=1}^{p}\sum_{s=0}^{\infty}\sum_{k=1}^{p}c_{ik,s}\Bigg{|}^{2}\;\Bigg{|}\sum_{j=1}^{p}\sum_{s^{\prime}=0}^{\infty}\sum_{k^{\prime}=1}^{p}c_{jk^{\prime},s^{\prime}}\Bigg{|}^{2}G\left(\frac{M_{T}\log(M_{T})}{T}\right)$
$\displaystyle\leq\Bigg{|}\Bigg{|}\sum_{s}C_{s}\Bigg{|}\Bigg{|}_{1,v}^{2}\;\Bigg{|}\Bigg{|}\sum_{s}C_{s}\Bigg{|}\Bigg{|}_{1,v}^{2}G\left(\frac{M_{T}\log(M_{T})}{T}\right)$
$\displaystyle\leq
G_{1}p^{2\delta^{\prime}}\left(\frac{M_{T}\log(M_{T})}{T}\right)=O\left(\frac{p^{2\delta^{\prime}}M_{T}\log(M_{T})}{T}\right).$
(23)

Since $\delta^{\prime}<\alpha$, from (A) we can derive that
$\mathrm{E}\max_{h\leq|M_{T}|}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{F}\leq
O\left(p^{\alpha}\sqrt{\frac{M_{T}\log(M_{T})}{T}}\right)$, from which the
thesis
$\mathrm{E}\max_{h\leq|M_{T}|}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}\leq
O\left(p^{\alpha}\sqrt{\frac{M_{T}\log(M_{T})}{T}}\right)$ follows.

#### Proof of Theorem 4.1

Following Luo, 2011a , we note that under Assumptions 3.4 and 3.6, setting
$\psi=\frac{p^{\alpha}}{\xi(T)}\frac{1}{\sqrt{T}}$ with $\rho=\gamma\psi$
(where $\gamma\in[9\xi(T),1/(6\mu(\Omega))]$), and further assuming that
$\underline{\delta}_{T}p^{2(\alpha-\underline{\delta})}<T<\overline{\delta}_{T}p^{6\delta}$
for some $\underline{\delta}_{T},\overline{\delta}_{T}$ such that
$0<\underline{\delta}_{T}<\overline{\delta}_{T}$ and that the minimum
eigenvalue of $L^{*}$, $\lambda_{r}(L^{*})$, is larger than
$G_{2}\frac{\psi}{\xi^{2}(T)}$, Propositions 12, 13, and 14 in Luo, 2011a can
be directly applied to our setting, proving that, for each
$\theta_{h}=\frac{h\pi}{M_{T}}$, $h\in-[M_{T}],\ldots,[M_{T}]$, the pair
$(\widehat{L}(\theta_{h}),\widehat{S}(\theta_{h}))$ minimizing (3) satisfies
the following theses:

  1. i)

$g_{\gamma}(\widehat{S}(\theta_{h})-{S}(\theta_{h}),\widehat{L}(\theta_{h})-{S}(\theta_{h}))$
is upper bounded;

  2. ii)

$\widehat{L}(\theta_{h})$ is rank-consistent:
$\text{rk}(\widehat{L}(\theta_{h}))=\text{rk}({L}(\theta_{h}))$;

with probability depending on the random loss
$\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}$. In particular, the
upper bound on
$g_{\gamma}(\widehat{S}(\theta_{h})-{S}(\theta_{h}),\widehat{L}(\theta_{h})-{S}(\theta_{h}))$
depends on $\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}$.
Considering the frequency grid $\theta_{h}=\pi\frac{h}{M_{T}}$,
$h\leq|M_{T}|$, Lemma 4.1 states that under Assumptions 3.8 and the kernel
assumption it holds
$\max_{h}\|\widehat{\Sigma}(\theta_{h})-\Sigma(\theta_{h})\|_{2}\leq
O\left(p^{\alpha}\sqrt{\frac{M_{T}\log(M_{T})}{T}}\right)$, which leads,
setting $\psi=\frac{p^{\alpha}}{\xi(\mathcal{T})}\sqrt{\frac{M_{T}\log
M_{T}}{T}}$, to

$g_{\gamma}(\widehat{S}(\theta_{h})-{S}(\theta_{h}),\widehat{L}(\theta_{h})-{L}(\theta_{h}))\leq
C\psi.$ (24)

From (24), all the claims of Theorem 4.1 follow, because

$\displaystyle||\widehat{L}(\theta_{h})-{L}(\theta_{h})||_{2}$
$\displaystyle\leq$ $\displaystyle C\psi;$
$\displaystyle||\widehat{S}(\theta_{h})-{S}(\theta_{h})||_{\infty}$
$\displaystyle\leq$ $\displaystyle C\gamma\psi\leq\xi(T)\psi;$
$\displaystyle||\widehat{S}(\theta_{h})-{S}(\theta_{h})||_{2}$
$\displaystyle\leq$ $\displaystyle Cq^{\prime}\gamma\psi\leq
q^{\prime}\xi(T)\psi;$
$\displaystyle||\widehat{\Sigma}(\theta_{h})-{\Sigma}(\theta_{h})||_{2}$
$\displaystyle\leq$ $\displaystyle C\psi+q^{\prime}\xi(T)\psi.$

Unlike Assumptions 3.4-3.6 and the lower bound on $\lambda_{r}(L^{*})$, if the
lower bound on $S_{\text{\tiny{min,off}}}$ does not hold, there is no
consequence on the identification of the two underlying algebraic varieties,
i.e. on parametric consistency and rank recovery. The only consequence lies in
the fact that some nonzero elements of $\widehat{S}(\theta_{h})$ are not
recovered. This fact can be appreciated by looking at the proofs of
Propositions 5.2 and 5.3 in Chandrasekaran et al.,  (2012), directly exploited
by Luo, 2011a .

If, instead, the condition
$\|S\|_{\text{\tiny{min,off}}}>G_{3}\frac{\psi}{\mu(\Omega)}$ holds in
addition to all the assumptions and conditions of Theorem 4.1, the same
Propositions in Chandrasekaran et al.,  (2012) allow to conclude that the
recovered sparsity pattern is also consistent:
$\text{sgn}(\mathrm{Re}(\widehat{S}(\theta_{h})))=\text{sgn}(\mathrm{Re}({S}(\theta_{h})))$.

#### Proof of Theorem 4.2

Under the assumptions of Theorem 4.2, the assumptions and conditions of
Theorem 4.1 are automatically satisfied. While the assumptions on kernel
window and temporal dependence are exactly the same as in the basic filter
setting of Assumption 3.2, the assumptions on latent eigenvalues and residual
sparsity pattern are reshaped to cope with the general filters prescribed by
Assumption 3.3. The assumption needed to ensure the identifiability of
underlying algebraic varieties is also reshaped accordingly. Therefore, all
the claims of Theorem 4.1 are still valid, in a much more general context
where the latent eigenvalues and the residual sparsity pattern are
intermediately spiked, and the latent coefficient matrices may have different
condition numbers across frequencies while the residual coefficient matrices
may have different sparsity patterns. The GDFM setting is a special case into
this context, where the latent eigenvalues are spiked with $p$ and the maximum
row-wise number of residual nonzeros is bounded by a constant. Therefore,
Theorem 2.1 is proved as a special case of Theorem 4.2 with $\alpha=1$.

#### Proof of Corollary 4.1

Let us define
$\phi=G(\frac{p^{\alpha}}{\xi(T)}\frac{1}{\sqrt{T}}+q^{\prime}\frac{p^{\alpha}}{\sqrt{T}})\sqrt{\frac{M_{T}\log
M_{T}}{T}}$ and
$\phi_{S}=G(q^{\prime}\frac{p^{\alpha}}{\sqrt{T}})\sqrt{\frac{M_{T}\log
M_{T}}{T}}$. Suppose that
$\widehat{L}=\widehat{L}_{\text{\tiny{UNALSE}}}(\theta_{h})$,
$\widehat{S}=\widehat{S}_{\text{\tiny{UNALSE}}}(\theta_{h})$,
$\widehat{\Sigma}=\widehat{\Sigma}_{\text{\tiny{UNALSE}}}(\theta_{h})$, for
each $\theta_{h}=\pi\frac{h}{M_{T}}$, $h\leq|M_{T}|$ and
$\Sigma=\Sigma(\theta_{h})$ for each $\theta\in[-\pi,\pi]$. Weyl’s Theorem
prescribes that, for any matrix $\Sigma$, we have

$|\widehat{\lambda}_{i}-\lambda|\leq\|\widehat{\Sigma}-\Sigma\|_{2},\;\forall
i=1,\ldots,p,$

where $\widehat{\lambda}_{i}$, $i=1,\ldots,p$, are the sample eigenvalues.
This result relates the rate of sample eigenvalues to the matrix spectral loss
rate. The triangular inequality gives

$|\lambda_{p}(\widehat{L}+\widehat{S})-\lambda_{p}(\Sigma)|\leq|\lambda_{p}(\widehat{L}+\widehat{S})|+|-\lambda_{p}(\Sigma)|=|\lambda_{p}(\widehat{L}+\widehat{S})|+\lambda_{p}(\Sigma),$
(25)

because $\Sigma$ is positive definite. Thus,

$|\lambda_{p}(\widehat{L}+\widehat{S})|\geq|\lambda_{p}(\widehat{L}+\widehat{S})-\lambda_{p}(\Sigma)|-\lambda_{p}(\Sigma).$

Since for Weyl’s theorem
$|\lambda_{p}(\widehat{L}+\widehat{S})-\lambda_{p}(\Sigma)|\leq\phi$ we have

$\lambda_{p}(\widehat{L}+\widehat{S})>0\Longleftrightarrow\lambda_{p}(\Sigma)>\phi.$
(26)

This proves the first part of the claim.

In order to achieve the same rate $\phi$ for the inverse spectral rate
$\|(\widehat{L}+\leavevmode\nobreak\ \widehat{S})^{-1}-\leavevmode\nobreak\
\Sigma^{-1}\|_{2}$, it is necessary that $\lambda_{p}(\Sigma)\geq 2\phi$. In
fact, the triangular inequality gives

$\|(\widehat{L}+\widehat{S})^{-1}-\Sigma^{-1}\|_{2}\leq\|(\widehat{L}+\widehat{S})^{-1}\|_{2}+\lambda_{p}(\Sigma)^{-1}$
(27)

By summing and subtracting $\Sigma$ and using the triangular inequality

$\|(\widehat{L}+\widehat{S})^{-1}\|_{2}=\|(\widehat{L}+\widehat{S}-\Sigma+\Sigma)^{-1}\|_{2}\leq\|(\widehat{L}+\widehat{S}-\Sigma)^{-1}\|_{2}+\|\Sigma^{-1}\|_{2}\leq\|(\widehat{L}+\widehat{S}-\Sigma)^{-1}\|_{2}+\lambda_{p}(\Sigma)^{-1}.$
(28)

For the triangular inequality, we have

$|\lambda_{p}((\widehat{L}+\widehat{S})^{-1})-\lambda_{p}(\Sigma)^{-1}|\leq|\lambda_{p}((\widehat{L}+\widehat{S})^{-1})|+|-\lambda_{p}(\Sigma)^{-1}|\leq|\lambda_{p}((\widehat{L}+\widehat{S})^{-1})|+\lambda_{p}(\Sigma)^{-1}$
(29)

since $\Sigma$ is positive definite.

At the same time, we want that

$\|(\widehat{L}+\widehat{S})^{-1}-\Sigma^{-1}\|_{2}\leq\phi.$

Hence, inequality (27) becomes

$\phi^{-1}\leq|\lambda_{p}((\widehat{L}+\widehat{S})^{-1})|+2\lambda_{p}(\Sigma)^{-1}.$

We can write

$|\lambda_{p}((\widehat{L}+\widehat{S})^{-1})|\geq\phi^{-1}-2\lambda_{p}(\Sigma)^{-1},$

which allows to conclude that

$\|\widehat{\Sigma}^{-1}-\Sigma^{-1}\|_{2}\leq\phi\Longleftrightarrow\phi^{-1}\geq
2\lambda_{p}(\Sigma)^{-1}.$ (30)

Using (30), it is possible to derive the rate for
$(\widehat{L}+\widehat{S})^{-1}$. By the property (see Luo, 2011a , pp.
31-32):

$\|(M+N)^{-1}-M^{-1}\|_{2}\leq\|M^{-1}\|_{2}\cdot\|N\|_{2}\cdot\|(M+N)^{-1}\|_{2}$
(31)

we obtain

$\displaystyle\|(\widehat{L}+\widehat{S})^{-1}-\Sigma^{-1}\|_{2}=\|(\widehat{L}+\widehat{S})^{-1}[\widehat{L}+\widehat{S}-\Sigma]\Sigma^{-1}\|_{2}\leq$
$\displaystyle\leq\|(\widehat{L}+\widehat{S})^{-1}\|_{2}\cdot\|[\widehat{L}+\widehat{S}-\Sigma]\|_{2}\cdot\|\Sigma^{-1}\|_{2}\leq\frac{2}{\lambda_{p}(\Sigma)^{2}}\|[\widehat{L}+\widehat{S}-\Sigma]\|_{2}.$

Hence, we have

$\|\widehat{\Sigma}^{-1}-\Sigma^{-1}\|_{2}\leq G(q^{\prime}\xi(T)+1)\psi=\phi$
(32)

The same reasoning can be carried out for $\widehat{S}$ by simply replacing
$\phi$ by $\phi_{S}$, thus obtaining

$\|\widehat{S}^{-1}-S^{-1}\|_{2}\leq G(q^{\prime}\xi(T))\psi=\phi_{S}.$ (33)

##  Appendix B Admissible sparsity regimes

######  Remark B.1 (Admissible sparsity regimes).

Finally, we analyze in detail the admissible sparsity regimes for the residual
spectral density matrix $S(\theta)$ at a given frequency. The parameters
involved in this analysis are the latent eigenvalues rate $\alpha$, the
maximum number of nonzero elements per row, $q^{\prime}$, the minimum absolute
nonzero off-diagonal element $\|S(\theta)\|_{min,off}$, the row-wise maximum
$l_{1}$ norm $\|S(\theta)\|_{1,v}$, and the sample size $T$.

First of all, we note that the following inequality holds:

$q^{\prime}\|S(\theta)\|_{min,off}\leq\|S(\theta)\|_{1,v}.$ (34)

At the same time, from Theorems 4.1 and 4.2 we know that sparsistency requires
that $\|S(\theta)\|_{min,off}>\frac{\psi}{\mu(\Omega)}$, from which we can
write
$\|S(\theta)\|_{min,off}>\frac{p^{\alpha}}{\sqrt{T}}\frac{1}{\xi(\mathcal{T})\mu(\Omega)}$,
that becomes

$\|S(\theta)\|_{min,off}\sqrt{T}\gtrsim p^{\alpha}$ (35)

because $\xi(\mathcal{T})\mu(\Omega)=O(1)$.

From the assumptions of Theorems 4.1 and 4.2, we know that $q^{\prime}\lesssim
p^{\delta}$, $\delta<\alpha$, and that
$\underline{\delta}_{T}p^{2(\alpha-\underline{\delta})}<T<\overline{\delta}_{T}p^{6\delta}$
for some $\underline{\delta}_{T},\overline{\delta}_{T}$ such that
$0<\underline{\delta}_{T}<\overline{\delta}_{T}$. For some $\delta_{2}>0$, we
require $\|S(\theta)\|_{1,v}\leq\delta_{2}p^{\underline{\delta}}$, with
$\underline{\delta}\leq\delta+0.5$, $\underline{\delta}\leq\delta^{\prime}$,
$\underline{\delta}<\alpha$. Letting
$\|S(\theta)\|_{min,off}=O(p^{\vartheta})$ and $\sqrt{T}=O(p^{\iota})$, by
combining (34) and (35) we obtain

$p^{\alpha-\iota}<p^{\vartheta}\leq p^{\underline{\delta}-\delta}.$ (36)

As $\iota\leq\delta$, sparsistency is not possible, because (36) leads to
$\underline{\delta}\geq\alpha$. As $\iota>\delta$, i.e., as $T$ grows,
sparsistency becomes possible. If $\iota=\frac{3}{2}\delta$, sparsistency
requires for instance that $\alpha<\underline{\delta}+\frac{1}{2}\delta$. This
condition is compatible, among others, with the setting $\alpha=\frac{3}{4}$,
$\underline{\delta}=\frac{2}{3}$, $\delta=\frac{1}{3}$, because by (36)
$\iota=\frac{3}{2}\delta$ leads to $\alpha<\frac{5}{6}$.

In order to understand the admissible relative scalings of $\vartheta$,
$\delta$, $\underline{\delta}$, we need to study the sign of the quantity
$\alpha-\iota$. If for instance $\iota=\frac{3}{2}\delta$ and
$\delta=\frac{1}{3}$, we can notice that
$\alpha-\iota=\alpha-\frac{3}{2}\delta=\alpha-\frac{3}{2}\frac{1}{3}<0$ if
$\alpha<\frac{1}{2}$. Therefore, from (36) we know that the condition
$\underline{\delta}-\delta<0$ is admissible in that case. As an example,
setting $\delta=\frac{1}{3}$ and $\alpha=\frac{3}{8}$, we obtain from (36)
that $\frac{3}{8}-\frac{1}{2}<\underline{\delta}-\frac{1}{3}$, which means
$\underline{\delta}>-\frac{1}{8}+\frac{1}{3}=\frac{5}{24}$, smaller than
$\delta=\frac{1}{3}$ and $\alpha=\frac{3}{8}$.

To sum up, if for some $\varepsilon>0$ $T\gtrsim p^{\delta+\varepsilon}$,
i.e., if $\iota>\delta$, sparsistency is possible. If $\alpha-\iota>0$, it is
necessary that $\vartheta>0$, with $\underline{\delta}$ large and $\delta$
small to respect the condition
$\alpha-\iota<\vartheta<\underline{\delta}-\delta$ with
$\underline{\delta}<\alpha$ and $\delta<\alpha$. This situation corresponds to
a sparsity pattern with few large nonzero entries.

As $\iota$ increases, i.e. as the sample size $T$ grows, it is more likely
that $\alpha-\iota<0$. In that case, $\vartheta$ may be smaller than $0$, and
it may also hold $\underline{\delta}<\delta$. This corresponds to a sparsity
pattern with a larger number of smaller nonzero entries (than before).

Note that the former is usually a case with large $\alpha$, while the latter
requires a smaller $\alpha$. This implies that spiked latent eigenvalues
require large residual entries to achieve sparsistency, while smaller latent
eigenvalues require small residual entries to ensure latent rank recovery.

##  Appendix C Simulation mechanism

Our purpose is to obtain data with a low rank plus sparse spectral density
matrix. First, we fix the dimension $p$, the sample size $T$, the latent rank
$r$ and the condition number of $L^{*}$, $c$. Concerning basic filters, the
first problem we encounter regards the generation of matrices $U_{L}$ and
$\Lambda_{u}$. We apply the generation algorithm in Farné,  (2016) to simulate
a latent multivariate spectrum with fixed condition number $c$ and trace equal
to $\tau\beta p$, where $\beta$ is the latent variance proportion (constant
across frequencies) and $\tau$ is a scale parameter. The procedure consists in
applying the Gram-Schmidt algorithm to a permutation matrix, randomly drawing
$r$ random eigenvectors (to become the orthonormal columns of the matrix
$U_{L}$), and then building $\Lambda_{u}$ as a diagonal $r\times r$ matrix
with ordered equidistant elements such that
${\Lambda_{u}}_{1,1}/{\Lambda_{u}}_{r,r}=c$. We thus obtain
$L^{*}=U_{L}\Lambda_{u}U_{L}^{\prime}$. For given time coefficients
$\lambda_{s}$ such that $\sum_{s=0}^{n_{l}}\lambda_{s}^{2}=1$ (where $n_{l}$
is the chosen number of lags), we obtain the coefficient matrices accordingly
as $B_{s}=U_{L}\sqrt{\Lambda_{u}}\lambda_{s}$, which allow us to generate the
common component $\chi_{t}$, $t=1,\ldots,T$, as a vector moving average (VMA)
with $n_{l}$ lags.

Following Assumption 3.2, we generate the sparse component in the following
way. First, we generate the diagonal of $S^{*}$ from a Dirichlet distribution
with parameter $(1-\beta)\tau p$. Then, we re-order these residual variances
matching the respective magnitude order of the variances in $L^{*}$.
Exploiting Cauchy-Schwartz inequality, we then randomly generate the off-
diagonal elements $S_{ij}^{*}$ from a uniform distribution
$Unif(0,\delta\sqrt{S_{ii}^{*}S_{jj}^{*}})$, where $\delta$ is a tuning
parameter. The generated residual off-diagonal elements are then ordered, and
the survival threshold is set to their maximum times a proportion parameter
$\delta_{bis}$. The residual coefficient matrices are then obtained
accordingly as $C_{s}=U_{S}\sqrt{\Lambda_{e}}\lambda_{s}$, where
$U_{S}\Lambda_{e}U_{S}^{\prime}$ is the spectral decomposition of $S^{*}$,
allowing to generate the idiosyncratic component $\epsilon_{t}$,
$t=1,\ldots,T$, as a VMA as well.

Concerning general filters, we follow Assumption 3.3. We start from the basic
filters $B_{s}=U_{L}\sqrt{\Lambda_{u}}\lambda_{s}$. We then replace the scalar
$\lambda_{s}$ by a $r\times r$ diagonal matrix $D_{L,s}$, built as follows. We
set a perturbation proportion $\kappa_{pert}=0.1$. We generate the diagonal of
$D_{L,s}$ as
$\lambda_{s}\mathbbm{1}_{r}-\kappa_{pert}\lambda_{s}\mathbbm{1}_{r}+2\kappa_{pert}\lambda_{s}W_{r}$,
where $\mathbbm{1}_{r}$ is a vector composed by $r$ ones, and $W_{r}$ is a
vector of $r$ random numbers drawn from a uniform between $0$ and $1$. The
general filters of the low rank component are then obtained as
$B_{s}=U_{L}D_{L,s}\sqrt{\Lambda_{u}}$. This procedure generates a varying
latent variance proportion $\beta(\theta_{h})$ across chosen frequencies.

Concerning the residual filters, our procedure generates the diagonals of
$\Gamma_{\epsilon}(s)$, $s=0,1,\ldots,n_{l}$, from a Dirichlet distribution
with parameter $(1-\beta)\tau p\lambda_{s}$. We then order the elements of
each of the $n_{l}+1$ diagonals matching the respective magnitude order of the
variances in $\Gamma_{\chi}(0)$. Starting from those diagonals, we repeat the
above exposed thresholding procedure $n_{l}+1$ times. For each of the
subsequent sparsified $\Gamma_{\epsilon}(s)$, we derive the spectral
decomposition $\Gamma_{\epsilon}(s)=U_{S,s}\Lambda^{2}_{S,s}U_{S,s}^{\prime}$,
and we set the generalized filters as $C_{s}=U_{S,s}\Lambda_{S,s}$. Note that
this passage is very delicate, as the matrices of eigenvectors may lead to
non-sparse residual spectra if the nonzeros are too many, due to rounding
errors. Anyway, as prescribed in Assumption 3.3, each generated
$\Gamma_{\epsilon}(s)$ has in principle a different sparsity pattern and a
different number of nonzeros $q_{s}$ due to the Cauchy-Schwartz inequality.

The spectral density matrices at each frequency $\theta_{h}$ are then
computed. We set $|h|\leq M_{T}$, $\theta_{h}=\frac{h\pi}{M_{T}}$, and we then
calculate the low rank and residual transfer function matrices as
$B(\theta_{h})=\sum_{s=0}^{n_{l}}B_{s}\mathrm{e}^{(-is\theta_{h})}$ and
$C(\theta_{h})=\sum_{s=0}^{n_{l}}C_{s}\mathrm{e}^{(-is\theta_{h})}$. The
latent and residual spectral density matrices are thus derived at each
frequency $\theta_{h}$ as
$L(\theta_{h})=B(\theta_{h})\Lambda_{u}B(\theta_{h})^{\prime}$ and
$S(\theta_{h})=C(\theta_{h})\Lambda_{e}C(\theta_{h})^{\prime}$. The basic
filter specifications lead to real spectra. On the contrary, the general
filters lead to complex spectra, as the residual spectral components are
complex.

Once we have generated the common component $\chi_{t}$ and the residual
component $\epsilon_{t}$ as in Section D, we can generate our series $x_{t}$,
$t=1,\ldots,T$, according to equation (5). A relevant choice to control the
spectral shape lies in the vector moving average coefficients $\lambda_{s}$,
$s=0,\ldots,n_{l}$. We observe that the spectral shapes across frequencies are
characterized by the following pattern:

  * •

$n_{l}=1$, a positive $\lambda_{0}$ and negative $\lambda_{1}$ lead to a
"reverse S-shape";

  * •

$n_{l}=1$, a negative $\lambda_{0}$ and positive $\lambda_{1}$ lead to a "
S-shape";

  * •

$n_{l}=2$, a positive $\lambda_{0}$, $\lambda_{1}=0$ and a negative
$\lambda_{2}$ lead to a "U-shape";

  * •

$n_{l}=2$, a negative $\lambda_{0}$, $\lambda_{1}=0$ and a positive
$\lambda_{2}$ lead to a "reverse U-shape";

We note that when the non-null coefficients are equal to $0.5$, the spectral
shape has the maximum variability possible across frequencies, i.e., the
spectral peak attains its maximum and the spectral drop is $0$. Any situation
with any of the coefficients equal to $1$ leads instead to a constant spectrum
across frequencies. In the simulation study of Section 6, we set
$\lambda_{0}=0.8$ and $\lambda_{1}=0.2$. We select as our target the reverse
S-shape without loss of generality, because the conditions of Theorems 4.1 and
4.2 must be satisfied frequency-wise, and the reverse S-shape across
frequencies is by far the most common in real time series.

##  Appendix D Additional simulation results

###  D.1 Scenario A

Figure 1: Estimated latent variance proportion $\widehat{\beta}(f_{h})$ \-
Scenario A.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 2: Positive predictive value $ppv(f_{h})$ \- Scenario A.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 3: Negative predictive value $npv(f_{h})$ \- Scenario A.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 4: $err_{\widehat{L}}(f_{h})$ \- Scenario A.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 5: $err_{ratio}(f_{h})$ \- Scenario A.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

###  D.2 Scenario B

Figure 6: Estimated latent variance proportion $\widehat{\beta}(f_{h})$ \-
Scenario B.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 7: Nonzero predictive values $nzpv(f_{h})$ \- Scenario B. Setting 3 | Setting 4 | Setting 5
---|---|---
|  |
Figure 8: $err_{\widehat{L}}(f_{h})$ \- Scenario B.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 9: $err_{ratio}(f_{h})$ \- Scenario B.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

###  D.3 Scenario C

Figure 10: Estimated latent variance proportion $\widehat{\beta}(f_{h})$ \-
Scenario C.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 11: $err_{\widehat{L}}(f_{h})$ \- Scenario C.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

Figure 12: $err_{ratio}(f_{h})$ \- Scenario C.

Setting 1 | Setting 2 | Setting 3
---|---|---
|  |

Setting 4 | Setting 5
---|---
|

## References

  * Altissimo et al., (2010) Altissimo, F., Cristadoro, R., Forni, M., Lippi, M., and Veronese, G. (2010).  New Eurocoin: Tracking economic growth in real time.  The Review of Economics and Statistics, 92(4):1024–1034.
  * Barigozzi and Brownlees, (2019) Barigozzi, M. and Brownlees, C. (2019).  NETS: Network estimation for time series.  Journal of Applied Econometrics, 34(3):347–364.
  * Barigozzi and Hallin, (2017) Barigozzi, M. and Hallin, M. (2017).  A network analysis of the volatility of high dimensional financial series.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 66(3):581–605.
  * Barigozzi et al., (2021) Barigozzi, M., Lippi, M., and Luciani, M. (2021).  Large-dimensional dynamic factor models: Estimation of impulse–response functions with $I(1)$ cointegrated factors.  Journal of Econometrics, 221(2):455–482.
  * Bickel and Levina, (2008) Bickel, P. J. and Levina, E. (2008).  Covariance regularization by thresholding.  The Annals of Statistics, 36(6):2577–2604.
  * Böhm and von Sachs, (2008) Böhm, H. and von Sachs, R. (2008).  Structural shrinkage of nonparametric spectral estimators for multivariate time series.  Electronic Journal of Statistics, 2:696–721.
  * Böhm and von Sachs, (2009) Böhm, H. and von Sachs, R. (2009).  Shrinkage estimation in the frequency domain of multivariate time series.  Journal of Multivariate Analysis, 100(5):913–935.
  * Breitung and Candelon, (2006) Breitung, J. and Candelon, B. (2006).  Testing for short-and long-run causality: A frequency-domain approach.  Journal of Econometrics, 132(2):363–378.
  * Brillinger, (2001) Brillinger, D. R. (2001).  Time Series: Data Analysis and Theory.  SIAM.
  * Cai et al., (2010) Cai, J.-F., Candès, E. J., and Shen, Z. (2010).  A singular value thresholding algorithm for matrix completion.  SIAM Journal on Optimization, 20(4):1956–1982.
  * Chandrasekaran et al., (2012) Chandrasekaran, V., Parrilo, P. A., and Willsky, A. S. (2012).  Latent variable graphical model selection via convex optimization.  The Annals of Statistics, 40(4):1935–1967.
  * Chandrasekaran et al., (2011) Chandrasekaran, V., Sanghavi, S., Parrilo, P. A., and Willsky, A. S. (2011).  Rank-sparsity incoherence for matrix decomposition.  SIAM Journal on Optimization, 21(2):572–596.
  * Chaudhuri and Lo, (2015) Chaudhuri, S. E. and Lo, A. W. (2015).  Spectral analysis of stock-return volatility, correlation, and beta.  In 2015 IEEE Signal Processing and Signal Processing Education Workshop (SP/SPE), pages 232–236. IEEE.
  * Corbae et al., (2002) Corbae, D., Ouliaris, S., and Phillips, P. C. (2002).  Band spectral regression with trending data.  Econometrica, 70(3):1067–1109.
  * (15) Dahlhaus, R. (2000a).  Graphical interaction models for multivariate time series.  Metrika, 51(2):157–172.
  * (16) Dahlhaus, R. (2000b).  A likelihood approximation for locally stationary processes.  The Annals of Statistics, 28(6):1762–1794.
  * Daubechies et al., (2004) Daubechies, I., Defrise, M., and De Mol, C. (2004).  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Communications on Pure and Applied Mathematics, 57(11):1413–1457.
  * Davis et al., (2016) Davis, R. A., Zang, P., and Zheng, T. (2016).  Sparse vector autoregressive modeling.  Journal of Computational and Graphical Statistics, 25(4):1077–1096.
  * Donoho, (2006) Donoho, D. L. (2006).  For most large underdetermined systems of linear equations the minimal $l_{1}$ norm solution is also the sparsest solution.  Communications on Pure and Applied Mathematics, 59(6):797–829.
  * Eichler, (2007) Eichler, M. (2007).  Granger causality and path diagrams for multivariate time series.  Journal of Econometrics, 137(2):334–353.
  * Fan et al., (2013) Fan, J., Liao, Y., and Mincheva, M. (2013).  Large covariance estimation by thresholding principal orthogonal complements.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 75(4):603–680.
  * Farné, (2016) Farné, M. (2016).  An algorithm to simulate VMA processes having a spectrum with fixed condition number.  Communications in Statistics-Simulation and Computation, 45(5):1664–1675.
  * Farnè and Montanari, (2020) Farnè, M. and Montanari, A. (2020).  A large covariance matrix estimator under intermediate spikiness regimes.  Journal of Multivariate Analysis, 176:104577.