The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
164,710 characters
Generalized Dynamic Factor Models and Volatilities: Consistency, Rates, and Prediction Intervals
\title {\textbf{\Large{Generalized Dynamic Factor Models and Volatilities:\\ \Large{Consistency, Rates, and Prediction Intervals}} } }
\author {{\sc Matteo Barigozzi$^{\dag}$ \qquad Marc Hallin$^{\ddag}$}\vspace{3mm} \\ \normalsize
$^{\dag}${\it LSE, Department of Statistics, Houghton Street, London WC2A 2AE, UK.} \\ \normalsize
E-Mail: [email removed]\vspace{2mm} \\ \normalsize
$^{\ddag}$ {\it ECARES, Universit\'e libre de Bruxelles CP114/4 B-1050 Bruxelles, Belgium. }\\ \normalsize
E-Mail: [email removed]
}
\date{\small{\today}}
\maketitle
\begin{abstract}
Volatilities, in high-dimensional panels of economic time series with a dynamic factor structure on the levels or returns,
typically also admit a dynamic factor decomposition. We consider a two-stage dynamic factor model method recovering the common and idiosyncratic components of both levels and log-volatilities. Specifically, in a first estimation step, we extract the common and idiosyncratic shocks for the levels, from which a log-volatility proxy is computed. In a second step, we estimate a dynamic factor model, which is equivalent to a multiplicative factor structure for volatilities, for the log-volatility panel. By exploiting this two-stage factor approach, we build one-step-ahead conditional prediction intervals for large $n\times T$ panels of returns. Those intervals are based on empirical quantiles, not on conditional variances; they can be either equal- or unequal-tailed. We provide uniform consistency and consistency rates results for the proposed estimators as both $n$ and~$T$ tend to infinity. We study the finite-sample properties of our estimators by means of Monte Carlo simulations. Finally, we apply our methodology to a panel of asset returns belonging to the S\&P100 index in order to compute one-step-ahead conditional prediction intervals for the period 2006-2013. A comparison with the componentwise GARCH benchmark (which does not take advantage of cross-sectional information) demonstrates the superiority of our approach, which is genuinely multivariate (and high-dimensional), nonparametric, and model-free. \\
\\
\noindent {\itshape JEL Classification}: C32, C38, C58.\\
\noindent
{\itshape Keywords}: Volatility, Dynamic Factor Models, Prediction intervals, GARCH.
\end{abstract}
\symbolfootnote[0]{\\ We thank Christian Brownlees, Christian Francq, and Haeran Cho for helpful comments. This paper was also presented at: ``Panel Data Forecasting Conference'', University of Southern California, Dornsife, Los Angeles, April, 2019; the ``6th Rimini Centre for Economic Analysis (RCEA) Time Series Econometrics Workshop'', University of Cyprus, Larnaca, June 2019, the ``International Association for Applied Econometrics (IAAE) 2019 Annual Conference'', University of Cyprus, Nicosia, June 2019, and the ``Workshop on High-Dimensional Data Analysis", Durham University, June
2019.
\noindent
}
\section{Introduction}
Data in high dimension unquestionably constitute one of the main challenges of contemporary statistics/econo\-metrics, and have become pervasive in most domains related with data sciences.
Time series have not escaped that evolution, and the analysis of high-dimensional time series---equivalently, large cross-sections of univariate time series or panels---today ranks among the most active topics in theoretical and applied econometrics.
The most successful methods so far in the analysis and prediction of high-dimensional time
series are based on the so-called factor model
approach. That approach, under its various forms, is based on a (non-observed) decomposition
of the observation (a large cross-section of time series with complex interrelations) into the sum of two mutually orthogonal (all
leads, all lags) components: the {\it common component}, driven by a small number of {\it factors} or
{\it common shocks}, and an {\it idiosyncratic component}, with some variations in the definitions of
``common'' and ``idiosyncratic,'' and the assumptions made. Regardless of the definition adopted, the common and idiosyncratic components typically are disentangled by means of adequate cross-sectional and/or temporal aggregation of the observed time series.
Those aggregation and factor model approaches are strongly rooted in the multivariate time-series methods developed in the eighties and nineties, of which George Tiao and his collaborators have been most influential and unremittable pioneers: see, for instance, \citet{tiao1972asymptotic}, \citet{tiao1978some}, \citet{tiao1980forecasting}, \citet{tiao1981modeling}, \citet{tsay1985use}, \citet{pena1987identifying}, and \citet{tiao1989model}.
The type of factor model we are considering here is the General or Generalized Dynamic Factor Model (GDFM) introduced by \citet{FHLR00}, which, by taking into account all leading and lagging linear dependencies among the data, encompasses most other models, as e.g. the static factor approaches by \citet{baing02}, \cite{stockwatson2002}, and \citet{FLM13}. Moreover, as emphasised in \citet{fornilippi01} and \citet{hallinlippi13}, beyond the usual assumptions of second-order stationarity and existence of spectral densities, the GDFM decomposition into a common and an idiosyncratic component basically does not place any structural constraints on the data-generating process. In this sense, contrary to static factor approaches, it is canonical, nonparametric and model-free. In this paper, we consider the one-sided GDFM estimation method recently described in \citet{FHLZ15,FHLZ17}.
Prediction, in classical univariate and moderately multivariate time series analysis, is an obvious and natural objective; it is certainly no less crucial in high dimension. Efficient prediction, however, should exploit the amount of information available, due to the complex cross-dependencies among the many cross-sectional components, in the present and lagged values of the whole cross-section; the larger the cross-section (i.e., the higher the dimension), the more crucial the role of that information, and the more delicate its recovering. Factor models naturally have been used in the construction of {\it point-predictors}, and quite successfully so: see, e.g., \citet{stockwatson2002}, \citet{baing08JoE}, \citet{FGLS18}, to quote only a very few. Those authors, however, are dealing, mostly, with macroeconomic data, while less attention has been given to factor model methods in the analysis and prediction of financial returns: see, e.g. \citet{chamberlainrotshild83}, \citet{CK93}, or \citet{ait2017}.
In particular, when dealing with returns, due to the presence of conditional distribution heterogeneity (of which conditional heteroskedasticity is only a very particular case), conditional volatility phenomenons are essential, and definitely should be taken into account when building conditional prediction limits or conditional prediction intervals.
Most multivariate methods available in the literature for the analysis of conditional heterogeneity are restricted to the study of conditional heteroskedasticity, and rely on parametrisations of the ARCH-GARCH or Stochastic Volatility type: see, for instance, the reviews by \citet{BLR06} and \citet{AMY06}. Because of the curse of dimensionality, however, only the very simplest models can be considered in high-dimensional panels, possibly inducing a nonnegligible loss of efficiency. Among those, the factor GARCH approach is the most popular, see e.g. \citet{DN89}, \citet{ENR92}, \citet{HRS92}, and \citet{SCF08}. {Static factor models directly based on volatilities have also been considered, but these fail to exploit the information contained in the idiosyncratic components of returns, see e.g. \citet{CKL06} and \citet{fan15}.} For these reasons, \citet{barigozzihallin15a} introduce a two-step GDFM approach by which the nonparametric and model-free virtues of factor models are used in a joint analysis of returns and volatilities. In \citet{barigozzihallin15b}, that two-step GDFM is combined with a GARCH strategy in order to produce point-forecasts for volatilities (see also \citealp{Trucios19} for a recent example), while \citet{barigozzihallin15c} and \citet{BHS18} apply the same methodology in a study of the dynamic interdependencies of US and international financial markets. A two-stage factor approach similar to ours but in a static factor model setting is proposed in \citet{CB15}.
The objective of this paper is to combine the same two-step GDFM approach with a quantile-based construction of conditional confidence limits producing conditional interval predictions rather than point-forecasts for returns. That objective requires nontrivial consistency results on the two-step GDFM estimation method, which are not provided in \citet{barigozzihallin15a,barigozzihallin15b,barigozzihallin15c}. The first part of this paper, therefore, is devoted to a careful asymptotic analysis of the two-step GDFM. We then describe the quantile-based construction of conditional confidence limits, which we apply to a dataset of S\&P100 daily returns.
The paper is organised as follows. In Section \ref{sec:mod}, we present the GDFM model for the stochastic processes of returns (levels) and log-volatilities, and give sufficient conditions for its existence and identification. Section~\ref{sec:est_summary} describes the estimation of the model, and Section \ref{sec:ap} establishes the consistency properties (with rates) of the proposed estimators. In Section \ref{eq:int}, we define the one-step-ahead conditional prediction confidence limits and intervals. In Section \ref{sec:sim}, we study the finite-sample properties of our estimators via simulations.
Section \ref{sec:emp} applies our methodology to a panel of daily returns of stocks listed in the S\&P100 index and investigates the resulting coverage performance. In Section \ref{sec:conc}, we conclude. Proofs are postponed to an Appendix.
\subsection*{Notation}
The sub-exponential norm of a scalar random variable~$X$ is defined as $\Vert X\Vert_{\psi_1}:=\sup_{p\ge 1} p^{-1}\mathrm E[|X|^p]^{1/p}$ (see e.g.~Definition 5.13 in \citealp{vershynin12}). The transposed complex conjugate of a complex vector $\bf p$ is denoted as ${\bf p}^\dag$ and $\Vert \bf p\Vert=\bf p^\dag\bf p$. For an hermitian complex $n\times n$ matrix $\mbf A$ with generic $(i,j)$ entry $a_{ij}$ and largest (in modulus) eigenvalue~$\mu^{\mbf A}_1$, let~$\Vert \mbf A\Vert_1:=\max_{j=1,\ldots,n} \sum_{i=1}^n |a_{ij}|$ and $\Vert \mbf A\Vert:={\mu_1^{\mbf A}}$.
As usual, $L$ stands for the lag operator, such that, given a stochastic vector process $\{\mbf Y_t | t\in\mathbb Z\}$, $L^k\mbf Y_t:=\mbf Y_{t-k}$ for any integer $k$ and any~$t\in\mathbb Z$. Last, we denote by $\mathbb I(\mathcal A)$ the indicator function of an event $\mathcal A$.
\section{A General Dynamic Factor Model for levels and volatilities}\label{sec:mod}
We throughout assume that all stochastic variables in this
paper belong to the Hilbert space $L_2(\Omega, \mathcal F , \mathrm P)$, \linebreak where~$(\Omega, \mathcal F , \mathrm P)$ is
some common probability space. We study double-indexed stochastic processes of the\linebreak form~$\mbf Y\!:= \{Y_{it} \vert i\in\mathbb{N} , \ t\in\mathbb{Z}\}$, with $n$-dimensional sub-processes $\mbf Y_n\!:= \{Y_{it} \vert i=~\!1,\ldots,n,~t\in~\!\mathbb{Z}\}$, $n\in\mathbb{N}$. In practice, we deal with the finite observed $n\times T$ realisation
$${\bf Y}_{n,T}:=\left(\begin{array}{cccc}
Y_{11}, & Y_{12}, & \ldots ,
& Y_{1T} \\
\vdots&\vdots&&\vdots \\
Y_{n1}, & Y_{n2}, & \ldots ,
& Y_{nT}
\end{array}\right)
$$
of $\mbf Y$. In the empirical application of Section \ref{sec:emp}, the $Y_{it}$'s are observed values of daily stock returns, and we therefore call $\mbf Y$ the ``levels'' process. The assumptions in Section \ref{sec:mod_level} are mainly taken from \citet{FHLZ17}, with some modifications, mostly concerning the idiosyncratic components.
On the other hand, the assumptions in Section \ref{sec:mod_vol} are new and are related to the log-volatility proxies originally introduced in \citet{barigozzihallin15a,barigozzihallin15b}.
\subsection{Model and assumptions for levels}\label{sec:mod_level}
The Generalized Dynamic Factor Model (GDFM) for the levels process $\mbf Y$ is a decomposition of $Y_{it}$ into
\begin{equation}
Y_{it}-\mathrm E[Y_{it}]= X_{it}+ Z_{it}, \quad i\in\mathbb{N} , \ t\in\mathbb{Z}\label{GDFM}
\end{equation}
with
\begin{equation}
X_{it} = \sum_{j=1}^q \sum_{k=0}^{\infty} b_{ijk} u_{jt-k}=\mbf b_i'(L)\mbf u_t\quad\text{and}\quad Z_{it}=\sum_{k=0}^{\infty} d_{ik} v_{it-k} = d_i(L) v_{it},\label{eq:idio_level}
\end{equation}
where $\mathrm E[Y_{it}]$ stands for the expected value of $Y_{it}$ and the processes $\mbf u:= \{u_{jt} \vert j=1,\ldots,q, \ t\in\mathbb{Z}\}$ \linebreak and~$\mbf v_{n}:= \{v_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\}$ are mutually orthogonal (at all leads and lags) $q$- and $n$-dimensional white noises, respectively. Call~$\mbf u$ the process of {\it common factors} or {\it common shocks} and~$\mbf v_n$ the process of {\it idiosyncratic shocks}; $X_{it}$ and~$Z_{it}$ are~$Y_{it}$'s {\it common} and {\it idiosyncratic components},
respectively.
Letting $\mbf X_n:= \{X_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\}$ and $\mbf Z_n:= \{Z_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\}$, equations~\eqref{eq:idio_level} in vector notation takes the form
\begin{align}\label{eq:common_lev_idio_vec}
\mbf X_{nt}=\mbf B_n(L)\mbf u_t,\qquad\mbf Z_{nt}=\mbf D_n(L)\mbf v_{nt}, \quad n\in\mathbb{N} , \ t\in\mathbb{Z}.
\end{align}
with $\mbf B_n(L):=(\mbf b_1(L)\ldots \mbf b_n(L))'$, and $\mbf D_n(L):=\text{\rm diag}(d_1(L)\ldots d_n(L))$. \medskip
More precisely, we assume that \eqref{GDFM}-\eqref{eq:idio_level} hold and satisfy the following assumptions:\vspace{-1mm}
\begin{assumption}[L1] $\,$
\begin{compactenum}[(i)]
\item the dimension $q$ of $\mbf u_t$ does not depend on $n$; the process $\mbf u :=\{\mbf u_t \vert t\in\mathbb{Z}\}$ is second-order white noise, with mean $\mbf 0_q$ and diagonal positive definite covariance $\bm{\Gamma}^{\rm u}$;
\item writing ${\mbf b}_{ik}:=(b_{i1k}\ldots b_{iqk})^\prime$ for the $q\times 1$ coefficient of $L^k$ in ${\mbf b}_i(L)$, there exists a constant~$M_1>0$ such that $\sum_{k=0}^{\infty}\Vert \mbf b_{ik}\Vert\, \vert k\vert\le M_1$ for all $i\in\mathbb{N}$;
\item the process $\mbf v := \{\mbf v_{nt} \vert t\in\mathbb{Z}\}$ is second-order white noise, with mean $\mbf 0_n$ and positive definite covariance~$\bm\Gamma_n^{\rm v}$; moreover, $\mathrm E[v_{it}|v_{is}]=0$ for all $i\in\mathbb N$ and $t,s\in\mathbb Z$ such that $t>s$;
\item there exists a constant $C_{\rm v}>0$ such that $\Vert \bm\Gamma_n^{\rm v}\Vert_1\le C_{\rm v}$ for all~$n\in\mathbb{N}$;
\item there exists a constant $M_2>0$ such that $\sum_{k=0}^{\infty}\vert d_{ik}\vert\, \vert k \vert\le M_2$ for all $i\in\mathbb{N}$;
\item $\text{\rm Cov}(u_{jt},v_{is})=0$ for all $i\in\mathbb N$, $j=1,\ldots, q$, and $t,s\in\mathbb Z$;
\item there exists a constant $M_3>0$ such that $\sum_{k_1,k_2,k_3\in\mathbb Z}
\big\vert\mathrm E[u_{j_1t}u_{j_2,t-k_1}u_{j_3,t-k_2}u_{j_4,t-k_3}]\big\vert\le M_3$ for all $j_1,j_2,j_3,j_4=1,\ldots,q$;
\item there exists a constant $M_4>0$ such that $\sum_{k_1,k_2,k_3\in\mathbb Z}
\big\vert\mathrm E[v_{i_1t}v_{i_2,t-k_1}v_{i_3,t-k_2}v_{i_4,t-k_3}]\big\vert \le M_4$ for all\linebreak $i_1,i_2,i_3,i_4\in\mathbb N$.\end{compactenum}
\end{assumption}
These assumptions are standard in the literature with the exception of part {\it (iv)} which imposes a mild form of sparsity on the covariance matrix of the idiosyncratic innovations. A similar condition can be found in \citet{FLM13} and is empirically verified by \citet{boivinng06} and \citet{baing08JoE} for US macroeconomic data, and by \citet{barigozzihallin15c} for stock returns. As a consequence of parts {\it (iv)} and {\it (v)}, the idiosyncratic components are allowed to be serially autocorrelated and mildly cross-correlated (see also Lemma \ref{lem:dyn_eval} below). Moreover, it is easy to check that such assumption is nesting other typical conditions on the cross-sectional dependence of idiosyncratic components (see e.g. \citealp{baing02}, and \citealp{stockwatson2002}, in the static factor model case). Parts {\it (ii)} and {\it (v)} imply absolute summability of the autocovariances and therefore the existence of a purely continuous spectral density.
Moreover these assumptions and existence of fourth-order moments in parts (vii) and (viii) are classical requirements for consistent estimation of the autocovariances and the spectral density (see e.g. Chapter IV, Theorem 6, in \citealp{hannan1970}, for the autocovariances, and the results in Section 6.2 in \citealp{priestley01}, and Theorem 5A in \citealp{parzen57}, for the spectral density). Last, in part {\it (iii)} we also make the typical assumption of martingale difference innovations used in the GARCH literature (see e.g. Definition 2.1 in \citealp{FZ11}).
It should be insisted, however, that the GDFM is not a {\it statistical model} in the usual sense, inasmuch as, beyond the requirement of second-order stationarity, the existence of a finite (but unspecified) $q$, and the existence of a spectrum, it does not really impose any restrictions on the data-generating process: as argued by \citet{fornilippi01} and \citet{hallinlippi13}, \eqref{GDFM}-\eqref{eq:idio_level} indeed constitute a representation result rather than a model equation.
On the filters $\mbf b_i(L)$ and $d_i(L)$ we furthermore impose the following assumptions:
\begin{assumption}[L2] $\,$
\begin{compactenum}[(i)]
\item $\mbf b_i(L)$ has rational entries, i.e. $b_{ij}(L)=\theta_{ij}(L)\phi^{-1}_{ij}(L)$, where $\phi_{ij}(z)$ and $\theta_{ij}(z)$, for all $i\in\mathbb N$ and \linebreak$j=1,\ldots, q$, are finite-order polynomials;
\item there exists a constant $\bar \phi>1$ such that $\phi_{ij}(z)\neq 0$ for all $i\in\mathbb N$, all $j=1,\ldots, q$, and all~$z\in\mathbb C$ such that $|z|\le \bar \phi$;
\item the coefficients $\theta_{ijk}$ of $\theta_{ij}(L)$ are such that $|\theta_{ijk}|\le B^X$ for some positive constant $B^X$, \linebreak all $k\in\mathbb N\cup \{0\}$, all~$i\in\mathbb N$, and $j=1,\ldots, q$;
\item $d_i(L)$ is of the form $c_i^{-1}(L)$ where $c_i(z)$, for all $i\in\mathbb N$, is a finite-order polynomial,~$c_{i}(0)=1$ and $c_{i}(z)\neq 0$ for all $z\in\mathbb C$ such that $|z|\le 1$.
\end{compactenum}
\end{assumption}
This latter assumption is not strictly needed and could be easily relaxed to allow for infinite order autoregressive dynamics---at the expense, however, of heavier notation and longer proofs; see also Section \ref{sec:ap} for a short discussion.
This assumption implies that both the common and idiosyncratic components have a rational spectral density. Rational filters for the common component are also assumed in \citet{FHLZ17}, while here we also assume that the idiosyncratic component admits a finite autoregressive representation. In particular, using part~{\it (iv)}, we can rewrite the second equation in \eqref{eq:idio_level} as
\begin{equation}
c_i(L) Z_{it} = v_{it}.\label{eq:var_lev_idio}
\end{equation}
Let $\bm\Sigma_n^Y(\theta)$, $\bm\Sigma_n^X(\theta)$ and $\bm\Sigma_n^Z(\theta)$, $\theta\in[-\pi,\pi]$, be the $n\times n$ spectral density matrices of the observed panel, the common, and the idiosyncratic components, respectively; the existence of those spectral densities is guaranteed by Assumption (L1). Denote by $\lambda_{nj}^Y(\theta)$, $\lambda_{nj}^X(\theta)$, and $\lambda_{nj}^Z(\theta)$ their respective~$j$-th largest eigenvalues---the {panel}, {common}, and idiosyncratic {\it dynamic eigenvalues}, on which we assume the following. {Hereafter, ``for\linebreak all~$\theta\in[-\pi,\pi]$'' or ``$\theta-a.e.$'' is to be understood as ``for all~$\theta$ but over a subset of values included in a set with Lebesgue measure zero.'' Similarly, $\sup_{\theta\in[-\pi,\pi]}$ in the sequel is an {\it essential} $\sup$, etc. \special{color cmyk 0 0 0 1.}
\begin{assumption}[L3] There exist a positive integer $\bar n$ and continuous functions $\alpha_{j}$ and $\beta_{j-1}$ from~$[-\pi,\pi]$ to $\mathbb R\,$, $j=1,\ldots,q$, independent of $n$, and such that
$$0< \beta_{j-1}(\theta) < \alpha_{j}(\theta)\le {\lambda_{nj}^X(\theta)}/{n}\le \beta_j(\theta)<\infty\quad\!\!\text{ $\theta$-a.e. in $[-\pi,\pi]$, all $j=1,\ldots, q$, and all $n>\bar n$. }$$
\end{assumption}
Under this assumption, the first $q$ common dynamic eigenvalues, irrespective of the frequency~$\theta$ (except possibly over a set of measure zero), are diverging linearly as~$n\to~\!\infty$. The following results then hold for the idiosyncratic dynamic eigenvalues
and those of the panel.
\begin{lemma} \label{lem:dyn_eval} Under Assumptions (L1) and (L3),
\begin{compactenum}[(i)]
\item there exists a constant $C^Z>0$ such that $\sup_{\theta\in[-\pi,\pi]} \lambda_{n1}^Z(\theta)\le C^Z$ for all $n\in\mathbb N$;
\item there exist a positive integer $\bar n$ and continuous functions $\alpha_{j}^Y$ and $\beta_{j-1}^Y$ from~$[-\pi,\pi]$ to $\mathbb R\,$,\linebreak $j=1,\ldots,q$, independent of $n$ and such that $0< \beta_{j-1}^Y(\theta) < \alpha_{j}^Y(\theta)\le {\lambda_{nj}^Y(\theta)}/{n}\le \beta_j^Y(\theta)\!<\infty$, $\theta$-a.e. in $[-\pi,\pi]$, all $j=1,\ldots, q$, and all $n>\bar n$;
\item there exists a constant $C^Y>0$ such that $\sup_{\theta\in[-\pi,\pi]} \lambda_{n,q+1}^Y(\theta)\le C^Y$ for all $n\in\mathbb N$.
\end{compactenum}
\end{lemma}
As a consequence of Lemma \ref{lem:dyn_eval}, identification of the model, i.e., consistently disentangling the unobserved common and idiosyncratic components, is possible, under the assumptions made in the limit, as $n\to\infty$, thanks to the behaviour of the dynamic eigenvalues.
Based on results by \citet{andersondeistler08} for singular vector processes with a rational spectrum, \citet{fornilippi11} and \citet{FHLZ15} prove that, for generic values of the coefficients of the filters $\mbf b_i(L)$ as defined in Assumption (L2), the space spanned by $u_{j,t-k}$ for~$j=1,\ldots, q$ and $k\ge 0$ is the same as the space spanned by any $(q+1)$-dimensional subvector of $\mbf X_{t}$ and its lags; moreover, those subvectors admit an autoregressive representation driven by the common shocks ${\bf u}_t$.
More precisely, any $(q+1)$-dimensional subvector $\mbf X^\ddag_t$ of ${\mbf X}_{nt}$ admits an autoregressive representation of the form
\begin{equation}\label{VAReq:lev}
\mbf A^\ddag(L)\mbf X^\ddag_t= \mbf H^\ddag\mbf u_t,
\end{equation}
where $\mbf A^\ddag(L)$ is a finite-order VAR operator such that $\mbf A^\ddag(0)=\mbf I_{q+1}$, $\mbf u_t$ is the vector of common shocks in \eqref{eq:common_lev_idio_vec}, and~$\mbf H^\ddag$ an appropriate~$(q+1)\times q$ matrix. On that representation, we make the following assumptions.
\begin{assumption}[L4]
Let $\mbf X^\ddag_t$ be an arbitrary $(q+1)$-dimensional subvector of ${\mbf X}_{nt}$: the autoregressive representation \eqref{VAReq:lev} is such that
\begin{compactenum}[(i)]
\item $\mbf A^\ddag(L)$ is uniquely defined;
\item the degree $S^\ddag$ of $\mbf A^\ddag(z)$ is uniformly bounded, that is, $S^\ddag\le S$ for some integer $S>0$ independent of $n$ and the choice of the subvector ${\mbf X}^\ddag_t$;
\item $\text{\rm det}[\mbf A^\ddag(z)]\neq 0$ for all $z\in\mathbb C$ such that $|z|\le 1$;
\item $\mbf H^\ddag$ is $(q+1)\times q$, with full rank $q$;
\item denoting by $\bm\Gamma_h^{X^\ddag}$ the lag-$h$ autocovariances of $\mbf X^\ddag:=\{ \mbf X^\ddag_t \vert t\in\mathbb{Z} \}$ and defining
\[
\bm{\mathcal C}^\ddag:=\left[
\begin{array}{cccc}
\bm\Gamma_0^{X^\ddag}& \bm\Gamma_1^{X^\ddag}& \cdots & \bm\Gamma_{S-1}^{X^\ddag}\\
\bm\Gamma_{-1}^{X^\ddag}&\bm\Gamma_0^{X^\ddag}& \cdots & \bm\Gamma_{S-2}^{X^\ddag}\\
\vdots&\vdots&\ddots&\vdots\\
\bm\Gamma_{-S+1}^{X^\ddag}& \bm\Gamma_{-S+2}^{X^\ddag}& \cdots & \bm\Gamma_{0}^{X^\ddag}\\
\end{array}
\right],
\]
$\text{\rm det}(\bm{\mathcal C}^\ddag) >d>0$, where $d$ is independent of the choice of the subvector~${\mbf X}^\ddag_t$.
\end{compactenum}
\end{assumption}
This assumption allows us to derive an alternative representation of the GDFM \eqref{eq:idio_level} which is particularly useful for estimation and for the construction, in Section \ref{sec:mod_vol} below, of a further GDFM for log-volatilities. Without loss of generality, let $n$ factorise into $n=m(q+1)$ for some positive integer~$m$, so that we can partition~${\mbf X}_n$ into~$m$ subprocesses, each of dimension $(q+1)$, of the form~${\mbf X}^{(k)}_t:=(X_{(k-1)(q+1),t}\ldots X_{k(q+1)-1,t})^\prime$, $k=1,\ldots,m$, with superscript $^{(k)}$ substituted for~$^\ddag$. Each ${\mbf X}^{(k)}$ satisfies \eqref{VAReq:lev} and Assumption~(L4). Defining the $n\times q$~ma\-trix~$\mbf H_n:=(\mbf H^{(k)'}\cdots \mbf H^{(m)'})'$, we thus have the VAR representation
\begin{equation}\label{comVAReq}
\mbf A_n(L) \mbf X_{nt} = \mbf H_n\mbf u_t,
\end{equation}
where $\mbf A_n(L)$ is $n\times n$ block-diagonal with diagonal blocks $\mbf A^{(1)}(L),\ldots, \mbf A^{(m)}(L)$. Moreover, in view of~\eqref{eq:common_lev_idio_vec}, we have $\left[\mbf A_n(L)\right]^{-1} \mbf H_n= \mbf B_n(L)$ (see Proposition 3 in \citealp{FHLZ17}).
Then, the following alternative and equivalent representation of the GDFM holds:
\begin{equation}\label{eq:gdfm_lev_static}
\mbf A_n(L)\left\{\mbf Y_{nt}-\mathrm E[\mbf Y_{nt}]\right\} = \mbf H_n\mbf u_t+\mbf A_n(L)\mbf Z_n.
\end{equation}
The advantage of this representation is that it is ``static'' in the sense that the common shocks $\mbf u$ now are loaded only contemporaneously and not via filters as in \eqref{eq:common_lev_idio_vec}.
To conclude with, note that the Yule-Walker equations
\begin{align}\label{eq:yw_pop}
\big(\mbf A_1^\ddag\cdots\mbf A_S^\ddag\big)=\big(\bm\Gamma_1^{X^\ddag}\cdots \bm\Gamma_S^{X^\ddag} \big)\big[\bm{\mathcal C}^\ddag\big]^{-1},
\end{align}
characterising the $S$ matrix coefficients of~$\mbf A^\ddag(L)$ in \eqref{VAReq}
are well defined in view of part~{\it (v)} of Assumption~(L4); the same conclusion holds, blockwise, for the $n$ -dimensional VAR \eqref{comVAReq}.
For ease of notation, define the filtered processes
\[
{\mbf Y}_n^*:=\mbf A_n(L)\left\{\mbf Y_n-\mathrm E[\mbf Y_n]\right\},\quad{\mbf X}^*_n:=\mbf A_n(L)\mbf X_n,\quad\text{and}\quad {\mbf Z}^*_n:=\mbf A_n(L)\mbf Z_n
\]
with traditional (static) covariance eigenvalues $\mu_{nj}^{Y^*}$, $\mu_{nj}^{X^*}$, and $\mu_{nj}^{Z^*}$, respectively.
Since~\eqref{eq:gdfm_lev_static} is a static factor model, it is natural to make the following assumption on the eigenvalues of the covariance of~${\mbf X}_n^*$ (see Assumption~4 in \citealp{FGLR09} or Assumption~6 in \citealp{FHLZ17}). Unless $q=1$, indeed, it does not even follow from Assumption~(L3) that $\mathrm E ({\mbf X}_n^*{\mbf X}_n^{*\prime})$ has rank $q$.
\begin{assumption}[L5] There exist a positive integer $\bar n$ and constants $a_{j}>b_{j-1}$, $j=1,\ldots, q$, independent of $n$ such that $0<a_j \le {\mu_{nj}^{X^*}}/{n}\le b_j<\infty$ for all $j=1,\ldots, q$ and all $n>\bar n$.
\end{assumption}
The following results then hold for the eigenvalues $\mu_{nj}^{Z^*}$ and $\mu_{nj}^{Y^*}$ of the covariance matrices of ${\mbf Z}^*_n$ and~${\mbf Y}^*_n$, respectively.
\begin{lemma}\label{lem:stat_eval}
Under Assumptions (L1), (L3), (L4), and (L5),
\begin{compactenum}[(i)]
\item there exists a constant $C^{Z^*}>0$ such that $\mu_{n1}^{Z^*}\le C^{Z^*}$ for all $n\in\mathbb{N}$;
\item there exist a positive integer $\bar n$ and constants $a_{j}^{Y^*}>b_{j-1}^{Y^*}$, $j=1,\ldots, q$, independent of $n$ such that \linebreak
$0<a_j^{Y^*}\le {\mu_{nj}^{Y^*}}/{n}\le b_j^{Y^*}<\infty$ for all $j=1,\ldots, q$ and all $n>\bar n$;
\item there exists a constant $C^{Y^*}>0$ such that $\mu_{n,q+1}^{Y^*}\le C^{Y^*}$ for all $n\in\mathbb{N}$.
\end{compactenum}
\end{lemma}
\subsection{Model and assumptions for volatilities}\label{sec:mod_vol}
We define the vector of common innovations (at time $t$) as the $n$-dimensional vector$$\mbf e_{nt}:=(e_{1t},\ldots,e_{nt})^\prime:=\mbf H_n\mbf u _t;$$
for $n>q$, the processes ${\mbf e}_n:=\{\mbf e_{nt}\vert t\in\mathbb{Z}\}$, $n\in\mathbb{N}$ clearly are singular. Then, letting $s_{it}:= e_{it}+ v_{it}$,
our log-volatility proxy is
\begin{equation}
h_{it} :=\log s_{it}^2= \log(e_{it}+ v_{it})^2,
\end{equation}
yielding the double-indexed stochastic process $\mbf h:= \{h_{it} \vert i\in\mathbb{N} , \ t\in\mathbb{Z}\}$, with $n$-dimensional sub-process\-es~$\mbf h_n:= \{h_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\}$. We call $\mbf h$ the ``log-volatilities'' process. Similar definitions are used in \citet{EM06} and our previous work (\citealp{barigozzihallin15a,barigozzihallin15b,barigozzihallin15c}, and \citealp{BHS18}). In order for such processes to be well defined we make the following assumption.
\begin{assumption}[V0] For all $i\in\mathbb N$ and $t\in\mathbb Z$, $\vert s_{it}\vert >0$ almost surely.
\end{assumption}
This assumption makes sure that no cancellation can happen between common and idiosyncratic innovations; it is required, since $e_i$ and $v_i$, although mutually orthogonal by Assumption~(L1.vi), need not be mutually independent (assuming, for instance, that $e_i$ and $v_i$ are absolutely continuous is not sufficient).
Assuming a GDFM with $Q$ factors for the log-volatilities, we obtain
\begin{align}
& h_{it}-\mathrm E[ h_{it}]= \chi_{it}+\xi_{it}\quad i\in\mathbb N ,\ t\in\mathbb Z\label{GDFMVol}\\
\text{with }\ & \chi_{it} = \sum_{j=1}^Q \sum_{k=0}^{\infty} f_{ijk} \varepsilon_{jt-k}=\mbf f_i'(L)\bm\varepsilon_t\quad\text{and}\quad \xi_{it}=\sum_{k=0}^{\infty} g_{ik} \nu_{it-k} = g_i(L) \nu_{it},\label{eq:idio_vol}
\end{align}
where $\mathrm E[h_{it}]$ is $h_{it}$'s expected value, $\chi_{it}$ and~$\xi_{it}$ are~$h_{it}$'s {\it common} and {\it idiosyncratic} components, and the pro\-cess\-es~$\bm\varepsilon:= \{\varepsilon_{jt} \vert j=1,\ldots,Q, \ t\in\mathbb{Z}\}$ and $\bm\nu_{n}:= \{\nu_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\}$, $n\in \mathbb N$ are mutually orthogonal (at all leads and lags) $Q$- and $n$-dimensional white noise, respectively. Note that a GDFM for log-volatilities implies a multiplicative GDFM representation \[
s_{it}^2=\exp(h_{it}) = \exp(\chi_{it}) \exp(\xi_{it}) \exp(\mathrm E[h_{it}]).
\]
for the volatilities themselves. Letting
$$\bm\chi_n:= \{\chi_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\}\quad\text{and}\quad\bm\xi_n:= \{\xi_{it} \vert i=1,\ldots,n, \ t\in\mathbb{Z}\},$$
equations \eqref{eq:idio_vol} in vector notation take the form
\begin{equation}\label{eq:common_vol_idio_vec}
\bm\chi_{nt}=\mbf F_n(L)\bm\varepsilon_t,\qquad \bm\xi_{nt}=\mbf G_n(L)\bm\nu_{nt}
\end{equation}
with $\mbf F_n(L):=(\mbf f_1(L)\ldots \mbf f_n(L))'$ and $\mbf G_n(L):=\text{\rm diag}(g_1(L)\ldots g_n(L))$.\medskip
The following assumptions then are the analogues, for log-volatilities and \eqref{GDFMVol}-\eqref{eq:idio_vol} , of Assumption~(L1).
\begin{assumption}[V1]$\,$
\begin{compactenum}[(i)]
\item The dimension $Q$ of $\bm\varepsilon_t$ does not depend on $n$; the process $\bm\varepsilon:=\{\bm\varepsilon_t \vert t\in\mathbb{Z}\}$ is second-order white noise, with mean $\mbf 0_Q$ and diagonal positive definite covariance~$\bm\Gamma^\varepsilon$;
\item writing ${\mbf f}_{ik}:=(f_{i1k}\ldots f_{iqk})^\prime$ for the $Q\times 1$ coefficient of $L^k$ in ${\mbf f}_i(L)$, there exists a constant $M_5>0$ such that $\sum_{k=0}^{\infty}\Vert \mbf f_{ik}\Vert\, \vert k\vert\leq M_5$ for all $i\in\mathbb{N}$;
\item the process $\{\bm\nu_{nt} \vert t\in\mathbb{Z}\}$ is second-order white noise, with mean $\mbf 0_n$ and positive definite covariance~$\bm\Gamma^\nu_n$; moreover, $\mathrm E[\nu_{it}|\nu_{is}]=0$ for all $i\in\mathbb N$ and and $t,s\in\mathbb Z$ such that $t>s$;
\item there exists a constant $C_{\nu}>0$ such that $\Vert \bm\Gamma_n^{\nu}\Vert_1 \leq C_{\nu}$ for all $n\in\mathbb{N}$;
\item there exists a constant $M_6>0$ such that $\sum_{k=0}^{\infty}\vert g_{ik}\vert\, \vert k \vert\le M_6$ for all $i\in\mathbb{N}$;
\item $\text{\rm Cov}(\varepsilon_{jt},\nu_{is})=0$ for all $i\in\mathbb N$, $j=1,\ldots, q$, and $t,s\in\mathbb Z$;
\item there exists a constant $M_7>0$ such that $\sum_{k_1,k_2,k_3\in\mathbb Z} \vert\mathrm E[\varepsilon_{j_1t-k_1}\varepsilon_{j_2t-k_2}\varepsilon_{j_3t-k_3}\varepsilon_{j_4t}]\vert\le M_7$ for all \linebreak $j_1,j_2,j_3,j_4=1,\ldots,Q$;
\item there exists a constant $M_8>0$ such that $\sum_{k_1,k_2,k_3\in\mathbb Z} \vert\mathrm E[\nu_{i_1t-k_1}\nu_{i_2t-k_2}\nu_{i_3t-k_3}\nu_{i_4t}]\vert\le M_8$ for all \linebreak$i_1,i_2,i_3,i_4\in\mathbb N$.
\end{compactenum}
\end{assumption}
The same comments made for Assumption (L1) apply here. Moreover, note that all moments of log-transforms of heavy-tailed variables exist and are finite, even for stable distributions (see e.g. Theorem 5.8.1 in \citealp{UZ11}). Pursuing with assumptions, the following one is the log-volatility counterpart of (L2).
\begin{assumption}[V2] $\,$
\begin{compactenum}[(i)]
\item $\mbf f_i(L)$ has rational entries $f_{ij}(L)=\tilde\theta_{ij}(L)\tilde\phi_{ij}^{-1}(L)$, where $\tilde\phi_{ij}(z)$ and $\tilde\theta_{ij}(z)$, for all $i\in\mathbb N$ and $j=1,\ldots, Q$, are finite-order polynomials;
\item there exists a constant $\underline \phi>1$ such that $\tilde\phi_{ij}(z)\neq 0$ for all $i\in\mathbb N$, all $j=1,\ldots, Q$, and all~$z\in\mathbb C$ such that~$|z|\le \underline \phi$;
\item the coefficients $\tilde\theta_{ijk}$ of $\tilde\theta_{ij}(L)$ are such that $|\tilde\theta_{ijk}|\le B^\chi$ for some constant $B^\chi>0$ and all $i\in\mathbb N$,\linebreak $j=1,\ldots, Q$, and~$k\in\mathbb N\cup \{0\}$;
\item $g_i(L)$ is of the form $p_i^{-1}(L)$ where $p_i(z)$, for all $i\in\mathbb N$, is a finite-order polynomial,~$p_{i}(0)=1$ and $p_{i}(z)\neq 0$ for all $z\in\mathbb C$ such that $|z|\le 1$.
\end{compactenum}
\end{assumption}
Assumptions (V2.iv) implies that we can rewrite \eqref{eq:idio_vol} also as
\begin{equation}
p_i(L) Z_{it} = \nu_{it}.\label{eq:var_vol_idio}
\end{equation}
As in the case of levels, this assumption could be relaxed to allow for an infinite autoregressive order.
Let $\bm\Sigma_n^h(\theta)$, $\bm\Sigma_n^\chi(\theta)$, and $\bm\Sigma_n^\xi(\theta)$, $\theta\in[-\pi,\pi]$ denote the $n\times n$ spectral density matrices of $\mbf h_n$, its common and its idiosyncratic components, with $j$-th largest eigenvalues $\lambda_{nj}^h(\theta)$, $\lambda_{nj}^\chi(\theta)$ and $\lambda_{nj}^\xi(\theta)$, respectively. As in~(L3), we assume the following.
\begin{assumption}[V3] There exist a positive integer $\bar n$ and continuous functions $\tilde\alpha_{j}(\theta)$ and $\tilde\beta_{j-1}(\theta)$ from~$[-\pi,\pi]$ to~$\mathbb R\,$, $j=1,\ldots,Q$, such that
$0<\tilde\beta_{j-1}(\theta)<\tilde\alpha_j(\theta)\le {\lambda_{nj}^\chi(\theta)}/{n}\le \tilde\beta_j(\theta)<\infty$, $\theta$-a.e. in $[-\pi,\pi]$, \linebreak all~$j=1,\ldots, Q$, and all $n>\bar n$.
\end{assumption}
Finally, the analogue (V4) of (L4) again is based on the representation results in \citet{FHLZ15}: \linebreak any~$(Q+1)$- dimensional subvector $\bm\chi^\ddag_t$ of $\bm\chi_{nt}$ admits an autoregressive representation of the form
\begin{equation}\label{VAReq}
\mbf M^\ddag(L)\bm\chi^\ddag_t= \mbf R^\ddag\bm\varepsilon_t,
\end{equation}
where $\mbf M^\ddag(L)$ is a finite-order VAR operator such that $\mbf M^\ddag(0)=\mbf I_{Q+1}$, $\bm\varepsilon_t$ is the vector of common shocks in~\eqref{eq:common_vol_idio_vec}, and~$\mbf R^\ddag$ an appropriate~$(Q+1)\times Q$ matrix. On that representation, we make the following assumptions:
\begin{assumption}[V4] $\,$
\begin{compactenum}[(i)]
\item $\mbf M^\ddag(L)$ is uniquely defined;
\item the degree $\tilde S^\ddag$ of $\mbf M^\ddag(z)$ is uniformly bounded, that is, $\tilde S^\ddag\le \tilde S$ for some integer $\tilde S>0$ independent of $n$ and the choice of the subvector ${\bm \chi}^\ddag_t$;
\item $\text{\rm det}[\mbf M^\ddag(z)]\neq 0$ for all $z\in\mathbb C$ such that $|z|\le 1$.
\item the $(Q+1)\times Q$ matrix $\mbf R^\ddag$ has full rank $Q$;
\item denoting by $\bm\Gamma_h^{\chi^\ddag}$ the lag-$h$ autocovariances of $\bm\chi^\ddag:=\{ \bm\chi^\ddag_t ,t\in\mathbb{Z} \}$ and defining $\bm{\mathcal V}^\ddag$ analogously to $\bm{\mathcal C}^\ddag$ in (L4),
$\text{\rm det}(\bm{\mathcal V}^\ddag) >\tilde d >0$, where $\tilde d$ is independent of the choice of the subvector~$\bm\chi^\ddag_t$.
\end{compactenum}
\end{assumption}
Now, Assumption (V4) implies $\left[\mbf M_n(L)\right]^{-1} \mbf R_n= \mbf F_n(L)$,
so that, assuming without loss of genera\-lity that~$n=\bar m(Q+1)$ (with $\bar m\neq m$ if $Q\ne q$) and defining a block-diagonal autoregressive operator $\mbf M_n(L)$ the way we defined $\mbf A_n(L)$ in the previous section, we can rewrite the GDFM for log-volatilities under the static form
\begin{equation}\label{eq:gdfm_vol_static}
\mbf M_n(L) \left\{\mbf h_{nt}-\mathrm E[\mbf h_{nt}]\right\} = \mbf R_n\bm\varepsilon_t+\mbf M_n(L)\bm\xi_{nt}.
\end{equation}
After defining, with obvious notation, the filtered processes $\vspace{1mm}{\mbf h}^*_n:=\mbf M_n(L)\left[\mbf h_n-\mathrm E[\mbf h_n]\right]$, ${\bm\chi}^*_n:=\mbf M_n(L)\bm\chi_n$, and~${\bm\xi}^*_n:=\mbf M_n(L)\bm\xi_n$, with (static) spectral eigenvalues $\mu_{nj}^{h^*}$, $\mu_{nj}^{\chi^*}$, and $\mu_{nj}^{\xi^*}$, we conclude with the analogues of~(L5) and Lemmas~\ref{lem:dyn_eval} and~\ref{lem:stat_eval} for the log-volatility panels.
\begin{assumption}[V5] There exist a positive integer $\bar n$ and constants $\tilde a_{j}>\tilde b_{j-1}>0$, $j=1,\ldots, Q$, independent of~$n$ such that $0< \tilde a_j \le {\mu_{nj}^{\chi^*}}/{n}\le \tilde b_j<\infty$ for all $j=1,\ldots, Q$ and all $n>\bar n$.
\end{assumption}
We then have the following.
\begin{lemma}\label{lem:eval_vol}
Under Assumptions (V0), (V1), (V3), (V4), and (V5),
\begin{compactenum}[(i)]
\item there exists a constant $C^\xi>0$ such that $\sup_{\theta\in[-\pi,\pi]} \lambda_{n1}^\xi(\theta)\le C^\xi$ for all $n\in\mathbb N$;
\item there exist a positive integer $\bar n$ and continuous functions $\alpha_{j}^h(\theta)$ and $\beta_{j-1}^h(\theta)$ from~$[-\pi,\pi]$ to $\mathbb R\,$, \linebreak $j=1,\ldots, Q$,~independent of~$n$ and such that~$0<\beta_{j-1}^h(\theta) <\alpha_j^h(\theta)\le {\lambda_{nj}^h(\theta)}/{n}\le \beta_j^h(\theta)\!<\infty$, $\theta$-a.e. in~$[-\pi,\pi]$, all $j=1,\ldots, Q$, and all $n>\bar n$;
\item there exists a constant $C^h>0$ such that $\sup_{\theta\in[-\pi,\pi]} \lambda_{n,Q+1}^h(\theta)\le C^h$ for all~$n\in\mathbb N$;
\item there exists a constant $C^{\xi^*}>0$ such that $\mu_{n1}^{\xi^*}\le C^{\xi^*}$ for all $n\in\mathbb{N}$;
\item there exist a positive integer $\bar n$ and constants $a_{j}^{h^*}>b_{j-1}^{h^*}$, $j=1,\ldots, Q$, independent of $n$ such \linebreak that~$0<a_j^{h^*}\le {\mu_{nj}^{h^*}}/{n}\le b_j^{h^*}<\infty$, for all~$j=1,\ldots, Q$ and all $n>\bar n$;
\item there exists a constant $C^{h^*}>0$ such that $\mu_{n,Q+1}^{h^*}\le C^{h^*}$ for all $n\in\mathbb{N}$.
\end{compactenum}
\end{lemma}
\setcounter{equation}{0}
\section{Estimation, consistency, and rates}\label{sec:est}
Hereafter, the terminology ``estimation'', ``estimator'', etc.\ is used, in an orthodox way, for data-driven quantities attempting at evaluating parameters (covariances, spectra, loadings, \ldots) but also, with a slight abuse, for data-driven quantities attempting at reconstructing unobserved variables (such as common factors, common and idiosyncratic components, \ldots ). All those ``estimators'', which are~${\bf Y}_{n,T}$-measurable random variables (hence depend both on $n$ and $T$) are carrying ``hats''.
\subsection{Summary of estimation}\label{sec:est_summary}
Estimation proceeds in two parts. The first part deals with the observed $n\times T$ panel ${\bf Y}_{n,T}$ of levels,
and follows along similar lines as in \citet{FHLZ17}, yielding estimated log-volatility proxies; the second part consists in repeating the same estimation steps, now based on those estimated log-volatility quantities.
Global consistency of the procedure is discussed in the next section, along with further necessary conditions.
To start with, we assume that $q$ and $Q$ are known---an assumption we are relaxing later on.
For simplicity of notation, we also assume $\mbf Y_n$ and $\mbf h_n$ to be centred, i.e., to have zero mean; in practice, sample means are to be subtracted in order to obtain centred variables---which has no impact on consistency nor consistency rates.
Here is a detailed list of the steps required for estimation. Further comments on the choice of the quantities needed for estimation and a schematic description of the procedure are given at the end of this section (see also Algorithms 1 and 2).
\begin{enumerate}
\item [(\textit {L.i})] To start with, compute the lag-window estimator
\[
\widehat{\bm\Sigma}_{n}^Y(\theta_h):=\frac 1{2\pi} \sum_{k=-T+1}^{T-1}\mathrm K\left(\frac{k}{B_T}\right)e^{-ik\theta_h}\widehat{\bm\Gamma}_{nk}^Y,\quad \theta_h=\frac{\pi h}{B_T}, \quad \vert h\vert \le B_T,
\]
of the spectral density matrix of returns, where $\widehat{\bm\Gamma}_{nk}^Y:=T^{-1}\sum_{t=|k|+1}^T {\mbf Y}_{nt}{\mbf Y}_{nt-|k|}'$ is the usual lag-$k$ sample autocovariance matrix of levels and $\mathrm K$ is a suitable kernel with bandwidth $B_T$.
We here adopt the common choice of a Bartlett kernel
\[
\mathrm K\left(x\right)=\left\{\begin{array}{cl}
1-|x|& \mbox{if } |x|\le 1\\
0& \mbox{otherwise},
\end{array}
\right.
\]
but other classical kernels are also possible.
\item [(\textit{L.ii})] Collect the $q$ normalised column eigenvectors associated with $\widehat{\bm\Sigma}^Y_{n}(\theta_h)$'s $q$ largest
eigenvalues into the $n\times q$ matrix $\widehat{\bf P}^Y_{n} (\theta_h)$, and collect the corresponding eigenvalues into the~$q\times q$ diagonal matrix $\widehat{\bm \Lambda}^Y_{n} (\theta_h)$. Take
\[
\widehat{\bm\Sigma}_{n}^X(\theta_h):=\widehat{\bf P}^Y_{n} (\theta_h) \widehat{\bm \Lambda}^Y_{n} (\theta_h) \widehat{\bf P}^{Y\dag}_{n} (\theta_h),
\]
as an estimate of the spectral density matrix of the level-common component process~${\bf X}_n$.
\item [(\textit{L.iii})] By inverse Fourier transform of $\widehat{\bm\Sigma}_{n}^X(\theta_h)$, estimate the autocovariance matrices of $\mbf X_n$:
\[
\widehat{\bm\Gamma}^X_{nk}:= \frac{\pi}{B_T}\sum_{h=-B_T}^{B_T} e^{ik\theta_h} \widehat{\bm\Sigma}_{n}^X(\theta_h), \qquad k\in\mathbb Z.
\]
\item [(\textit{L.iv})] Assuming, for simplicity,
\footnote{{In practice, the last $n-\lfloor n/(q+1) \rfloor (q+1)$ cross-sectional items can be added to the last block in the analysis which will then have size larger than $(q+1)$. Since the arguments in \citet{FHLZ17} used in the next section apply to any partition of blocks of size $(q+1)$ or larger, nothing changes in what follows. }}
that $n=m(q+1)$, consider the $m$ diagonal $(q+1)\times (q+1)$ blocks of the
$\widehat{\bm\Gamma}^X_{nk}$'s. For each block, estimate, via Yule-Walker methods, the coefficients of a~$(q+1)$-dimensional VAR model (order determined via AIC or BIC). In other words, compute the sample analogue of \eqref{eq:yw_pop}. This yields, for the $\ell$-th diagonal block, an estimator~$\widehat{\mbf A}^{(\ell)}(L)$ of the autoregressive filter $\mbf A^{(\ell)}(L)$ appearing in Assumption~(L4), hence an estimator $\widehat{\mbf A}_n(L)$ of the VAR filter $\mbf A_n(L)$. The resulting estimated filtered process and its estimated covariance matrix are~$
\widehat{{\mbf Y}}^*_{nt}:=\widehat{\mbf A}_n(L)\mbf Y_{nt}$ and
$\widehat{\bm\Gamma}^{\widehat{{Y}}^*}_n:=T^{-1}\sum_{t=1}^T\widehat{{\mbf Y}}^*_{nt}\widehat{{\mbf Y}}^{*'}_{nt}$,
respectively.
\item [(\textit{L.v})] Collect the $q$ normalised (column) eigenvectors corresponding to $\widehat{\bm\Gamma}^{\widehat{{Y}}^*}_n$'s $q$ largest
eigenvalues into the~$n\times q$ matrix $\widehat{\bf Q}^{\widehat{{Y}}^*}_{n}$. Projecting $\widehat{{\mbf Y}}^*_{nt}$ onto the space spanned by the columns of $\widehat{\bf Q}^{\widehat{{Y}}^*}_{n}$ provides an estimate $\widehat{\mbf e}_n$ of the innovation process~$\mbf e_n$. Taking into account the set of identifying restrictions described in Assumption~(I) below, we obtain the estimators
\[
\widehat{\mbf H}_n:= \sqrt n \widehat{\bf Q}^{\widehat{{Y}}^*}_{n},\qquad \widehat{\mbf u}_t:=\frac{1}{n}{\widehat{\mbf H}_n' \widehat{{\mbf Y}}^*_{nt}},\quad\text{and}\quad \widehat{\mbf e}_{nt}:=\widehat{\mbf H}_n\widehat{\mbf u}_t=\widehat{\bf Q}^{\widehat{{Y}}^*}_{n}\widehat{\bf Q}^{\widehat{{Y}}^{*'}}_{n}\widehat{{\mbf Y}}^*_{nt}.
\]
Our estimator of the dynamic loadings then is
$
\widehat{\mbf B}_n(L):= \widehat{\mbf A}_n^{-1}(L)\widehat{\mbf H}_n$,
where we truncate the filter $\widehat{\mbf A}_n^{-1}(L)$ at some finite lag $\bar k_1$. From this we obtain an estimator
$
\widehat{\mbf X}_{nt} := \widehat{\mbf B}_n(L)\widehat{\mbf u}_t$
of the common component.
\item [(\textit{L.vi})] The resulting estimator of the idiosyncratic component is $\widehat{\mbf Z}_{nt} := \mbf Y_{nt}-\widehat{\mbf X}_{nt}$. Fitting a univariate AR model (order determined via AIC or BIC), either by least squares or via Yule-Walker methods, to each of the $n$ components of~$\widehat{\mbf Z}_{nt}$ yields estimators $\widehat{\mbf v}_n$ of the residuals and $\widehat{\mbf C}_n(L)$ of the diagonal matrix of coefficients from which we also obtain $\widehat{\mbf D}_n(L)\!:=\!\widehat{\mbf C}_n^{-1}(L)$ with $\widehat{\mbf C}_n^{-1}(L)$ truncated at some finite lag~$\bar k_2$.
\item [(\textit{R})] For all $i=1,\ldots, n$ and $t=1,\ldots,T$, let $\widehat s_{it}:= \widehat e_{it}+\widehat v_{it}$ and define the estimated log-volatility proxies as capped values of $\log (\widehat s_{it}^{\ 2})$:
\[
\widehat h_{it}:=\log (\widehat s_{it}^{\ 2})\ \mathbb I(\vert \widehat s_{it}\vert\ge \kappa_T) +\log (\kappa_T^2)\ \mathbb I(\vert \widehat s_{it}\vert< \kappa_T),
\]
where $\kappa_T> 0$ is a sequence of constants to be chosen in order to make our proxy robust to the log-transform. Note that consistency of our estimation procedure requires an adaptive choice of~$\kappa_T$, depending on the sample size as explained in Assumption (R) below.
In particular, $\kappa_T$ must be strictly positive for consistency to hold.
\item [(\textit{V.i})] Denote by $\widehat{\mbf h}_{nt}:=\big(\widehat h_{1t} \ldots \widehat h_{nt}\big)^\prime$, $t=1,\ldots,T$ the $n$-dimensional vector of log-volatility proxies and compute the lag-window estimator
\[
\widehat{\bm\Sigma}_{n}^{\widehat{h}}(\theta_\ell):=\frac 1{2\pi} \sum_{k=-T+1}^{T-1}\mathrm K\left(\frac{k}{M_T}\right)e^{-ik\theta_\ell}\widehat{\bm\Gamma}_{nk}^{\widehat h},\quad \theta_\ell=\frac{\pi \ell}{M_T}, \quad \vert \ell\vert \le M_T,
\]
of its spectral density matrix, where $\widehat{\bm\Gamma}_{nk}^{\widehat h}:=T^{-1}\sum_{t= |k|+1}^T {\widehat{\mbf h}}_{nt}{\widehat{\mbf h}}_{n,t-|k|}'$ is the lag-$k$ sample autocovariance matrix of estimated log-volatilities. Here again we adopt the Bartlett kernel, with bandwidth $M_T$, which could be different from $B_T$ in step (\textit{L.i}).
\item [(\textit{V.ii})-(\textit{V.vi})] Repeat steps (\textit{L.ii})-(\textit{L.vi}) for $\widehat{\mbf h}_n$. In particular, steps (\textit{V.ii})-(\textit{V.v}) yield the estimators~$\widehat{\mbf M}_n(L)$ and $\widehat{\mbf R}_n$, from which we compute
\begin{align}
\widehat{\mbf h}_{nt}^* := \widehat{\mbf M}_n(L)\widehat{\mbf h}_{nt},\qquad
\widehat{\bm\varepsilon}_t:=\frac{1}{n}{\widehat{\mbf R}_n'\widehat{\mbf h}_{nt}^*},\quad \text{and}\quad
\widehat{\mbf F}_n(L):= \widehat{\mbf M}_n^{-1}(L)\widehat{\mbf R}_n,\nonumber
\end{align}
while from (\textit{V.vi}) we obtain $\widehat{\bm\nu}_{n}$ and $\widehat{\mbf P}_n(L)$, hence $\widehat{\mbf G}_n(L):=\widehat{\mbf P}_n^{-1}(L)$. As before, $\widehat{\mbf M}_n^{-1}(L)$ and~$\widehat{\mbf P}_n^{-1}(L)$ are truncated at finite lags $\bar k_1^{*}$ and $\bar k_2^{*}$.
\end{enumerate}
\vskip .2cm
\begin{algorithm}\label{tab:alg1}
\DontPrintSemicolon
\footnotesize{
\KwInput{data in levels $\mathbf Y$ of dimension $n\times T$, number of factors $q$, bandwidth for estimating spectral density $B_T$, number of lags for impulse responses $\bar k_1$ and $\bar k_2$, number of permutations for estimating the common component $nrep$}
\KwOutput{common component $\widehat{\mathbf X}$, idiosyncratic component $\widehat{\mathbf Z}$, common shocks $\widehat {\mathbf u}$ and $\widehat{\mathbf e}$, common impulse responses $\widehat{\mathbf B}(L)$, idiosyncratic shocks $\widehat{\mathbf v}$, idiosyncratic impulse responses $\widehat{\mathbf D}(L)$}
\medskip
Compute autocovariance matrices of data $\widehat{\bm \Gamma}_{k}^Y$ for $|k|\le B_T$\smallskip
Compute the lag-window estimator of the spectral density matrix of data $\widehat{\bm\Sigma}^Y(\theta_h)$ for $\theta_h=\pi h/B_T$ and $|h|\le B_T$, using $\widehat{\bm \Gamma}_{k}^Y$ and the Bartlett kernel \smallskip
\For{$h\leftarrow -B_T$ \KwTo $B_T$}
{Compute the $q$ largest eigenvalues $\widehat \lambda_{1}^Y(\theta_h),\ldots,\widehat \lambda_{q}^Y(\theta_h)$ of $\widehat{\bm\Sigma}^Y(\theta_h)$ and collect the corresponding eigenvectors into the columns of
$\widehat{\mathbf P}^Y(\theta_h)$. Let $\widehat{\bm\Sigma}^X(\theta_h)=\widehat{\mathbf P}^Y(\theta_h)\mbox{diag}(\widehat \lambda_{1}^Y(\theta_h),\ldots,\widehat \lambda_{q}^Y(\theta_h))\widehat{\mathbf P}^{Y\dag}(\theta_h)$}\smallskip
Compute the autocovariance matrices of the common component $\widehat{\bm \Gamma}_{nk}^X$ for $|k|\le B_T$ by inverse Fourier transform of $\widehat{\bm\Sigma}^X(\theta_h)$\smallskip
\For{$\mathcal P\leftarrow 1$ \KwTo $nrep$}
{Choose a random partition $\mathcal P(1),\ldots, \mathcal P(m(q+1))$ of the $n$ series into $m=\lfloor n/(q + 1)\rfloor(q + 1)$ blocks
such that the first $q$ series are always included
and let $\mathbf Y_{\mathcal P}=(Y_{\mathcal P(1)},\ldots ,Y_{\mathcal P(m(q+1))})^\prime$\,\smallskip
\If({Add the last $n-m$ series to the last block}\,){$m(q+1)<n$}\smallskip
\For{$\ell \leftarrow 1$ \KwTo $m$}{Obtain the coefficients $\widehat{\mathbf A}^{(\ell)}_{\mathcal P}(L)$ fitting a VAR($p_1^{(\ell)}$) on $\mathbf Y^{(\ell)} := (Y_{\mathcal P((\ell-1)(q+1))}\ldots Y_{\mathcal P(\ell(q+1)-1)})^\prime$ via Yule Walker equations using $\widehat{\bm \Gamma}_{k}^X$ for $k=0,\ldots, \ell$, with $p_1^{(\ell)}\le B_T$ and determined via BIC}\,\smallskip
Let $\widehat{\mathbf A}_{\mathcal P}(L)=\mbox{diag}(\widehat{\mathbf A}^{(1)}_{\mathcal P}(L),\ldots, \widehat{\mathbf A}^{(m)}_{\mathcal P}(L))$ and let $\widehat{\mathbf Y}_{t,\mathcal P}^{*}=\widehat{\mathbf A}_{\mathcal P}(L)\mathbf Y_{t,\mathcal P}$ for $t=1,\ldots, T$\,\smallskip
Compute $\widehat{\mathbf H}_{\mathcal P}$ as $\sqrt n$ times the $q$ leading eigenvectors of the sample covariance matrix of $\widehat{\mathbf Y}_{\mathcal P}^{*}$ \,\smallskip
Compute $\widetilde{\mathbf B}_{\mathcal P}(L)=\widehat{\mathbf A}^{-1}_{\mathcal P}(L)\widehat{\mathbf H}_{\mathcal P}$ truncating at lag $\bar k_1$\,\smallskip
Compute $\widehat{\mathbf B}_{\mathcal P}(L)=\widetilde{\mathbf B}_{\mathcal P}(L)\bm{\mathcal R}_{\mathcal P}$ with $\bm{\mathcal R}_{\mathcal P}$ is $q\times q$ orthogonal and such that the $q\times q$ block of $\widehat{\mathbf B}_{\mathcal P}(0)$ obtained by isolating the rows corresponding to the first $q$ series in $\mathbf Y$ is lower triangular\,\smallskip
Compute $\widehat{\mathbf u}_{t,\mathcal P}=n^{-1}\bm{\mathcal R}_{\mathcal P}^\prime\widehat{\mathbf H}_{\mathcal P}^\prime\widehat{\mathbf Y}_{t,\mathcal P}^{*}$ for $t=1,\ldots, T$\,\smallskip
}
Compute the common shocks as $\widehat{\mathbf u}_t=(nrep)^{-1}\sum_{\mathcal P=1}^{nrep}\widehat{\mathbf u}_{t,\mathcal P}$ for $t=1,\ldots, T$\,\smallskip
Compute $\widehat{\mathbf e}_t=(nrep)^{-1}\sum_{\mathcal P=1}^{nrep}\widehat{\mathbf H}_{\mathcal P}\bm{\mathcal R}_{\mathcal P}\widehat{\mathbf u}_{t,\mathcal P}$ \,\smallskip
Compute the impulse response functions $\widehat{\mathbf B}(L)=(nrep)^{-1}\sum_{\mathcal P=1}^{nrep}\widehat{\mathbf B}_{\mathcal P}(L)$\,\smallskip
Compute the common component as $\widehat{\mathbf X}_{t}=\widehat{\mathbf B}(L)\widehat{\mathbf u}_{t}$ for $t=1,\ldots, T$\,\smallskip
Compute the idiosyncratic component as $\widehat{\mathbf Z}={\mathbf Y}-\widehat{\mathbf X}$ such that $\widehat{\mathbf Z}=(\widehat Z_1\ldots \widehat Z_n)^\prime$\,\smallskip
\For{$i\leftarrow 1$ \KwTo $n$}{Obtain the coefficients $\widehat{c}_i(L)$ fitting a VAR($s_{1i}$) on $\widehat{Z}_i$ via least squares, with $s_{1i}$ determined via BIC\,\smallskip
Let $\widehat{v}_{it}=\widehat{c}_i(L)\widehat{Z}_{it}$ for $t=1,\ldots, T$}\,\smallskip
Compute the impulse response functions as $\widehat{\mathbf D}(L)=\mbox{diag}(\widehat c_1^{\,-1}(L),\ldots, \widehat c_n^{\, -1}(L))$ truncating at lag $\bar k_2$\,\smallskip
Let the idiosyncratic shocks be $\widehat{\mathbf v}_t=(\widehat{v}_{1t}\ldots \widehat{v}_{nt})^\prime$ for $t=1,\ldots, T$
\smallskip
}
\caption{\small Estimation of dynamic factor model for levels}
\end{algorithm}
\begin{algorithm}\label{tab:alg2}
\DontPrintSemicolon
\footnotesize{
\KwInput{from Algorithm 1: common and idiosyncratic shocks $\widehat{\mathbf e}$ and $\widehat{\mbf v}$ both of dimension $n\times T$\\
number of factors $Q$, capping constant $\kappa_T$, bandwidth for estimating spectral density $M_T$, number of lags for impulse responses~$\bar k_1^*$ and $\bar k_2^*$, number of permutations for estimating the common component $nrep$}
\KwOutput{common component $\widehat{\bm\chi}$, idiosyncratic component $\widehat{\bm \xi}$, common shocks $\widehat {\bm\varepsilon}$ and $\widehat{\bm\eta}$, common impulse responses $\widehat{\mathbf F}(L)$, idiosyncratic shocks $\widehat{\bm\nu}$, idiosyncratic impulse responses $\widehat{\mathbf G}(L)$}
\medskip
\For{$i\leftarrow 1$ \KwTo $n$}{
\For{$t\leftarrow 1$ \KwTo $T$}{Compute log-volatility proxy $\widehat h_{it}$\,\smallskip
\If({$\widehat h_{it}=\log(\widehat e_{it}+\widehat v_{it})^2$}\,){$|\widehat e_{it}+\widehat v_{it}|\ge \kappa_T$\,}\smallskip
\Else($\widehat h_{it}=\kappa_T$)\,
}}
Compute autocovariance matrices of log-volatility $\widehat{\bm \Gamma}_{k}^{\widehat h}$ for $|k|\le M_T$\smallskip
Compute the lag-window estimator of the spectral density matrix of log-volatility $\widehat{\bm\Sigma}^{\widehat h}(\theta_h)$ for $\theta_h=\pi h/M_T$ and $|h|\le M_T$, using $\widehat{\bm \Gamma}_{k}^{\widehat h}$ and the Bartlett kernel \smallskip
\For{$h\leftarrow -M_T$ \KwTo $M_T$}
{Compute the $Q$ largest eigenvalues $\widehat \lambda_{1}^{\widehat h}(\theta_h),\ldots,\widehat \lambda_{Q}^{\widehat h}(\theta_h)$ of $\widehat{\bm\Sigma}^{\widehat h}(\theta_h)$ and collect the corresponding eigenvectors into the columns of
$\widehat{\mathbf P}^{\widehat h}(\theta_h)$. Let $\widehat{\bm\Sigma}^{\widehat \chi}(\theta_h)=\widehat{\mathbf P}^{\widehat h}(\theta_h)\mbox{diag}(\widehat \lambda_{1}^{\widehat h}(\theta_h),\ldots,\widehat \lambda_{Q}^{\widehat h}(\theta_h))\widehat{\mathbf P}^{\widehat h \dag}(\theta_h)$}\smallskip
Compute the autocovariance matrices of the common component $\widehat{\bm \Gamma}_{nk}^{\widehat \chi}$ for $|k|\le M_T$ by inverse Fourier transform of $\widehat{\bm\Sigma}^{\widehat \chi}(\theta_h)$\smallskip
\For{$\mathcal P\leftarrow 1$ \KwTo $nrep$}
{Choose a random partition $\mathcal P(1),\ldots, \mathcal P(m(Q+1))$ of the $n$ series into $m=\lfloor n/(Q + 1)\rfloor(Q + 1)$ blocks
such that the first $Q$ series are always included
and let $\widehat {\mbf h}_{\mathcal P}=({\widehat h}_{\mathcal P(1)},\ldots ,{\widehat h}_{\mathcal P(m(Q+1))})^\prime$\,\smallskip
\If({Add the last $n-m$ series to the last block}\,){$m(Q+1)<n$}\smallskip
\For{$\ell \leftarrow 1$ \KwTo $m$}{Obtain the coefficients $\widehat{\mathbf M}^{(\ell)}_{\mathcal P}(L)$ fitting a VAR($p_2^{(\ell)}$) on $\widehat{\mathbf h}^{(\ell)} := ({\widehat h}_{\mathcal P((\ell-1)(Q+1))}\ldots {\widehat h}_{\mathcal P(\ell(Q+1)-1)})^\prime$ via Yule Walker equations using $\widehat{\bm \Gamma}_{k}^{\widehat \chi}$ for $k=0,\ldots, \ell$, with $p_2^{(\ell)}\le M_T$ and determined via BIC}\,\smallskip
Let $\widehat{\mathbf M}_{\mathcal P}(L)=\mbox{diag}(\widehat{\mathbf M}^{(1)}_{\mathcal P}(L),\ldots, \widehat{\mathbf M}^{(m)}_{\mathcal P}(L))$ and let $\widehat{\mathbf h}_{t,\mathcal P}^{*}=\widehat{\mathbf M}_{\mathcal P}(L)\widehat{\mathbf h}_{t,\mathcal P}$ for $t=1,\ldots, T$\,\smallskip
Compute $\widehat{\mathbf R}_{\mathcal P}$ as $\sqrt n$ times the $q$ leading eigenvectors of the sample covariance matrix of $\widehat{\mathbf h}_{\mathcal P}^{*}$ \,\smallskip
Compute $\widetilde{\mathbf F}_{\mathcal P}(L)=\widehat{\mathbf M}^{-1}_{\mathcal P}(L)\widehat{\mathbf R}_{\mathcal P}$ truncating at lag $\bar k_1^*$\,\smallskip
Compute $\widehat{\mathbf F}_{\mathcal P}(L)=\widetilde{\mathbf F}_{\mathcal P}(L)\bm{\mathcal R}_{\mathcal P}$ with $\bm{\mathcal R}_{\mathcal P}$ is $Q\times Q$ orthogonal and such that the $Q\times Q$ block of $\widehat{\mathbf M}_{\mathcal P}(0)$ obtained by isolating the rows corresponding to the first $Q$ series in $\widehat{\mathbf h}$ is lower triangular\,\smallskip
Compute $\widehat{\bm \varepsilon}_{t,\mathcal P}=n^{-1}\bm{\mathcal R}_{\mathcal P}^\prime\widehat{\mathbf R}_{\mathcal P}^\prime\widehat{\mathbf h}_{t,\mathcal P}^{*}$ for $t=1,\ldots, T$\,\smallskip
}
Compute the common shocks as $\widehat{\bm\varepsilon}_t=(nrep)^{-1}\sum_{\mathcal P=1}^{nrep}\widehat{\bm\varepsilon}_{t,\mathcal P}$ for $t=1,\ldots, T$\,\smallskip
Compute $\widehat{\bm\eta}_t=(nrep)^{-1}\sum_{\mathcal P=1}^{nrep}\widehat{\mathbf R}_{\mathcal P}\bm{\mathcal R}_{\mathcal P}\widehat{\bm\varepsilon}_{t,\mathcal P}$ \,\smallskip
Compute the impulse response functions $\widehat{\mathbf F}(L)=(nrep)^{-1}\sum_{\mathcal P=1}^{nrep}\widehat{\mathbf F}_{\mathcal P}(L)$\,\smallskip
Compute the common component as $\widehat{\bm\chi}_{t}=\widehat{\mathbf F}(L)\widehat{\bm\varepsilon}_{t}$ for $t=1,\ldots, T$\,\smallskip
Compute the idiosyncratic component as $\widehat{\bm\xi}=\widehat{\mathbf h}-\widehat{\bm\chi}$ such that $\widehat{\bm\xi}=(\widehat {\xi}_1\ldots \widehat {\xi}_n)^\prime$\,\smallskip
\For{$i\leftarrow 1$ \KwTo $n$}{Obtain the coefficients $\widehat{p}_i(L)$ fitting a VAR($s_{2i}$) on $\widehat{\xi}_i$ via least squares, with $s_{2i}$ determined via BIC\,\smallskip
Let $\widehat{\nu}_{it}=\widehat{p}_i(L)\widehat{\xi}_{it}$ for $t=1,\ldots, T$}\,\smallskip
Compute the impulse response functions as $\widehat{\mathbf G}(L)=\mbox{diag}(\widehat p_1^{\,-1}(L),\ldots, \widehat p_n^{\, -1}(L))$ truncating at lag $\bar k_2^*$\,\smallskip
Let the idiosyncratic shocks be $\widehat{\bm\nu}_t=(\widehat{\nu}_{1t}\ldots \widehat{\nu}_{nt})^\prime$ for $t=1,\ldots, T$
\smallskip
}
\caption{\small Estimation of dynamic factor model for log-volatilities}
\end{algorithm}
An important remark needs to be made here. The cross-sectional ordering of the panel has an impact on the selection of the diagonal blocks in steps (\textit{L.iv}) and (\textit{V.iv}). Each cross-sectional permutation of the panel, thus, would lead to distinct estimators---all sharing the same asymptotic properties. A Rao-Blackwell argument (see \citealp{FHLZ17} for details) suggests aggregating these estimators into a unique one by simple averaging (after obvious reordering of the cross-section) of the resulting estimated shocks.~Although averaging over all~$n!$ permutations is clearly unfeasible, as stressed by \citet{FHLZ17} and verified empirically also in \citet{FGLS18}, a few of them are enough, in practice, to deliver stable averages (which therefore are matching the infeasible average over all~$n!$ permutations).
Implementation of the above estimation steps is described in Algorithms 1 and 2. Those algorithms require setting bandwidths $B_T$ and $M_T$ for the estimation of the spectral densities, a capping constant $\kappa_T$, and the number of factors $q$ and $Q$. Concerning the bandwidths and the capping constant, we refer to Section \ref{sec:ap} for the required asymptotic properties (see Assumptions (K) and (R), respectively), while a numerical assessment of the impact of these quantities is provided in Section \ref{sec:sim} on simulated data (see also the results in Appendix \ref{app:sim}) and in Section~\ref{sec:emp} on real data. Overall, our numerical analysis shows that low levels of capping or even no capping at all are preferable, as they avoid inducing too much bias in the log-volatility distributions. As for the bandwidths, large values of $T$ are required to construct reliable estimates, since they allow setting $M_T$ large enough to capture the high persistence of log-volatility series. Our results are quite insensitive to the choice of $B_T$, due to the fact that financial returns typically are only weakly autocorrelated.
Finally, we can determine the numbers $q$ and $Q$ of common shocks by means of the information criteria proposed by \citet{hallinliska07} and applied on the panels $\mbf Y_n$ and $\widehat{\mbf h}_n$, respectively. The resulting data-driven estimators $\widehat{q}$ and $\widehat{Q}$ converge in probability to $q$ and $Q$, respectively. Since $q$ and $Q$ are integers, this means that, for any $\epsilon >0$, there exist
$n(\epsilon)$ and $T(\epsilon)$ such that, for all $n>n(\epsilon)$ and~$T>T(\epsilon)$, $\widehat{q}=q$ and $\widehat{Q}=Q$ with probability larger than $1-\epsilon$. Hence, in Section~3.2 below, we safely can assume that $q$ and~$Q$ are known.
\subsection{Consistency and rates}\label{sec:ap}
Consistency of the estimators of the GDFM model for levels is proved in \citet{FHLZ17}. Some differences exist, though, between their approach and ours.
First, \citet{FHLZ17} make slightly weaker assumptions on idiosyncratic serial dependence and, by exploiting results in \citet{WZ18} on spectral density estimation, they derive their consistency results under the constraint that~${B_T\log B_T}/T\to 0$ as $T\to\infty$. A more classical approach is adopted here, based on Assumptions~(L1) and~(V1), which as a consequence requires mildly stronger constraints on the range of admissible values for the bandwidths $B_T$ and $M_T$. Specifically, we require the following.
\begin{assumption}[K] As $T\to\infty$, ${B_T}=o(\sqrt T)$ and ${M_T}=o(\sqrt T)$.
\end{assumption}
Note that for $T\simeq 1000$ as in our empirical study, the range of admissible bandwidths is still such that most of the serial dependence in the data is captured when estimating the spectral density (see Section \ref{sec:emp} for more details on the choice of the bandwidths).
Second, the results in \citet{FHLZ17} hold pointwise in $t$, which is not sufficient for our needs when it comes to prove consistency in the second part of the estimation procedure. Indeed, we need uniform (over all~$t\in\{1,\ldots,T\}$) consistency of the estimators of the common and idiosyncratic components. For this reason, we make additional assumptions on the distribution of common and idiosyncratic components.
\begin{assumption}[T] There exist constants $K_u>0$, $K_{\varepsilon}>0$, $K_Z>0$, and $K_\xi>0$, such that, for any~$t=1,\ldots ,T$,
\begin{compactenum}[(i)]
\item $\max_{j=1,\ldots, q}\Vert u_{jt}\Vert_{\psi_1}\le K_u$;
\item $\max_{j=1,\ldots, Q}\Vert \varepsilon_{jt}\Vert_{\psi_1}\le K_\varepsilon$;
\item $\sup_{\bm w_n: \Vert\bm w_n\Vert = 1}\Vert \bm w_n'\mbf Z_{nt} \Vert_{\psi_1}\le K_Z$, for all $n\in\mathbb N$;
\item $\sup_{\bm w_n: \Vert\bm w_n\Vert = 1}\Vert \bm w_n'\bm \xi_{nt} \Vert_{\psi_1}\le K_\xi$, for all $n\in\mathbb N$.
\end{compactenum}
\end{assumption}
This assumption is equivalent to an assumption of sub-exponential tails of the common factors and the normed linear combinations of idiosyncratic components. Specifically, it can be shown that (T{\it i}) is equivalent to requiring for any $j=1,\ldots,q$, that $\mathrm{P}(|u_{jt}|>\epsilon)\le K_u^* \exp\left(- {\epsilon}/K_{u}^{**}\right)$ for any $\epsilon>0$ and some finite $K_u^*,K_u^{**}>0$ (see also \citealp{vershynin12}, and Appendix \ref{app:prop1} for details). The same holds also for (T{\it ii}), (T{\it iii}), and (T{\it iv}). See Remark 1 at the end of this section for a discussion of the implications and possible relaxations of this assumption.
Two remarks on (T{\it iii}) and (T{\it iv}) are in order here (see Sections 5.2.4 and 5.2.5 in \citealp{vershynin12} for details). First, note that by letting $\bm w_n =(0\ldots w_i \ldots 0)'$, with $w_i=1$ for a given $i$, those assumptions imply that each idiosyncratic component has marginal sub-exponential distribution. Second, an implication of Lemmas \ref{lem:stat_eval} and \ref{lem:eval_vol} is that vectors of the form~$\bm w_n'\mbf Z_n$ and $\bm w_n'\bm\xi_n$ have finite variance for all $n$, a necessary condition for pointwise consistency. However,
(T{\it iii}) and (T{\it iv}) are stricter on idiosyncratic cross-sectional dependence, since they control all moments of normed linear combinations of idiosyncratic components. Indeed, since the common components~$\mbf X_n$ and $\bm\chi_n$ are recovered by aggregation across the~$n$ elements of $\mbf Y_n$ and $\widehat{\mbf h}_n$, respectively, uniform consistency requires limiting the contribution of the tails of the distribution of cross-sectional averages of idiosyncratic components.
Finally, since factors and factor loadings are not separately identified, we can, without loss of generality, impose the following assumptions, which are just identification constraints (see \citealp{FGLR09} for similar conditions).
\begin{assumption}[I]
\begin{compactenum}[(i)]
\item Denoting by $\mbf P^{X^*}_n$ the $n\times q$ matrix of normalized column eigenvectors corresponding to the $q$ largest eigenvalues of the covariance matrix of $\mbf X_n^*$, put $\mbf H_n := \sqrt n\mbf P^{X^*}_n$ and~$\mbf u_t:={{\mbf P^{X^*}_n}'\mbf X_n^*}/{\sqrt{n}}$;
\item denoting by $\mbf P^{\chi^*}_n$ the $n\times Q$ matrix of normalized eigenvectors corresponding to the $Q$ largest eigenvalues of the covariance matrix of $\bm\chi_n^*$, put $\mbf R_n := \sqrt n\mbf P^{\chi^*}_n$ and $\bm\varepsilon_t:= {{\mbf P^{\chi^*}_n}'\bm\chi_n^*}/{\sqrt{n}}$.
\end{compactenum}
\end{assumption}
In other words, Assumption (I) requires the common factors $\mbf u_t$ ($\bm\varepsilon_t$) to be the (non-normalised) principal components of $\mbf X_n^*$ ($\bm\chi_n^*$). Note that, under Assumption~(I), both the factors and their loadings depend on $n$; their product, however, does not, which is particularly convenient and simplifies the proofs. Other identification constraints are commonly used in principal component analysis (see e.g. \citealp{FLM13}); they do not affect the results below, but lead to much heavier notation.
The consistency properties of the estimated GDFM for the levels as described in steps (\textit{L.i})-(\textit{L.vi}) are as follows.
\begin{proposition} \label{prop:level} Let $\rho_{nT}:=\max\big({B_T}/{\sqrt T}, 1/{B_T}, 1/{\sqrt n}\big)$.
Then, under Assumptions (L1)-(L5), (K), (T), and (I), there exists a $q\times q$ diagonal matrix $\mbf J$ with entries~$\pm 1$ such that
\begin{compactenum}[(a)]
\item $\max_{i=1,\ldots, n}\Vert \widehat{\mbf b}_{ik}'- \mbf b_{ik}'\mbf J\Vert=O_{\rm P}(\rho_{nT})$, for all $k\le \bar k_1$;
\item $\max_{t=1,\ldots, T}\Vert \widehat{\mbf u}_t-\mbf J\mbf u_t\Vert=O_{\rm P}(\rho_{nT}\log T)$;
\item $\max_{i=1,\ldots, n}\vert \widehat{d}_{ik}- d_{ik}\vert=O_{\rm P}(\rho_{nT}\log^2 T)$, for all $k\le \bar k_2$;
\item $\max_{i=1,\ldots, n}\max_{t=1,\ldots, T}\vert \widehat{v}_{it}-v_{it}\vert=O_{\rm P}(\rho_{nT}\log^2 T)$.
\end{compactenum}
\end{proposition}
The proof of parts {\it (a)} and {\it (b)} of Proposition \ref{prop:level} follows directly from \citet{FHLZ17} together with Assumptions (T{\it i}) and (T{\it iii}). However, parts {\it (c)} and {\it (d)} concerning the idiosyncratic components are new results and provide uniform consistency over both time and the cross-section (see also Remark 1 below). In particular, notice that parts {\it (c)} and {\it (d)} of Proposition \ref{prop:level} are proved under Assumption (L2{\it iv}) of a finite-order autoregressive representation for the idiosyncratic component. Relaxing that assumption into possibly infinite-order autoregressive repressentations would require addressing, in the proofs of parts {\it (c)} and {\it (d)}, the issue of truncation errors related to finite-order~AR fitting. Consistency still could be proved, but with rates depending on the rate of decay of the autocovariances of idiosyncratic components, as shown, for example, in \citet{denhaan97}. For simplicity, we do not consider this here.
As for the global consistency properties (after the second estimation step), we need a final condition on the choice of the capping sequence $\kappa_T$ in step (\textit{R}).
\begin{assumption}[R] The sequence $\kappa_T>0$ is such that the sets ${\mathcal T}_{i;nT} :=\big\{t \in\{1,\ldots ,T\}\, \big\vert \, \vert \widehat s_{it}\vert <~\! \kappa_T\big\}$ satisfy $\max_{i=1,\ldots, n}\vert{\mathcal T}_{i;nT}\vert =o_{\rm P}(\sqrt T)$ uniformly in $n$ as $T\to\infty$. Moreover, there exist a positive integer~$\bar T$ and constants~$\varphi>1$ and $0< \underline c \leq \overline c$, independent of $n$, such that $\underline c\le\kappa_T \log^{\varphi} T\le \overline c$ for all~$T>\bar T$.
\end{assumption}
The intuition behind this assumption is as follows.~As shown in Appendix \ref{app:prop2}, an immediate consequence of Proposition \ref{prop:level} is that the volatility proxies are consistently estimated, namely,
\[
\max_{i=1,\ldots, n} \max_{t=1,\ldots,T} \vert\widehat s_{it}- s_{it}\vert = O_{\rm P}(\rho_{nT}\log^2 T),\quad\text{ as } n,T\to\infty.
\]
Now, setting $\kappa_T=0$ in step (\textit{R}), then, due to the log-transform, uniform consistency of $\widehat h_{it}$ becomes problematic when $\widehat s_{it}$ gets ``close to zero''. For this reason, we need $\kappa_T>0$. The set ${\mathcal T}_{i;nT}$ is that of all time points $\{1,\ldots, T\}$ at which~$\widehat s_{it}$ is close to zero, and uniform consistency of $\widehat h_{it}$ for $t\in{\mathcal T}_{i;nT}^c$ straightforwardly follows from uniform consistency of $\widehat s_{it}$. On the other hand, the sets ${\mathcal T}_{i;nT}$ should not contain too many time points, and have cardinality going to zero at appropriate rate---whence Assumption~(R). In particular, we suggest to choose $\kappa_T$ of the order of~$\log ^{-\varphi} T$ for all $i$.
Although we do not have theoretical results justifying this choice of $\kappa_T$ in practice, simulation-based results (see Appendix \ref{sec:app_R}) indicate that, the condition on the cardinality of the sets ${\mathcal T}_{i;nT}$ is indeed satisfied for $\kappa_T$ decreasing logarithmically in $T$.
Consistency of the estimated GDFM for log-volatilities as described in steps (\textit{R}) and (\textit{V.i})-(\textit{V.vi}) then follows.
\begin{proposition}\label{prop:vol}
Let $\tau_{nT}:=\max\big( {B_TM_T}/{\sqrt T}, {M_T}/{\sqrt n}\big)$ and assume that $B_T\ge c T^{1/4}$ for some finite $c>0$. \linebreak Then, under Assumptions (L1)-(L5), (V1)-(V5), (K), (T), (I), and (R), there exists a $Q\times Q$ diagonal matrix $\mbf S$ with entries $\pm 1$ such that
\begin{compactenum}[(a)]
\item $\max_{i=1,\ldots, n}\Vert \widehat{\mbf f}_{ik}'- \mbf f_{ik}'\mbf S\Vert=O_{\rm P}(\tau_{nT} \log^{3+\varphi} T)$ for all $k\le \bar k_1^*$;
\item $\max_{t=1,\ldots, T}\Vert \widehat{\bm \varepsilon}_t-\mbf S\bm \varepsilon_t\Vert=O_{\rm P}(\tau_{nT}\log^{4+\varphi} T)$;
\item $\max_{i=1,\ldots, n}\vert \widehat{g}_{ik}- g_{ik}\vert=O_{\rm P}(\tau_{nT}\log^{5+\varphi} T)$ for all $k\le \bar k_2^*$;
\item $\max_{i=1,\ldots, n}\max_{t=1,\ldots, T}\vert \widehat{\nu}_{it}-\nu_{it}\vert=O_{\rm P}(\tau_{nT}\log^{5+\varphi} T)$.
\end{compactenum}
\end{proposition}
This result, which is new, provides the theoretical foundation for the consistency of the estimators used in \citet{barigozzihallin15a,barigozzihallin15b,barigozzihallin15c} and in this paper. Note that parts {\it (c)} and {\it (d)}, just as parts {\it (c)} and {\it (d)} of Proposition \ref{prop:level}, are proved under Assumption (V2{\it iv}) of a finite-order autoregressive representation for the idiosyncratic components; the same comments as for Proposition \ref{prop:level} apply.
Our results show that, up to logarithmic factors and the bandwidth-related ones, the rates of consistency of our estimators are of order $\min(\sqrt T,\sqrt n)$ as in classical one-step factor models.
The following three technical remarks discuss how our assumptions, in particular Assumptions (T) and (K), affect the consistency rates, and how the effect of those logarithmic and bandwidth-related factors could be controlled further if we were willing to make additional assumptions.
\begin{remark}[Serial dependence of idiosyncratic components]\upshape{Inspection of the proof of part {\it (c)} of Proposition \ref{prop:level} shows that the extra (with respect to part {\it (b)}) $\log T$ factor there is due to terms of the type $T^{-1}\sum_{t=1}^T Z_{it}$. Now, while the cross-sectional dependence of idiosyncratic components is controlled via Assumption (T{\it iii}), we do not impose (beyond weak stationarity) any specific assumption on their serial dependence. However, it is worth noting that, if we made some mild additional mixing assumption controlling that serial dependence, then those terms could be bounded by a Bernstein-type inequality, as for example in Theorem 1 by \citet{MPR11}. Similar comments apply to Proposition \ref{prop:vol} and bounds on the idiosyncratic sums $T^{-1}\sum_{t=1}^T \xi_{it}$. If such additional assumptions were made, the rates in Proposition \ref{prop:level} parts {\it (c)} and {\it (d)} would change to $O_{\rm P}(\rho_{nT} \log T)$, those in Proposition \ref{prop:vol} part {\it (a)} to $O_{\rm P}(\tau_{nT} \log^{1+\varphi} T)$, those in part {\it (b)} to $O_{\rm P}(\tau_{nT}\log^{2+\varphi} T)$, and those in parts {\it (c)} and {\it (d)} to $O_{\rm P}(\tau_{nT}\log^{2+\varphi} T)$.}
\end{remark}
\begin{remark}[Tail behavior]\upshape{
In Section \ref{sec:emp}, we analyze a panel of stock returns, and it is therefore worth discussing how our assumptions relate to the distributional properties of financial data. First, let us stress that it is common, in the financial econometrics literature, to assume Gaussianity of log-volatility proxies \citep[see e.g.][]{ABD02}. This is in agreement with the tail Assumptions~(T{\it ii}) and~(T{\it iv}) since sub-Gaussians tails are lighter than sub-exponentials.
In the Gaussian case, the rates in Proposition \ref{prop:vol} part {\it (a)} would change to $O_{\rm P}(\tau_{nT}\log^{5/2+\varphi} T)$, those in
part {\it (b)} to $O_{\rm P}(\tau_{nT}\log^{3+\varphi} T)$, and those in parts {\it (c)}, and {\it (d)} to $O_{\rm P}(\tau_{nT}\log^{7/2+\varphi} T)$.
Second, Assumption (T{\it i}) straightforwardly generalizes to more general classes of distributions such that, for some finite constants $K_u^*>0$, $K_u^{**}>0$, and~$\vartheta>0$, $\mathrm{P}(|u_{jt}|>\epsilon)\le K_u^* \exp\left(- {\epsilon}^\vartheta/K_{u}^{**}\right)$ for any $\epsilon>0$\linebreak and~$j=1,\ldots,q$; (T{\it iii}) can be generalized similarly for level idiosyncratic components. These distributions are studied in the literature under the name of {\it sub-Weibull distributions} (\citealp{KC18}, and \citealp{VA19}) or {\it semi-exponential} (\citealp{borovkov00}).\footnote{Note that the assumption of a sub-Weibull tail decay is equivalent to the moment condition $(\mathrm E[|u_{jt}|^k])^{1/k}\le Ck^{1/\vartheta}$ for all $k\ge 1$ and some finite $C>0$ (see \citealp[Theorem 2.1]{VA19}); fourth-order moments in that case always exist.} By letting $\vartheta<1$, we could allow for tails, which, although still exponentially decaying, could be heavier than assumed in Assumption~(T), thus accounting for moderately extreme events.
Following the same steps as in Appendix \ref{app:prop1}, it is easily seen that in this case the rates in Proposition \ref{prop:level} part {\it (b)}
would change to~$O_{\rm P}(\rho_{nT}\log^{1/\vartheta} T)$ and those in parts~{\it (c)} and {\it (d)}) to~$O_{\rm P}(\rho_{nT}\log^{2/\vartheta} T)$. As for Proposition \ref{prop:vol}, would we assume a sub-Weibull distribution also in (T{\it ii}) and (T{\it iv}) (with the same value of $\vartheta$), then rates would change to $O_{\rm P}(\tau_{nT}\log^{3/\vartheta+\varphi} T)$ in part {\it (a)}, to $O_{\rm P}(\tau_{nT}\log^{4/\vartheta+\varphi} T)$ in part {\it (b)}, and to~$O_{\rm P}(\tau_{nT}\log^{5/\vartheta+\varphi} T)$ in parts {\it (c)} and {\it (d)}. To conclude, assuming sub-Gaussian tails in (T{\it ii}) and~(T{\it iv}) modifies the rates in part {\it (a)} of Proposition \ref{prop:vol}
into $O_{\rm P}(\tau_{nT}\log^{2/\vartheta+1/2+\varphi} T)$, those in part {\it (b)} in\-to~$O_{\rm P}(\tau_{nT}\log^{2/\vartheta+1+\varphi} T)$, and those in parts {\it (c)} and {\it (d)} into $O_{\rm P}(\tau_{nT}\log^{2/\vartheta+3/2+\varphi} T)$.
Finally, in principle, we also could assume {\it power-law decay}---that is, the existence of finite constants $K_u^*>0$ and $\beta>0$ such that $\mathrm{P}(|u_{jt}|>\epsilon)\le K_u^* \epsilon^{-\beta}$ for any $j=1,\ldots,q$ and $\epsilon>0$; we similarly could generalize~(T{\it iii}) for level idiosyncratic components. We do not explore this possibility in detail, but we notice that , in order to have consistency under this setting, we would need at least $\beta>2$; moreover, the smaller $\beta$, the smaller the range of admissible choices for the bandwidths $B_T$ and $M_T$. Notice however that, in practice, determining the actual values of~$\vartheta$ and $\beta$ is very tricky, and that small values of $\vartheta$ can generate a tail behavior which is comparable to the power-law behavior (see Figure \ref{fig:powerlaw}).}
\end{remark}
\begin{figure}[t!]\caption{\small Comparison of the tails of log-normal (black), power-law (red) and sub-Weibull (blue) probability density functions $f(x)$ (log-scales on both axis).}\label{fig:powerlaw}
\centering \smallskip\noindent
\setlength{\tabcolsep}{.01\textwidth}
\begin{tabular}{@{}c}
\includegraphics[width=.7\textwidth,trim=1cm 1cm 2cm 0cm,clip]{powerlaw.eps} \\
\end{tabular}
\end{figure}
\begin{remark}[Bandwidths and estimation of spectral densities]\upshape{The results in Propositions \ref{prop:level} and \ref{prop:vol} require uniform consistency of the estimated spectral density over all frequencies. For this reason, we have stronger than usual asymptotic constraints on the bandwidths. These could be relaxed if we made stronger assumptions on the shocks. First, notice that in our setting the level shocks are just uncorrelated (see Assumptions (L1{\it i}) and~(L1{\it iii})), and are by no means independent. However, if we are willing to assume the existence, for level shocks, of moments of all orders, then we could apply Theorem 7.7.4 in \cite{brillinger2001}, which would allow us to replace~$B_{T}$ with~$B_T^\epsilon \sqrt {B_T}$ for any~$\epsilon>0$ in the definition of $\rho_{nT}$ in Proposition \ref{prop:level}. Second, we could, in principle, allow for independent shocks on log-volatilities (e.g. assuming Gaussianity, see Remark 1) and therefore make use of Theorem 4 and Section 4.2 in \cite{WZ18}, which would allow us to replace $M_{T}$ with $\sqrt {M_T\log M_T}$ in the definition of $\tau_{nT}$ in Proposition \ref{prop:vol}.
}\end{remark}
}
\section{Conditional prediction intervals}\label{eq:int}
Before describing our prediction intervals, let us summarise here the main notation developed in the previous sections. Given an observed dataset of size $n\times T$, we have, for the levels,
\begin{align}
Y_{it} &= X_{it} + Z_{it} +\mathrm E[Y_{it}], \label{eq:summary1}\\
X_{it} &= \mbf b_{i0}' \mbf u_t + \sum_{k=1}^\infty \mbf b_{ik}' \mbf u_{t-k}:=e_{it} + X_{it|t-1},\quad Z_{it} = d_{i0} v_{it} + \sum_{k=1}^\infty d_{ik} v_{it-k} := v_{it} + Z_{it|t-1}, \nonumber\\
s_{it} &:= e_{it} + v_{it},\ \ i=1,\ldots, n,\ t=1,\ldots, T\nonumber
\end{align}
where $d_{i0}=1$ because of \eqref{eq:var_lev_idio} and, for the log-volatilities,
\begin{align}
h_{it} &:= \log s_{it}^2 = \chi_{it}+ \xi_{it} + \mathrm E[h_{it}], \label{eq:summary2}\\
\chi_{it}&= \mbf f_{i0}' \bm\varepsilon_t + \sum_{k=1}^\infty \mbf f_{ik}' \bm\varepsilon_{t-k}:=\eta_{it} + \chi_{it|t-1},\quad \xi_{it}= g_{i0} \nu_{it} + \sum_{k=1}^\infty g_{ik} \nu_{it-k} := \nu_{it} + \xi_{it|t-1}, \nonumber\\
\omega_{it} &:= \eta_{it} + \nu_{it},\ \ i=1,\ldots, n,\ t=1,\ldots, T\nonumber
\end{align}
where $g_{i0}=1$ because of \eqref{eq:var_vol_idio}.
The optimal one-step-ahead linear predictors of level $Y_{it}$ and log-volatility $h_{it}$ are thus
\begin{equation}\label{predpop}
Y_{it|t-1} :=X_{it|t-1} + Z_{it|t-1}+ \mathrm E[Y_{it}] \quad\text{and }\quad h_{it|t-1} :=\chi_{it|t-1} + \xi_{it|t-1}+ \mathrm E[h_{it}],
\end{equation}
with innovations $s_{it}$ and $\omega_{it}$, respectively. As a consequence, the level innovations are
\[
s_{it} = \exp\big( {h_{it}}/{2}\big) \mbox{sign}(s_{it})= \exp\big({h_{it|t-1}}/{2}\big) \exp\big({\omega_{it}}/{2}\big) \mbox{sign}(s_{it}).
\]
We therefore define a one-step-ahead predictor of the volatilities as
\[
s_{it|t-1} := \exp\big({h_{it|t-1}}/2\big),
\]
with associated ``multiplicative innovations''
\[
w_{it} := \exp\big({\omega_{it}}/{2}\big) \mbox{sign}(s_{it}).
\]
Note, however, that, due to the nonlinear nature of the exponential transformation from $h_i$ to $s_i$, this multiplicative decomposition of volatilities into a predictor and an ``innovation'' does not enjoy (in the space of volatilities) the traditional $L^2$ optimality properties, which only hold for their logarithms (in the space of log-volatilities). This, however, will not be a concern in the quantile-based construction we now describe, due to the fact that the coverage probabilities of a interquantile interval are invariant under continuous monotone transformations: the quantile of $w_{it}$ .
Denoting by $q(\alpha;w_{i})$ the (unconditional) $\alpha$-quantile of $w_{i}:=\{w_{it}\vert t=1,\ldots, T\}$, $i=1,\ldots, n$ (which, by stationarity, does not depend on~$t$), theoretical
lower and upper prediction bounds with confidence level $(1-\alpha )$ and $\alpha\in(0,1)$ are
\begin{equation}\mathcal L_{it|t-1}(\alpha):=Y_{it|t-1} + s_{it|t-1}\, q(\alpha;w_{i}) \label{eq:LUtt1}
\ \text{ and }\
\mathcal U_{it|t-1}(\alpha):=Y_{it|t-1} + s_{it|t-1}\, q(1-\alpha;w_{i}),
\end{equation}
respectively. Note that $Y_{it|t-1}$ lies above $\mathcal L_{it|t-1}(\alpha)$ for $\alpha <{\rm P}[w_{it}\leq 0]$ and lies below $\mathcal U_{it|t-1}(\alpha)$ \linebreak for~$\alpha <1-{\rm P}[w_{it}\leq 0]$.
Prediction intervals with coverage probability $(1-\alpha )$ can be constructed as
\begin{equation}\label{eq:Itt1}
\mathcal I_{it|t-1}(\alpha) :=\big[\,\mathcal L_{it|t-1}(\alpha^-) ,\, \mathcal U_{it|t-1}(\alpha^+)\; \big]
\end{equation}
with $\alpha^\pm<1/2$ \special{color cmyk 0 0 0 1.} and $\alpha^- + \alpha^+ =\alpha $, covering $Y_{it|t-1}$ (see \eqref{eq:Htt1}) provided~that
\begin{equation}
\alpha^-<{\rm P}[w_{it}\leq 0]\quad\text{ and }\quad
\alpha^+ <1-{\rm P}[w_{it}\leq 0].\label{unqtails}
\end{equation}
Clearly, the lower bound $\mathcal L_{it|t-1}(\alpha)$ provides a measure of the Value-at-Risk of level $\alpha$ at time~$t$, which we denote as ${\rm VaR}_{it}(\alpha):=-\mathcal L_{it|t-1}(\alpha)$ (see Section 12.3.1 in \citealp{FZ11} for a review).\footnote{Usually, a Value-at-Risk is reported as a positive quantity. That will be the case with ${\rm VaR}_{it}(\alpha)$ for $\alpha$ small enough. Positive values of~$\mathcal L_{it|t-1}(\alpha)$ are possible, though: in such cases, ${\rm VaR}_{it}(\alpha)$ is defined to be zero by convention (see \citealp{FZ11}, Definition~12.1).}
The advantage of quantile-based prediction intervals of the form \eqref{eq:Itt1} over their conditional heteroske\-dasticity-based competitors stems from the fact that, irrespective of the way $w_{i}$ has been obtained, the conditional $\alpha$-quantiles of $Y_{it}$ (conditional on $Y_{i,t-1}, Y_{i,t-2}, \ldots$) {\it are} of the form \eqref{eq:LUtt1}. This quantile-based approach moreover allows for unequal tails ($\alpha^- \neq \alpha^+ $ in \eqref{unqtails}---hence, distinct attitudes towards losses and gains) and automatically takes into account the typical skewness of financial data distributions.
In practice, the model is estimated from a $n\times T$ observed panel; the empirical counterparts of~${Y}_{i,T+1|T} $ and~${h}_{i,T+1|T}$ for $i=1,\ldots,n$ are\vspace{-3mm}
\[
\widehat{Y}_{i,T+1|T} = \widehat{X}_{i,T+1|T}+\widehat{Z}_{i,T+1|T}+ \frac 1 T\sum_{t=1}^T Y_{it}= \sum_{k=1}^{\bar k_1} \widehat{\mbf b}_{ik}'\widehat{\mbf u}_{T-k+1} + \sum_{k=1}^{\bar k_2} \widehat{d}_{ik}\widehat{v}_{i,T-k+1}+ \frac 1 T\sum_{t=1}^T Y_{it}\vspace{-3mm} \]
and\vspace{-3mm}
\[\widehat{h}_{i,T+1|T}=\widehat{\chi}_{i,T+1|T}+\widehat{\xi}_{i,T+1|T}+ \frac 1 T\sum_{t=1}^T \widehat h_{it}= \sum_{k=1}^{\bar k_1^*} \widehat{\mbf f}_{ik}'\widehat{\bm\varepsilon}_{T-k+1} + \sum_{k=1}^{\bar k_2^*} \widehat{g}_{ik}\widehat{\nu}_{i,T-k+1}+ \frac 1 T\sum_{t=1}^T \widehat h_{it},
\vspace{-1mm}\]
and we accordingly define
$
\widehat{s}_{i,T+1|T}:=\exp\big({\widehat{h}_{i,T+1|T}}/{2}\big)
$;
based on the estimates $
\widehat s_{it}$ and~$
\widehat{\omega}_{it}$ of $s_{it}$ and~$\omega_{it}$, let~$\widehat w_{it}:=\exp\big({\widehat {\omega}_{it}}/2\big)\text{sign}(\widehat s_{it})$.
For any $i$, denote by $\widehat w_{i(1)},\ldots, \widehat w_{i(T)}$ the order statistic of $\widehat w_{i1},\ldots, \widehat w_{iT}$; the empirical quantile~$w_{i(\lceil T\alpha\rceil)}$ then can be used as an estimator of $q(\alpha;w_{i})$. Empirical versions of the prediction limits and intervals \eqref{eq:LUtt1} and \eqref{eq:Itt1} are
\begin{equation}\nonumber
\widehat{\mathcal L}_{i,T+1|T}(\alpha):=\widehat Y_{i,T+1|T} + \widehat s_{i,T+1|T}\, \widehat w_{i(\lceil T\alpha\rceil)},\quad
\widehat{\mathcal U}_{i,T+1|T}(\alpha):=\widehat Y_{i,T+1|T} + \widehat s_{i,T+1|T}\, \widehat w_{i(\lceil T(1-\alpha)\rceil)}
\end{equation}
and
\begin{equation}\label{eq:Itt1hat}
\widehat{\mathcal I}_{i,T+1|T}(\alpha) :=\big[\widehat{\mathcal L}_{i,T+1|T}(\alpha^-),\widehat{\mathcal U}_{i,T+1|T}(\alpha^+)\big]
\end{equation}
with $\alpha^\pm<1/2$ and $\alpha^-+\alpha^+=\alpha\in(0,1)$. A schematic description of this procedure is given in Algorithm~3.
If the $w_{it}$'s were i.i.d.\ instead of weak white noise, the convergence (for given $\alpha^-$ and $\alpha^+$, without rates) of~\eqref{eq:Itt1hat} to~\eqref{eq:Itt1} would follow from the fact that, as a consequence of the consistent estimation of the GDFMs for levels and volatilities, for any~$n_0$ and $T_0$, $\max_{1\leq i\leq n_0}\max_{1\leq t\leq T_0}\vert\widehat w_{it}-w_{it}\vert$ converges to zero as $n$ and $T$ tend to infinity.
Then, the difference between the empirical quantile of order~$\alpha$ computed from $\{\widehat w_{1t},\ldots,\widehat w_{iT_0}\}$ and the empirical quantile of order~$\alpha$ computed from the unobservable $\{w_{i1},\ldots,w_{iT_0}\}$ is~$o_{\rm P}(1)$ for given $1\leq i\leq n_0$ as $n$ and~$T$ tend to infinity. Now, for given $i$, were the $w_{it}$'s i.i.d., the empirical $\alpha$-quantile computed from $\{w_{i1},\ldots,w_{iT_0}\}$ is, for~$T_0$ large enough, arbitrarily close to its theoretical counterpart $q(\alpha; w_i)$ with probability arbitrarily close to one. The same conclusion extends to the present case where the $w_{it}$'s are stationary and uncorrelated provided that they satisfy some additional mild ergodicity or mixing assumption. The literature on Glivenko-Cantelli and quantile consistency under ergodicity and mixing is abundant, and we will not proceed with imposing any specific mixing conditions here which anyway hardly can be checked from the data. The reader may like to refer to Theorem 3.1 in \citet{FZ19} for details.
Once prediction regions have been constructed, it is important to evaluate their actual coverage performance. For this, it is useful to define the conditional coverage indicators---namely, for prediction intervals $ \widehat{\mathcal I}_{i,T+1|T}(\alpha)$,
\begin{equation}\label{eq:Htt1}
\widehat{\mathcal H}_{i,T+1|T}(\alpha):=\mathbb I\big(Y_{i,T+1}\in \widehat{\mathcal I}_{i,T+1|T}(\alpha)\big).
\end{equation}
For a given $i$, we say that $ \widehat{\mathcal I}_{i,T+1|T}(\alpha)$ provides the correct coverage if
\begin{equation}\nonumber
\mathrm{P} (Y_{i,T+1}\in \widehat{\mathcal I}_{i,T+1|T}(\alpha) | Y_{i,T},\ldots ,Y_{i1} )=\mathrm E[\widehat{\mathcal H}_{i,T+1|T}(\alpha) | Y_{i,T},\ldots ,Y_{i1}]=(1-\alpha),
\end{equation}
which is equivalent (see e.g. Lemma 1 in \citealp{christoffersen1998}) to the hypothesis that
\begin{equation}\label{eq:null}
\widehat{\mathcal H}_{i,T+1|T}(\alpha)\stackrel{iid}{\sim}\mbox{Bernoulli}(1-\alpha).
\end{equation}
That hypothesis can be tested against alternatives of insufficient coverage probability values, against non-sharp prediction limits, or against alternatives of serial dependence. We refer to Section~\ref{backtestSec} for details and implementation.
\begin{algorithm}[t]
\DontPrintSemicolon
\footnotesize{
\KwInput{
data in levels $\mbf Y$, $\alpha^-\in[0,1/2]$ and $\alpha^+\in[0,1/2]$ such that the confidence level is $\alpha=\alpha^++\alpha^-\in(0,1)$\\
from Algorithm 1: common level shocks $\widehat{\mathbf u}$ of size $q\times T$ and $\widehat{\mathbf e}$ of size $n\times T$, idiosyncratic level shocks $\widehat{\mbf v}$ of size $n\times T$, common level impulse responses $\widehat{\mbf B}(L)$ of size $n\times q\times \bar k_1$, idiosyncratic level impulse responses $\widehat{\mbf D}(L)$ of size $n\times n\times \bar k_2$\\
from Algorithm 2: log-volatility proxy $\mbf h$ of size $n\times T$, common log-volatility shocks $\widehat{\bm\varepsilon}$ of size $Q\times T$ and $\widehat{\bm\eta}$ of size $n\times T$, idiosyncratic log-volatility shocks $\widehat{\bm \nu}$ of size $n\times T$, common level impulse responses $\widehat{\mbf F}(L)$ of size $n\times Q\times\bar k_1^*$, idiosyncratic level impulse responses $\widehat{\mbf G}(L)$ of size $n\times n\times \bar k_2^*$}
\KwOutput{lower bounds of conditional prediction interval $\widehat{\mathcal L}_{1,T+1|T}(\alpha^-),\ldots,\widehat{\mathcal L}_{n,T+1|T}(\alpha^-)$\\
upper bounds of conditional prediction interval $\widehat{\mathcal U}_{1,T+1|T}(\alpha^+),\ldots,\widehat{\mathcal U}_{n,T+1|T}(\alpha^+)$}
\medskip
Compute $\bar {\mbf Y}$ the sample mean of levels $\mbf Y$\,\smallskip
Compute the one-step-ahead prediction of common and idiosyncratic components of levels $\widehat{\mbf X}_{T+1|T}=\sum_{k=1}^{\bar k_1} \widehat{\mbf B}_k\widehat{\mathbf u}_{T-k+1}$\,\smallskip
Compute the one-step-ahead prediction of idiosyncratic component of levels $\widehat{\mbf Z}_{T+1|T}=\sum_{k=1}^{\bar k_2} \widehat{\mbf D}_k\widehat{\mathbf v}_{T-k+1}$\,\smallskip
Compute the one-step-ahead prediction of levels $\widehat{\mbf Y}_{T+1|T}=\widehat{\mbf X}_{T+1|T}+\widehat{\mbf Z}_{T+1|T}+\bar {\mbf Y}$ such that $\widehat{\mbf Y}_{T+1|T}=(\widehat{Y}_{1,T+1|T}\ldots \widehat{Y}_{n,T+1|T})^\prime$\,\smallskip
Compute $\bar {\widehat{\mbf h}}$ the sample mean of log-volatilities $\widehat{\mbf h}$\,\smallskip
Compute the one-step-ahead prediction of common component of log-volatilities $\widehat{\bm\chi}_{T+1|T}=\sum_{k=1}^{\bar k_1^*} \widehat{\mbf F}_k\widehat{\bm\varepsilon}_{T-k+1}$\,\smallskip
Compute the one-step-ahead prediction of idiosyncratic component of log-volatilities $\widehat{\bm\xi}_{T+1|T}=\sum_{k=1}^{\bar k_2^*} \widehat{\mbf G}_k\widehat{\bm\nu}_{T-k+1}$\,\smallskip
Compute the one-step-ahead prediction of log-volatilities $\widehat{\mbf h}_{T+1|T}=\widehat{\bm\chi}_{T+1|T}+\widehat{\bm\xi}_{T+1|T}+\bar {\widehat{\mbf h}}$\,\smallskip
Compute the one-step-ahead prediction of volatilities $\widehat{\mbf s}_{T+1|T}=\exp(\widehat{\mbf h}_{T+1|T}/2)$ such that $\widehat{\mbf s}_{T+1|T}=(\widehat{s}_{1,T+1|T}\ldots \widehat{s}_{n,T+1|T})^\prime$\,\smallskip\,\smallskip
Compute the log-volatility innovations $\widehat{\bm\omega}_t=\widehat{\bm\eta}_t+\widehat{\bm\nu}_t$ for $t=1,\ldots, T$\,\smallskip
Compute the volatility proxy $\widehat{\mbf s}_{t}=\exp(\mbf h_t/2)$ or equivalently $\widehat{\mbf s}_{t}=\widehat{\mbf e}_t+\widehat{\mbf v}_t$ for $t=1,\ldots, T$\,\smallskip
Compute the volatility innovations $\widehat{\mbf w}_{t}=\exp(\widehat{\bm\omega}_t/2)\mbox{sign}(\widehat{\mbf s}_t)$
such that $\widehat{\mbf w}_{t}=(\widehat{w}_{1t}\ldots \widehat{w}_{nt})^\prime$ for $t=1,\ldots, T$\,\smallskip
\,\smallskip
\For{$i\leftarrow 1$ \KwTo $n$}{Compute the order statistics $w_{i(\lceil T\alpha^-\rceil)}$ and $w_{i(\lceil T(1-\alpha^+)\rceil)}$ of $w_i$\,\smallskip
Compute the lower bound $\widehat{\mathcal L}_{i,T+1|T}(\alpha^-)=\widehat Y_{i,T+1|T} + \widehat s_{i,T+1|T}\, \widehat w_{i(\lceil T\alpha^-\rceil)}$\,\smallskip
Compute the upper bound $\widehat{\mathcal U}_{i,T+1|T}(\alpha^+)=\widehat Y_{i,T+1|T} + \widehat s_{i,T+1|T}\, \widehat w_{i(\lceil T(1-\alpha^+)\rceil)}$\,\smallskip
}
}
\caption{\small Estimation of conditional prediction intervals}
\end{algorithm}
\setcounter{equation}{0}
\section{Simulation study}\label{sec:sim}
\subsection{Setup}
To study the performance of our estimator on finite samples, we simulate data ($\mathcal M$ replications) according to the model described in \eqref{eq:summary1}-\eqref{eq:summary2}.
For each Monte Carlo replication $m=1,\ldots, \mathcal M$ and for given values of $n,T,q$, and $Q$, we first simulate a multiplicative factor model for the volatilities which in turn implies a factor structure also for the levels. The common component of the log-volatilities is generated as
\[
\bm\chi_{nt,m}:=(\mbf M_{n,m}(L))^{-1} \mbf R_{n,m}\bm\varepsilon_{t,m}, \quad t=1,\ldots, T,
\]
where $\bm\varepsilon_{t,m}\stackrel{iid}{\sim} N(\mbf 0_Q,\mbf I_Q)$, $\mbf R_{n,m}$ is $n\times Q$ with entries $[\mbf R_{n,m}]_{ij}\stackrel{iid}{\sim} N(0,1)$ and rescaled such that $\mbf R_{n,m}^\prime \mbf R_{n,m}=n$,\linebreak and $\mbf M_{n,m}(L)=\mbf I_n-\sum_{k=1}^3\mbf M_{kn,m}L^k$ where the coefficients $\mbf M_{kn,m}$ are diagonal $n\times n$ matrices with entries~$[\mbf M_{kn,m}]_{ij}\stackrel{iid}{\sim} N(0,1)$ and rescaled in such a way that $\det(\mbf M_{n,m}(z))\ne 0$ for $|z|\le 1$.\footnote{In particular, when looking at simulated data 25\% of the total $3n^2$ roots are found to be in the range $(0.7,1)$, thus accounting for high persistence in log-volatilities, see also Table \ref{tab:dgp1} below.} Then, we generate the process
\[
\bm\xi_{nt,m}^*:=(\mbf P_{n,m}^*(L))^{-1} \bm\nu_{nt,m}^*, \quad t=1,\ldots, T,
\]
where $\bm\nu_{nt,m}^*\stackrel{iid}{\sim} N(\mbf 0_n,\bm \Sigma_{n,m})$, with $\bm \Sigma_{n,m}$ a Toeplitz matrix with entries $[\bm \Sigma_{n,m}]_{ij}:=0.5^{|i-j|}$, if $|i-j|\le 2$ and zero otherwise, and $\mbf P_{n,m}^*(L)$ generated in the same way as $\mbf M_{n,m}(L)$. Denoting by~$\xi^*_{it,m}$ the $i$th element of~$\bm\xi^*_{it,m}$, we rescale it into $\xi^{**}_{it,m}:= \xi^*_{it,m} [\text{\rm Var}( \chi_{it,m})/\{2\text{\rm Var}( \xi^*_{it,m})\}]^{1/2}$ so that the signal-to-noise ratio is~2.
Define\vspace{-1mm}
\begin{align}
e_{it,m}^*&:=\exp(\chi_{it,m}/2)\ \pi_{it,m},\;\mbox{ and }\; v_{it,m}^*:=\exp(\chi_{it,m}/2)\ \exp(\xi^{**}_{it,m}/2) \pi_{it,m},\quad t=1,\ldots, T,\ i=1,\ldots, n,\nonumber
\end{align}
where $\pi_{it,m}=\pm1$ with equal probabilities 0.5 and
$\chi_{it,m}$ is the $i$th element of~$\bm\chi_{t,m}$. The volatility and log-volatility proxies then are
\begin{align}
s_{it,m}^2&:=(e_{it,m}^*+v_{it,m}^*)^2 = \exp(\chi_{it,m})\left[1+\exp(\xi^{**}_{it,m})+2\exp(\xi^{**}_{it,m}/2) \right],\nonumber\\
h_{it,m}&:=\log(s_{it,m}^2)= \chi_{it,m}+\log \left[1+\exp(\xi^{**}_{it,m})+2\exp(\xi^{**}_{it,m}/2) \right],\; t=1,\ldots, T,\ i=1,\ldots, n,\nonumber
\end{align}
from which we see that, since each $\chi_{i,m}$ is driven by the $Q$-dimensional vector of shocks $\bm\varepsilon_m$, it has the role of common log-volatility, while the $n$ shocks $\bm\nu^*_{n,m}$ have only an idiosyncratic role.
Letting $\mbf V$ be the $q$ normalized eigenvectors corresponding to the $q$ largest eigenvalues of the sample covariance of the vector $\mbf e_{nt,m}^*:=(e_{1t,m}^*\ldots e_{nt,m}^*)^\prime$, we build the level shocks as
\begin{align}
\mbf e_{nt,m}&:=\mbf V\mbf V^\prime\mbf e^*_{nt,m},\;\mbox{ and }\;\mbf v_{nt,m}:=\mbf V_\perp\mbf V_\perp^\prime\mbf e^*_{nt,m}+\mbf v_{nt,m}^*,\quad t=1,\ldots, T,\nonumber
\end{align}
where $\mbf V_\perp$ is $n\times (n-q)$ such that $\mbf V_\perp^\prime \mbf V=\mbf 0_{(n-q)\times q}$, and $\mbf v_{nt,m}^*:=(v_{1t,m}^*\ldots v_{nt,m}^*)^\prime$. Note that, by construction, the elements $e_{it,m}$ and $v_{it,m}$ of the vectors $\mbf e_{nt,m}$ and $\mbf v_{nt,m}$ are such that $(e_{it,m}+v_{it,m})= (e_{it,m}^*+ v_{it,m}^*)$: therefore, we can also write $s_{it}^2=(e_{it,m}+v_{it,m})^2$.
Finally, we generate the vectors of common and idiosyncratic components of the levels as
\begin{align}
\mbf X_{nt,m} &:= (\mbf I_n-\mbf A_{n,m}L)^{-1}\mbf e_{nt,m}, \;\mbox{ and }\;\mbf Z_{nt,m}:= (\mbf I_n-\mbf C_{n,m}L)^{-1}\mbf v_{nt,m}, \quad t=1,\ldots, T\nonumber
\end{align}
where $\mbf A_{n,m}$ is a diagonal $n\times n$ matrix with entries $[\mbf A_{n,m}]_{ij}\stackrel{iid}{\sim} U[-0.3,0.7]$, and $\mbf C_{n,m}$ is generated in the same way but with entries from a uniform distribution over $[\mbf C_{n,m}]_{ij}\stackrel{iid}{\sim} U[-0.5,0.5]$; since these matrices are diagonal, the autoregressive models for $\mbf X_{n,m}$ and $\mbf Z_{n,m}$ are causal. The panel of
levels
then is generated \linebreak as~$\mbf Y_{nt,m}:=\mbf X_{nt,m}+\mbf Z_{nt,m}$.
In our numerical study, we let $n\in\{100,200\}$, $T\in\{200,500,1000\}$, and either $q=1$ and $Q=1$, $q=3$ and $Q=2$ (as in the empirical application of the next section), or $q=2$ and $Q=3$. For each configuration considered, we simulate and estimate the model $\mathcal M=200$ times.
It has to be noticed that the data-generating process we are considering is similar to a stochastic volatility model. To illustrate the properties of the generated data, we report in Table \ref{tab:dgp1} the autocorrelations up to lag 10 of $h_{i,m}$, $s_{i,m}$, $e_{i,m}$, $v_{i,m}$, $X_{i,m}$, $X_{i,m}^2$, $Z_{i,m}$, $Z_{i,m}^2$, $Y_{i,m}$, and $Y_{i,m}^2$, averaged over all $\mathcal M$ replications and over all~$n$ series, and when $n=200$, $T=1000$.
It can be seen that log-volatilities $h_{i,m}$ and volatilities $s_{i,m}$ have high persistence, while, due to the way they are generated, the shocks $e_{i,m}$ and $v_{i,m}$ display no linear serial dependence, i.e. are weak white noises. Turning to the kurtosis of the level shocks reported in the left panel of Table \ref{tab:kurt}, these display heavy tails (especially the common ones) for the case $q=1$ and $Q=1$, while the kurtosis tends to decrease when increasing $Q$, possibly due to the aggregation of shocks in generating the common components of the log-volatility~$\chi_{i,m}$. Similar comments apply to the absolute values of skewness reported in the right panel of Table~\ref{tab:kurt}: especially in the case $q=1$ and $Q=1$, the common shocks display a high degree of asymmetry. Because of these features of the simulated data the case $q=1$ and $Q=1$ is particularly interesting to study to assess the performance of our estimators when dealing with heavy-tailed and skewed data.
\begin{table}[t!]
\centering
\caption{\small Autocorrelations of simulated variables. Average values over all $n$ series and all $\mathcal M$ replications for $n=200$, $T=1000$, and $\mathcal M=200$. Values outside the $[\pm 1.96/\sqrt T]=[\pm 0.0620]$ interval are starred.}\label{tab:dgp1}
\vskip .2cm
\footnotesize
\begin{tabular}{l | ccc ccc }
\hline\hline
$q=1$ & \multicolumn{6}{|c}{lag}\\
$Q=1$& 1&2&3&4&5&6\\
\hline\hline
$h_{i,m}$ & 0.2967* & 0.2856* & 0.1178* & 0.1691* & 0.0528 & 0.1142* \\
$s_{i,m}$ & 0.3082* & 0.2743* & 0.1201* & 0.1426* & 0.0552 & 0.0758* \\
\hline
$e_{i,m}$ & -0.0378 & 0.0696* & -0.0565 & 0.0089 & -0.0099 & -0.0283 \\
$e^2_{i,m}$ & 0.0326 & 0.0495 & 0.0020 & 0.0035 & 0.0031 & 0.0042 \\
\hline
$v_{i,m}$ & -0.0023 & -0.0017 & -0.0035 & 0.0005 & -0.0015 & -0.0028 \\
$v^2_{i,m}$ & 0.2961* & 0.2674* & 0.1196* & 0.1445* & 0.0565 & 0.0771* \\
\hline
$X_{i,m}$ & 0.1753* & 0.1770* & 0.0096 & 0.0302 & -0.0032 & -0.0262 \\
$X^2_{i,m}$ & 0.0989* & 0.0791* & 0.0068 & 0.0023 & 0.0016 & -0.0037 \\
\hline
$Z_{i,m}$ & 0.0020 & 0.0783* & -0.0033 & 0.0118 & -0.0018 & -0.0002 \\
$Z^2_{i,m}$ & 0.2877* & 0.2260* & 0.1054* & 0.1200* & 0.0522 & 0.0622* \\
\hline
$Y_{i,m}$ & 0.1174* & 0.1439* & 0.0051 & 0.0219 & -0.0031 & -0.0172 \\
$Y^2_{i,m}$ & 0.0994* & 0.0798* & 0.0073 & 0.0023 & 0.0018 & 0.0029 \\
\hline
\hline
$q=3$ & \multicolumn{6}{|c}{lag}\\
$Q=2$& 1&2&3&4&5&6\\
\hline\hline
$h_{i,m}$ & 0.2654* & 0.2826* & 0.1183* & 0.1611* & 0.0508 & 0.1272* \\
$s_{i,m}$ & 0.2757* & 0.2607* & 0.1237* & 0.1305* & 0.0511 & 0.0913* \\
$e_{i,m}$ & -0.0089 & -0.0690* & -0.0028 & 0.0304 & -0.0353 & -0.0071 \\
$e^2_{i,m}$ & 0.1939* & 0.0721* & 0.0034 & 0.0166 & 0.0077 & 0.0113 \\
$v_{i,m}$ & -0.0010 & -0.0011 & -0.0030 & 0.0002 & -0.0038 & -0.0030 \\
$v^2_{i,m}$ & 0.2635* & 0.2515* & 0.1212* & 0.1283* & 0.0503 & 0.0907* \\
$X_{i,m}$ & 0.1977* & 0.0592 & 0.0465 & 0.0492 & -0.0142 & 0.0002 \\
$X^2_{i,m}$ & 0.2545* & 0.0822* & 0.0161 & 0.0222 & 0.0073 & 0.0101 \\
$Z_{i,m}$ & -0.0172 & 0.0814* & -0.0052 & 0.0119 & -0.0045 & -0.0013 \\
$Z^2_{i,m}$ & 0.2708* & 0.2133* & 0.1076* & 0.1059* & 0.0502 & 0.0764* \\
$Y_{i,m}$ & 0.1227* & 0.0671* & 0.0285 & 0.0374 & -0.0123 & -0.0003 \\
$Y^2_{i,m}$ & 0.2166* & 0.0909* & 0.0237 & 0.0242 & 0.0079 & 0.0009 \\
\hline
\hline
$q=2$ & \multicolumn{6}{|c}{lag}\\
$Q=3$& 1&2&3&4&5&6\\
\hline\hline
$h_{i,m}$ & 0.2730* & 0.2692* & 0.1254* & 0.1717* & 0.0494 & 0.1293* \\
$s_{i,m}$ & 0.2375* & 0.2348* & 0.1036* & 0.1481* & 0.0238 & 0.0885* \\
\hline
$e_{i,m}$ & -0.0221 & 0.0127 & -0.0204 & -0.0001 & -0.0071 & -0.0067 \\
$e^2_{i,m}$ & 0.0025 & 0.0640* & 0.0222 & 0.0616 & 0.0001 & 0.0134 \\
\hline
$v_{i,m}$ & -0.0026 & -0.0002 & 0.0000 & 0.0026 & -0.0011 & -0.0035 \\
$v^2_{i,m}$ & 0.2330* & 0.2313* & 0.1015* & 0.1456* & 0.0229 & 0.0879* \\
\hline
$X_{i,m}$ & 0.1835* & 0.1254* & 0.0364 & 0.0277 & 0.0079 & 0.0008 \\
$X^2_{i,m}$ & 0.1159* & 0.1020* & 0.0388 & 0.0377 & 0.0131 & 0.0057 \\
\hline
$Z_{i,m}$ & 0.0317 & 0.0872* & 0.0032 & 0.0156 & -0.0008 & -0.0007 \\
$Z^2_{i,m}$ & 0.2509* & 0.1918* & 0.0969* & 0.1227* & 0.0293 & 0.0711* \\
\hline
$Y_{i,m}$ & 0.1313* & 0.1121* & 0.0253 & 0.0220 & 0.0040 & 0.0012 \\
$Y^2_{i,m}$ & 0.1047* & 0.0869* & 0.0426 & 0.0425 & 0.0043 & 0.0133 \\
\hline
\hline
\end{tabular}
\end{table}
\begin{table}[t!]
\centering
\caption{\small Kurtosis and absolute value of skewness of simulated common level shocks $e_{i,m}$ and idiosyncratic level shocks~$v_{i,m}$. Maximum and average values over all $n$ series and all $\mathcal M$ replications for $n=200$, $T=1000$, and~$\mathcal M=200$.}\label{tab:kurt}
\vskip .2cm
\footnotesize
\begin{tabular}{l | cc | cc | cc|| cc | cc | cc}
\hline
\hline
& \multicolumn{6}{c||}{kurtosis}& \multicolumn{6}{c}{skewness}\\
&\multicolumn{2}{c|}{$q=1$, $Q=1$}&\multicolumn{2}{c|}{$q=3$, $Q=2$}&\multicolumn{2}{c||}{$q=2$, $Q=3$}&\multicolumn{2}{|c|}{$q=1$, $Q=1$}&\multicolumn{2}{c|}{$q=3$, $Q=2$}&\multicolumn{2}{c}{$q=2$, $Q=3$}\\
& max.& aver.&max.& aver.&max.& aver.& max.& aver.&max.& aver.&max.& aver.\\
\hline
$e_{i,m}$& 161.40 & 83.53&67.60& 10.94&31.86& 5.60&7.88& 0.28&4.08& 0.03&2.53& 0.02\\
$v_{i,m}$& 15.03 & 3.02& 12.68& 3.02&9.36& 3.01&1.14& 0.01&1.07& 0.01&0.85& 0.01\\
\hline
\hline
\end{tabular}
\end{table}
Furthermore, notice that $\mbf e_{n,m}$, by construction, is a singular vector (as it should be) and has the role of a common level innovation. Moreover, the elements of $\mbf v_{n,m}$, in general, are cross-sectionally dependent. As a consequence, both $\mbf Y_{n,m}$ and $\mbf h_{n,m}$ have an approximate dynamic factor structure. In Figure \ref{fig:eval_dgp} we show scree-plots with the ten largest eigenvalues of the zero-frequency sample spectral density matrices of $\mbf Y_{n,m}$ (blue crosses), and $\mbf h_{n,m}$ (red circles), normalized by the largest zero-frequency eigenvalue, averaged over all $\mathcal M$ realisations, when $n=200$, $T=1000$, $q=3$, and $Q=2$.
\begin{figure}[t!]\caption{\small Normalized eigenvalues of the zero-frequency spectral densities of simulated data for $n=200$ and $T=1000$. Blue crosses: levels, $\mbf Y_{n,m}$;
red circles: log-volatilities, $\mbf h_{n,m}$.}\label{fig:eval_dgp}
\centering \smallskip\noindent
\setlength{\tabcolsep}{.01\textwidth}
\begin{tabular}{@{}c}
\includegraphics[width=.7\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{eval_spec0_32.eps} \\
\end{tabular}
\end{figure}
\subsection{Results}
For each replication, we estimate the model as described in Section \ref{sec:est}. The capping constants $\kappa_T$ and the bandwidths $B_T$ and $M_T$ involved in the estimation of the spectral density are chosen as in the empirical analysis of the next section. Specifically, we let $\kappa_T\in\{0, 0.2,0.4\}$, while the bandwidths values are~$B_T=2$ and $M_T=10$ for $T=200$, $B_T=2$ and $M_T=15$ for $T=500$, $B_T=2$ and $M_T=20$ for~$T=1000$ (see Appendix~\ref{app:sim} for results based on other values). Once we obtain estimated common components $\widehat{X}_{i,m}$ for the levels and $\widehat{\chi}_{i,m}$ for the log-volatilities, we compute the global error measures
\begin{align}
&MSE^X =\frac1{\mathcal MnT}\sum_{m=1}^{\mathcal M}{\sum_{i=1}^n\sum_{t=1}^T ({X}_{it,m}-\widehat{X}_{it,m})^2}, \quad
MSE^\chi=\frac1{\mathcal MnT}\sum_{m=1}^{\mathcal M}\sum_{i=1}^n\sum_{t=1}^T ({\chi}_{it,m}-\widehat{\chi}_{it,m})^2,\nonumber\\
&MAD^X =\frac1{\mathcal MnT}\sum_{m=1}^{\mathcal M}{\sum_{i=1}^n\sum_{t=1}^T \vert{X}_{it,m}-\widehat{X}_{it,m}\vert}, \quad
MAD^\chi=\frac1{\mathcal MnT}\sum_{m=1}^{\mathcal M}\sum_{i=1}^n\sum_{t=1}^T \vert{\chi}_{it,m}-\widehat{\chi}_{it,m}\vert\nonumber
\end{align}
and the maximal errors over all realizations:
\begin{align}
&MAX^X ={ \max_{i=1,\ldots, n}\max_{t=1,\ldots, T}\max_{m=1,\ldots, \mathcal M} \vert{X}_{it,m}-\widehat{X}_{it,m}\vert}, \nonumber\\
&MAX^\chi={ \max_{i=1,\ldots, n}\max_{t=1,\ldots, T}\max_{m=1,\ldots,\mathcal M} \vert{\chi}_{it,m}-\widehat{\chi}_{it,m}\vert}.\nonumber
\end{align}
Notice that the error in the estimation of the common component $X_{i,m}$ of the levels (first step of the estimation procedure) has already been studied in \citet{FHLZ17} and \citet{FGLS18}. We therefore consider it as the benchmark error with respect to which the performance of the second estimation step, which is the novelty of this paper, is to be compared. Results are provided in Table~\ref{tab:sim1}. We note that MSE and MAD in the second step tend to be about 1.5 times higher than in the first step, which is not unexpected as first- and second- step errors typically cumulate in a two-stage procedure. However, when turning to MAX, this is no longer the case, since levels in our data-generating process display heavier tails than log-volatilities---in line with the typical behavior of daily stock returns and their volatilities. Increasing $n$ and $T$ improves the performance of all estimators; the role of $n$, in that respect, seems to be the main one---a manifestation of the ``blessing of dimensionality". On the other hand increasing $Q$ the number of common log-volatility shocks, tends to make estimation of the second step harder, but still results are in line with the case $Q=1$. Capping has an effect in controlling the maximum error but does not affect the MSE and MAD results much.
To illustrate the good performances of our method, in Figure~\ref{fig:sim1014} we show, for one replication, the estimated (in red) and simulated (in blue) common components of levels, and of volatilities, respectively, for $n=200$, $T=1000$, $q=1$, and $Q=1$ (which is the case exhibiting the heaviest tails), setting $\kappa_T=0.2$. The choice of bandwidths adopted seems to work quite well, and, comparing to alternative choices considered in Appendix \ref{app:sim}, it can be shown that $M_T$ must be large enough to capture the persistence in log-volatilities, while lower values of $B_T$ are enough for levels and do not affect much the second step of estimation.
\begin{table}[t!]
\centering
\caption{\small Simulation results. MSEs and MADs for common components. Bandwidths are $B_T=2$ and $M_T=10$ for~$T=200$; $B_T=2$ and $M_T=15$ for $T=500$; $B_T=2$ and $M_T=20$ for $T=1000$. }\label{tab:sim1}
\vskip .2cm
\footnotesize
\begin{tabular}{l l | cc | cc | cc }
\multicolumn{8}{c}{$q=1$, $Q=1$}\\
\hline\hline
&& \multicolumn{2}{|c|}{$T=200$}& \multicolumn{2}{|c|}{$T=500$}& \multicolumn{2}{|c}{$T=1000$}\\
&& $n=100$ & $n=200$& $n=100$ & $n=200$& $n=100$ & $n=200$\\
\hline\hline
$MSE^X$&& 0.215 & 0.219 & 0.168 & 0.164 & 0.125 & 0.154 \\
$MSE^\chi$&$\kappa_T=0$& 0.368 & 0.321 & 0.302 & 0.247 & 0.251 & 0.241 \\
$MSE^\chi$&$\kappa_T=0.2$ & 0.369 & 0.324 & 0.276 & 0.247 & 0.240 & 0.230 \\
$MSE^\chi$ & $\kappa_T=0.4$ & 0.377 & 0.334 & 0.279 & 0.238 & 0.238 & 0.230 \\
\hline
$MAD^X$&& 0.278 & 0.245 & 0.255 & 0.226 & 0.234 & 0.228 \\
$MAD^\chi$&$\kappa_T=0$& 0.442 & 0.401 & 0.394 & 0.346 & 0.360 & 0.342 \\
$MAD^\chi$&$\kappa_T=0.2$ & 0.436 & 0.395 & 0.367 & 0.346 & 0.344 & 0.323 \\
$MAD^\chi$ & $\kappa_T=0.4$ & 0.437 & 0.397 & 0.364 & 0.324 & 0.337 & 0.317 \\
\hline
$MAX^X$&& 10.797 & 13.560 & 17.113 & 16.352 & 14.405 & 19.757 \\
$MAX^\chi$&$\kappa_T=0$& 6.587 & 6.709 & 8.996 & 6.270 & 7.462 & 7.912 \\
$MAX^\chi$&$\kappa_T=0.2$ & 8.259 & 7.247 & 7.967 & 6.270 & 9.328 & 8.431 \\
$MAX^\chi$ & $\kappa_T=0.4$ & 8.987 & 8.402 & 8.295 & 10.223 & 9.821 & 8.960 \\
\hline\hline
\\
\multicolumn{8}{c}{$q=3$, $Q=2$}\\
\hline\hline
&& \multicolumn{2}{|c|}{$T=200$}& \multicolumn{2}{|c|}{$T=500$}& \multicolumn{2}{|c}{$T=1000$}\\
&& $n=100$ & $n=200$& $n=100$ & $n=200$& $n=100$ & $n=200$\\
\hline\hline
$MSE^X$ && 0.143 & 0.155 & 0.101 & 0.112 & 0.086 & 0.085 \\
$MSE^\chi$&$\kappa_T=0$ & 0.291 & 0.284 & 0.237 & 0.216 & 0.209 & 0.185 \\
$MSE^\chi$ & $\kappa_T=0.2$ & 0.262 & 0.261 & 0.197 & 0.199 & 0.179 & 0.163 \\
$MSE^\chi$ & $\kappa_T=0.4$ & 0.250 & 0.252 & 0.182 & 0.176 & 0.162 & 0.141 \\
\hline
$MAD^X$ && 0.265 & 0.261 & 0.227 & 0.228 & 0.210 & 0.205 \\
$MAD^\chi$&$\kappa_T=0$ & 0.412 & 0.402 & 0.370 & 0.348 & 0.346 & 0.322 \\
$MAD^\chi$ & $\kappa_T=0.2$ & 0.389 & 0.381 & 0.336 & 0.327 & 0.317 & 0.299 \\
$MAD^\chi$ & $\kappa_T=0.4$ & 0.378 & 0.372 & 0.321 & 0.307 & 0.300 & 0.275 \\
\hline
$MAX^X$ && 6.634 & 8.693 & 9.682 & 12.606 & 10.943 & 19.277 \\
$MAX^\chi$&$\kappa_T=0$ & 5.532 & 5.624 & 5.128 & 6.725 & 4.488 & 5.209 \\
$MAX^\chi$ & $\kappa_T=0.2$ & 5.688 & 5.800 & 5.561 & 6.204 & 4.813 & 4.675 \\
$MAX^\chi$ & $\kappa_T=0.4$ & 5.713 & 6.180 & 5.630 & 7.317 & 4.978 & 5.235 \\
\hline\hline
\\
\multicolumn{8}{c}{$q=2$, $Q=3$}\\
\hline\hline
&& \multicolumn{2}{|c|}{$T=200$}& \multicolumn{2}{|c|}{$T=500$}& \multicolumn{2}{|c}{$T=1000$}\\
&& $n=100$ & $n=200$& $n=100$ & $n=200$& $n=100$ & $n=200$\\
\hline\hline
$MSE^X$ & & 0.151 & 0.161 & 0.112 & 0.129 & 0.082 & 0.085 \\
$MSE^\chi$ & $\kappa_T=0$ & 0.356 & 0.323 & 0.290 & 0.277 & 0.247 & 0.224 \\
$MSE^\chi$ & $\kappa_T=0.2$ & 0.324 & 0.299 & 0.259 & 0.247 & 0.210 & 0.190 \\
$MSE^\chi$ & $\kappa_T=0.4$ & 0.311 & 0.287 & 0.248 & 0.231 & 0.193 & 0.173 \\
\hline
$MAD^X$ & & 0.279 & 0.271 & 0.241 & 0.248 & 0.212 & 0.213 \\
$MAD^\chi$ & $\kappa_T=0$ & 0.461 & 0.434 & 0.414 & 0.398 & 0.382 & 0.360 \\
$MAD^\chi$ & $\kappa_T=0.2$ & 0.437 & 0.415 & 0.387 & 0.372 & 0.350 & 0.327 \\
$MAD^\chi$ & $\kappa_T=0.4$ & 0.426 & 0.404 & 0.376 & 0.356 & 0.334 & 0.311 \\
\hline
$MAX^X$ & & 5.395 & 9.180 & 5.900 & 8.936 & 6.210 & 11.411 \\
$MAX^\chi$ & $\kappa_T=0$ & 4.833 & 5.066 & 5.009 & 5.654 & 5.111 & 5.571 \\
$MAX^\chi$ & $\kappa_T=0.2$ & 4.988 & 5.325 & 5.269 & 5.998 & 5.355 & 5.654 \\
$MAX^\chi$ & $\kappa_T=0.4$ & 5.058 & 5.748 & 5.660 & 6.208 & 5.599 & 5.413 \\
\hline\hline
\end{tabular}
\end{table}
\begin{figure}[t!]\caption{\small Simulation results. True (blue) and estimated (red) common components of levels, $\widehat{X}_{it,m}$, and of volatilties, $\exp(\widehat{\chi}_{it,m})$, when $n=200$, $T=1000$, $q=1$, $Q=1$, and $\kappa_T=0.2$. One series and one realisation.}\label{fig:sim1014}
\centering \smallskip\noindent
\setlength{\tabcolsep}{.01\textwidth}
\begin{tabular}{@{}cc}
\includegraphics[width=.5\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{X14.eps}&
\includegraphics[width=.5\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{chi14.eps} \\
\end{tabular}
\end{figure}
Finally, for $T=1000$, we estimated the model using the first 900 observations, then ran a recursive pseudo-out-of-sample forecasting exercise constructing one-step-ahead prediction intervals for the remaining $100$ observations (from $901$ to $1000$), as described in Section~\ref{eq:int}. The $\alpha/2$-upper and $\alpha/2$-lower bounds $\widehat{\mathcal U}_{i,\tau+1|\tau,m}(\alpha/2)$ and~$\widehat{\mathcal L}_{i,\tau+1|\tau,m}(\alpha/2)$ of prediction intervals with coverage probability $(1 - \alpha)$ are then computed for each series and replication and each out-of-sample observation. From the latter, we compute the observed coverage frequencies across all series and replications
\[
C(\alpha) := \frac 1 {\mathcal Mn100}\sum_{m=1}^\mathcal M\sum_{i=1}^n \sum_{\tau=900}^{999}
\mathbb I\Big(\widehat{\mathcal L}_{i,\tau+1|\tau,m}(\alpha/2)\le Y_{i,\tau+1,m}\! \le \widehat{\mathcal U}_{i,\tau+1|\tau,m}(\alpha/2)\Big)\nonumber
\]
and the proportions of coverage violations in the upper and lower tails,
\begin{align}
&V_{+}(\alpha/2) := \frac 1 {\mathcal Mn100}\sum_{m=1}^\mathcal M\sum_{i=1}^n \sum_{\tau=900}^{999} \mathbb I\Big(Y_{i,\tau+1,m}\! > \widehat{\mathcal U}_{i,\tau+1|\tau,m}(\alpha/2)\Big),\nonumber
\end{align}
and
\begin{align}
&V_{-}(\alpha/2) :=\frac 1 {\mathcal Mn100}\sum_{m=1}^\mathcal M\sum_{i=1}^n \sum_{\tau=900}^{999} \mathbb I\Big(Y_{i,\tau+1,m}\! < \widehat{\mathcal L}_{i,\tau+1|\tau,m}(\alpha/2)\Big),\nonumber
\end{align}
respectively. Results are shown in Table \ref{tab:covsim}. Overall performances look reasonably good---the larger $n$ and $T$, the better. We note that capping has a clear effect on the empirical coverage; too much capping seems to affect mostly the cases in which $\alpha=0.32$ and $0.2$. No capping at all works quite well in practice, despite the fact that theoretical results require $\kappa_T>0$. Moreover, the same comments apply to empirical coverage as for the choice of bandwidths, with the additional finding that higher values of $B_T$ yield more reliable prediction performances (see Appendix \ref{app:sim}).
\begin{table}[t!]
\centering
\caption{
\small Simulation results. Empirical coverage and frequencies of prediction bounds violations, averaged over all $n$ series and all $\mathcal M$ replications, for $T=1000$ and $\mathcal M=200$. Bandwidths values: $B_T=2$ and $M_T=10$ for~$T=200$; $B_T=2$ and $M_T=15$ for~$T=500$; $B_T=2$ and $M_T=20$ for~$T=1000$.}\label{tab:covsim}
\vskip .2cm
\footnotesize
\begin{tabular}{l l | ccccc | ccccc }
\multicolumn{12}{c}{$q=1$, $Q=1$}\\
\hline
\hline
&&\multicolumn{5}{c|}{$n=100$}&\multicolumn{5}{c}{$n=200$}\\
\hline
&&\multicolumn{5}{c|}{$\alpha$}&\multicolumn{5}{c}{$\alpha$}\\
&&0.32&0.2&0.1&0.05&0.01&0.32&0.2&0.1&0.05&0.01\\
\hline
$C(\alpha)$& $\kappa_T=0$&0.6409& 0.7637& 0.8667& 0.9312& 0.9869&0.6765& 0.7992& 0.9082& 0.9573& 0.9926 \\
$V_+(\alpha/2)$&&0.1810& 0.1195& 0.0674& 0.0342& 0.0057&0.1635& 0.1026& 0.0470& 0.0221& 0.0040\\
$V_-(\alpha/2)$&&0.1781& 0.1168& 0.0659& 0.0346& 0.0074&0.1601& 0.0983& 0.0449& 0.0206& 0.0035\\
\hline
$C(\alpha)$& $\kappa_T=0.2$&0.6691& 0.7769& 0.8685& 0.9197& 0.9681 &0.7201& 0.8285& 0.9226 &0.9628& 0.9934\\
$V_+(\alpha/2)$&&0.1636& 0.1124& 0.0682& 0.0403& 0.0153&0.1422& 0.0876& 0.0392& 0.0194& 0.0034\\
$V_-(\alpha/2)$&&0.1673& 0.1107& 0.0633& 0.0400& 0.0166& 0.1378& 0.0840& 0.0383& 0.0179& 0.0033\\
\hline
$C(\alpha)$& $\kappa_T=0.4$& 0.7072& 0.7987& 0.8799 &0.9238 &0.9699&0.7119& 0.7957& 0.8763& 0.9257& 0.9703\\
$V_+(\alpha/2)$&&0.1453& 0.1023& 0.0617& 0.0384& 0.0145&0.1429& 0.1007& 0.0600& 0.0360& 0.0150\\
$V_-(\alpha/2)$&&0.1475& 0.0990& 0.0584& 0.0378& 0.0156&0.1453& 0.1037& 0.0638& 0.0383& 0.0147\\
\hline
\hline
\\
\multicolumn{12}{c}{$q=3$, $Q=2$}\\
\hline
\hline
&&\multicolumn{5}{c|}{$n=100$}&\multicolumn{5}{c}{$n=200$}\\
\hline
&&\multicolumn{5}{c|}{$\alpha$}&\multicolumn{5}{c}{$\alpha$}\\
&&0.32&0.2&0.1&0.05&0.01&0.32&0.2&0.1&0.05&0.01\\
\hline
$C(\alpha)$& $\kappa_T=0$& 0.6350& 0.7523& 0.8543& 0.9131& 0.9688&0.6718& 0.7917& 0.8929& 0.9457& 0.9895\\
$V_+(\alpha/2)$&&0.1810& 0.1252& 0.0749& 0.0449& 0.0157&0.1619& 0.1037& 0.0530& 0.0277& 0.0053\\
$V_-(\alpha/2)$&&0.1840& 0.1225& 0.0708& 0.0420& 0.0155& 0.1664& 0.1047& 0.0542& 0.0266& 0.0053\\
\hline
$C(\alpha)$& $\kappa_T=0.2$&0.6767& 0.7775& 0.8665& 0.9206& 0.9701&0.7031& 0.8081& 0.8986& 0.9469& 0.9916\\
$V_+(\alpha/2)$&&0.1601& 0.1145& 0.0691& 0.0422& 0.0162&0.1458&0.0935& 0.0503& 0.0250& 0.0039\\
$V_-(\alpha/2)$&&0.1632& 0.1080& 0.0644& 0.0372& 0.0137& 0.1512& 0.0985& 0.0512& 0.0282& 0.0046\\
\hline
$C(\alpha)$& $\kappa_T=0.4$&0.7129& 0.7993& 0.8803& 0.9267& 0.9724&0.7565& 0.8447& 0.9222& 0.9610& 0.9923 \\
$V_+(\alpha/2)$&&0.1445& 0.1033& 0.0614& 0.0384& 0.0147&0.1209& 0.0780& 0.0382& 0.0196& 0.0043\\
$V_-(\alpha/2)$&&0.1426& 0.09740& 0.0583& 0.0349& 0.0129&0.1227& 0.0774& 0.0397& 0.0195& 0.0035\\
\hline
\hline
\\
\multicolumn{12}{c}{$q=2$, $Q=3$}\\
\hline
\hline
&&\multicolumn{5}{c|}{$n=100$}&\multicolumn{5}{c}{$n=200$}\\
\hline
&&\multicolumn{5}{c|}{$\alpha$}&\multicolumn{5}{c}{$\alpha$}\\
&&0.32&0.2&0.1&0.05&0.01&0.32&0.2&0.1&0.05&0.01\\
\hline
$C(\alpha)$& $\kappa_T=0$&0.6888& 0.8045& 0.8981& 0.9500& 0.9874&0.6391& 0.7563& 0.8623& 0.9215& 0.9784 \\
$V_+(\alpha/2)$&&0.1568& 0.0995& 0.0527& 0.0258& 0.0065&0.1786& 0.1212& 0.0678& 0.0387& 0.0107\\
$V_-(\alpha/2)$&& 0.1544& 0.0960& 0.0492& 0.0242& 0.0061&0.1824& 0.1226& 0.0700& 0.0399& 0.0110\\
\hline
$C(\alpha)$& $\kappa_T=0.2$&0.7335&0.8290& 0.9105& 0.9539& 0.9890&0.6770& 0.7814& 0.8752& 0.9277& 0.9791\\
$V_+(\alpha/2)$&&0.1345& 0.0863& 0.0459& 0.0239& 0.0056&0.1602& 0.1077& 0.0613& 0.0355& 0.0103\\
$V_-(\alpha/2)$&&0.1320& 0.0847& 0.0436& 0.0222& 0.0054&0.1629& 0.1110& 0.0636& 0.0368& 0.0106\\
\hline
$C(\alpha)$& $\kappa_T=0.4$&0.7733&0.8539& 0.9213& 0.9595& 0.9898&0.7167& 0.8066& 0.8879& 0.9352& 0.9809 \\
$V_+(\alpha/2)$&&0.1154& 0.0747& 0.0408& 0.0207& 0.0055&0.1395& 0.0953& 0.0551& 0.0326& 0.0093\\
$V_-(\alpha/2)$&&0.1113& 0.0714& 0.0379& 0.0198& 0.0047&0.1439& 0.0981& 0.0570& 0.0322& 0.0099\\
\hline
\hline
\end{tabular}
\end{table}
\setcounter{equation}{0}
\section{Interval prediction for S\&P100 returns }\label{sec:emp}
In this section, we apply our methodology to a panel of $n=90$ daily returns of stocks from the Standard \& Poor's 100 Index. Data are observed from January 4, 2000 through September 30, 2013, for a total of~$T=3456$ observations. We run a pseudo-out-of-sample forecasting exercise by estimating the model using data over the period $t=1,\ldots, \tau$, with $\tau=(T-M),\ldots, (T-1)$ and $M=~\!1948$, corresponding to an evaluation period running from January 3, 2006 through September 27, 2013. For each value of $\tau$, we estimate the $n=90$ one-step-ahead prediction intervals as defined in \eqref{eq:Itt1hat}. The data cover the following sectors (in parentheses, the number of series in each sector): Consumer Discretionary~(11), Consumer Staples (10), Energy (12), Financials (13), Health Care (11), Industrials~(14), Information Technology~(12), Materials (3), Telecommunications Services (2), Utilities (2) (see Appendix \ref{app:data} for the names of individual stocks).
Although we should, in principle, fully re-estimate the whole model at each of the $M$ iterations, some quantities were kept fixed throughout the exercise. In particular,
when applied to the full~$n\times~\!T$ panel, the \citet{hallinliska07} criterion returns $\widehat q=3$ common factors for the level panel and~$\widehat Q=~\!2$ common factors for log-volatility panel: those values are used in all subsequent analyzes. We also choose the bandwidths by minimizing, over a grid of possible bandwidth values, the mean-squared errors
\begin{align}
\frac 1 {nT} \sum_{i=1}^n \sum_{t=1}^T(Y_{it}-\widehat X_{it|t-1})^2 \quad\text{and}\quad \frac 1 {nT} \sum_{i=1}^n \sum_{t=1}^T(\widehat h_{it}-\widehat \chi_{it|t-1})^2,\nonumber
\end{align}
respectively, leading to possibly distinct bandwidths for $\widehat X_{it|t-1}$ and $\widehat \chi_{it|t-1}$. More precisely, we first determine $B_T$ and then determine $M_T$ using the chosen $B_T$ to compute $\widehat h_{it}$. As a result we throughout use $B_T=2$ and~$M_T=17$. The VAR orders and the orders of their truncated inverse MA representations needed to compute impulse responses are set as follows:
\begin{inparaenum}[(i)]
\item $\text{deg}[\mbf A_n(L)]=1$, with inverse MA truncated at lag $\bar k_1=20$;
\item $\text{deg}[\mbf C_n(L)]=1$, with inverse MA truncated at lag $\bar k_2=20$;
\item $\text{deg}[\mbf M_n(L)]=5$, with inverse MA truncated at lag~$\bar k_1^*=~\!100$;
\item $\text{deg}[\mbf P_n(L)]=1$, with inverse MA truncated at lag $\bar k_2^*=100$.
\end{inparaenum}
The estimation of the GDFM is based on 10 cross-sectional permutations, as explained at the end of Section~\ref{sec:est_summary}. Finally, regarding the choice of the capping constant $\kappa_T$, we choose $\kappa_T\in\{0,\, 0.1,\, 0.25,\, 0.5\}$ irrespective of~$i$; note that, with reference to Assumption (R), we have $\log^{-1} T= 0.12$. Also note that, on the average across the~$M$ iterations, 6\%, out of the total $n\tau$ observations, are capped when $\kappa_T=0.1$, 14\% when~$\kappa_T=~\!0.25$, and 27\% when $\kappa_T=0.5$.
For any given sample size $\tau$, we compute the quantiles of $\widehat{w}_i$ using $(\widehat w_{i,\tau-\ell+1},\ldots,\widehat w_{i,\tau})$, where we\linebreak set~$\ell\in\{126, 252, 504, \tau\}$, hence using either the past six months, one year, or two years of available data, or using all available past observations. Denoting by $\widehat {\mbf w}^{(\ell)}_{i}$ the vector of the most recent $\ell$ observations (so that $\widehat {\mbf w}^{(\tau)}_{i}$ coincides with $\widehat {\mbf w}_{i}$), for levels $\alpha\in\{0.32,0.2,0.1,0.05,0.01\}$ and window sizes $\ell$, and for $\tau=(T-M),\ldots, (T-1)$, we obtain the estimates
\begin{align}
&\widehat{\mathcal U}_{i,\tau+1|\tau}^{(\ell)}(\alpha):=\widehat Y_{i,\tau+1|\tau} + \widehat s_{i,\tau+1|\tau}\, \widehat w^{(\ell)}_{i(\lceil \ell(1-\alpha)\rceil)},&\quad \quad \quad
\widehat{\mathcal L}_{i,\tau+1|\tau}^{(\ell)}(\alpha):=\widehat Y_{i,\tau+1|\tau} + \widehat s_{i,\tau+1|\tau}\, \widehat w^{(\ell)}_{i(\lceil \ell\alpha\rceil)},\nonumber\\
&\widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)}(\alpha):=\big[\widehat{\mathcal L}^{(\ell)}_{i,\tau+1|\tau}(\alpha^-),\widehat{\mathcal U}^{(\ell)}_{i,\tau+1|\tau}(\alpha^+)\big], \quad\text{ and}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!&\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha):=\mathbb I\big(Y_{i,\tau+1}\in \widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)}(\alpha)\big).\ \ \ \ \,\nonumber
\end{align}
\subsection{Coverage performance: qualitative analysis}
For each of the $n=90$ series considered we compute the coverage frequency
\[
C_i^{(\ell)}(\alpha) := \frac 1 {M} \sum_{\tau=T-M}^{T-1} \widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha)=\frac 1M\sum_{\tau=T-M}^{T-1} \mathbb I\Big(\widehat{\mathcal L}_{i,\tau+1|\tau}^{(\ell)}(\alpha^-)\le Y_{i,\tau+1}\! \le \widehat{\mathcal U}_{i,\tau+1|\tau}^{(\ell)}(\alpha^+)\Big),
\]
the proportions
\begin{align}
V^{(\ell)}_{i,+}(\alpha^+) \! := \frac 1 {M}\!\! \sum_{\tau=T-M}^{T-1}\! \! \!\!\mathbb I\Big(Y_{i,\tau+1}\! > \widehat{\mathcal U}_{i,\tau+1|\tau}^{(\ell)}(\alpha^+)\Big)
\text{ and }
V^{(\ell)}_{i,-}(\alpha^-)\! := \frac 1 {M}\!\! \sum_{\tau=T-M}^{T-1} \! \! \!\!\mathbb I\Big(Y_{i,\tau+1}\! < \widehat{\mathcal L}_{i,\tau+1|\tau}^{(\ell)}(\alpha^-)\Big)\nonumber
\end{align}
of coverage violations in the upper and lower tails, and the average interval length
\[
L_{i}^{(\ell)}(\alpha):=\frac 1{M} \sum_{\tau=T-M}^{T-1} \left( \widehat{\mathcal U}_{i,\tau+1|\tau}^{(\ell)}(\alpha)-\widehat{\mathcal L}_{i,\tau+1|\tau}^{(\ell)}(\alpha)\right).
\]
Table \ref{tab:viol} reports, for $\alpha^+=\alpha^-=\alpha/2$ with $\alpha\in\{0.32,0.2,0.1,0.05,0.01\}$ (corresponding to coverage levels 68\%, 80\%, 90\%, 95\% and 99\%) and $\kappa_T\in\{0,0.1,0.25,0.5\}$, the cross-sectional average $C^{(\ell)}(\alpha)$ of the empirical coverage frequencies $C_i^{(\ell)}(\alpha)$, the cross-sectional averages $V^{(\ell)}_+(\alpha/2)$ and $V^{(\ell)}_-(\alpha/2)$ of the proportions of coverage violations $V^{(\ell)}_{i,+}(\alpha^+)$ and $V^{(\ell)}_{i,-}(\alpha^-)$, and the cross-sectional average ${L}^{(\ell)}(\alpha)$ of the average interval lengths~$L_{i}^{(\ell)}(\alpha)$.
\begin{table}[t!]
\centering
\caption{\small Standard \& Poor's 100 Index data ($n=90$ daily returns). Empirical coverage, frequency of prediction bounds violations, and average length of prediction intervals for GDFM, averaged over the cross-section.}\label{tab:viol}
\vskip .2cm
\footnotesize
\begin{tabular}{l | ccccc | ccccc}
\hline
\hline
&\multicolumn{5}{c|}{$\kappa_T=0$}&\multicolumn{5}{c}{$\kappa_T=0.1$}\\
\hline
&\multicolumn{5}{c|}{$\alpha$}&\multicolumn{5}{c}{$\alpha$}\\
&0.32&0.2&0.1&0.05&0.01 &0.32&0.2&0.1&0.05&0.01 \\
\hline
$C^{(126)}(\alpha)$ &0.6709 & 0.7894 & 0.8887 & 0.9400 & 0.9812 &0.6874 & 0.7985& 0.8931& 0.9416& 0.9813 \\
$V^{(126)}_+(\alpha/2)$&0.1641 & 0.1048 & 0.0552 & 0.0299 & 0.0094&0.1559& 0.1002& 0.0533& 0.0291& 0.0095 \\
$V^{(126)}_-(\alpha/2)$&0.1650 & 0.1058 & 0.0561 & 0.0301 & 0.0094&0.1566& 0.1013& 0.0536& 0.0292& 0.0091 \\
${L}^{(126)}(\alpha)$&3.3934 & 4.5156 & 6.1305 & 7.7681 & 12.4174&3.4726 & 4.5726 & 6.1553 & 7.7698 & 12.3130\\
\hline
$C^{(252)}(\alpha)$ &0.6708 & 0.7903 & 0.8902 & 0.9415 & 0.9848 &0.6882& 0.7999& 0.8940& 0.9424& 0.9846 \\
$V^{(252)}_+(\alpha/2)$&0.1647 & 0.1044 & 0.0544 & 0.0289 & 0.0077&0.1560& 0.0998& 0.0526& 0.0287& 0.0078 \\
$V^{(252)}_-(\alpha/2)$&0.1644 & 0.1053 & 0.0554 & 0.0296 & 0.0075&0.1558& 0.1003& 0.0534& 0.0289& 0.0076 \\
${L}^{(252)}(\alpha)$&3.3621 & 4.4794 & 6.0949 & 7.7240 & 12.5008&3.4351 & 4.5290 & 6.1078 & 7.7074 &12.3767\\
\hline
$C^{(504)}(\alpha)$ &0.6711 & 0.7895 & 0.8895 & 0.9412 & 0.9846 &0.6886& 0.7995& 0.8929& 0.9419& 0.9843 \\
$V^{(504)}_+(\alpha/2)$&0.1651 & 0.1057 & 0.0551 & 0.0290 & 0.0078&0.1561& 0.1005& 0.0536& 0.0288& 0.0081 \\
$V^{(504)}_-(\alpha/2)$&0.1638 & 0.1047 & 0.0554 & 0.0298 & 0.0076&0.1553& 0.1000& 0.0535& 0.0292& 0.0076 \\
${L}^{(504)}(\alpha)$&3.3034 & 4.4179 & 6.0266 & 7.6643 & 12.1190&3.3786 & 4.4708 & 6.0439 & 7.6539 & 12.0462\\
\hline
$C^{(\tau)}(\alpha)$& 0.7010 & 0.8142 & 0.9049 & 0.9506 & 0.9881&0.7187& 0.8244& 0.9096& 0.9523 & 0.9881 \\
$V^{(\tau)}_+(\alpha/2)$&0.1516 & 0.0933 & 0.0474 & 0.0247 & 0.0061&0.1424& 0.0879& 0.0452& 0.0237& 0.0062 \\
$V^{(\tau)}_-(\alpha/2)$&0.1474 & 0.0925 & 0.0477 & 0.0248 & 0.0058&0.1389& 0.0877& 0.0452& 0.0239& 0.0057 \\
${L}^{(\tau)}(\alpha)$&3.4523 & 4.6632 & 6.4305 & 8.2802 & 13.3895&3.5562 & 4.7560 & 6.5201 & 8.3747 & 13.5115\\
\hline
\hline
&\multicolumn{5}{c|}{$\kappa_T=0.25$}&\multicolumn{5}{c}{$\kappa_T=0.5$}\\
\hline
&\multicolumn{5}{c|}{$\alpha$}&\multicolumn{5}{c}{$\alpha$}\\
&0.32&0.2&0.1&0.05&0.01 &0.32&0.2&0.1&0.05&0.01 \\
\hline
$C^{(126)}(\alpha)$ &0.7126 & 0.8141 & 0.8997 & 0.9452 & 0.9821&0.7552& 0.8391& 0.9119& 0.9507& 0.9836 \\
$V^{(126)}_+(\alpha/2)$&0.1435 & 0.0926 & 0.0500 & 0.0274 & 0.0091&0.1222& 0.0800& 0.0436& 0.0243& 0.0082 \\
$V^{(126)}_-(\alpha/2)$&0.1439 & 0.0932 & 0.0504 & 0.0274 & 0.0088&0.1226& 0.0809& 0.0446& 0.0251& 0.0081 \\
${L}^{(126)}(\alpha)$&3.6203 & 4.6949 & 6.2419 & 7.8371 & 12.3547&3.9076& 4.9426 & 6.4443 & 8.0189 & 12.5330\\
\hline
$C^{(252)}(\alpha)$ & 0.7138 & 0.8143 & 0.9009 & 0.9452 & 0.9851&0.7556 & 0.8398 & 0.9127& 0.9512& 0.9866\\
$V^{(252)}_+(\alpha/2)$&0.1433 & 0.0928& 0.0491& 0.0271 & 0.0077&0.1224& 0.0798 & 0.0432& 0.0242& 0.0069 \\
$V^{(252)}_-(\alpha/2)$&0.1428 & 0.0929 & 0.0500 & 0.0277 & 0.0072&0.1219& 0.0804& 0.0440& 0.0246& 0.0065 \\
${L}^{(252)}(\alpha)$&3.5796& 4.6500& 6.1923& 7.7700& 12.4405&3.8737& 4.9105 & 6.4108 & 7.9706 & 12.7086\\
\hline
$C^{(504)}(\alpha)$ &0.7149 & 0.8150 & 0.9002 & 0.9449 & 0.9846&0.7588& 0.8422 & 0.9132& 0.9514& 0.9861 \\
$V^{(504)}_+(\alpha/2)$& 0.1430 & 0.0927 & 0.0495 & 0.0274 & 0.0080&0.1204 & 0.0782& 0.0426& 0.0236& 0.0072 \\
$V^{(504)}_-(\alpha/2)$&0.1420 & 0.0923 & 0.0502 & 0.0277 & 0.0074&0.1208& 0.0795& 0.0442 & 0.0250& 0.0067 \\
${L}^{(504)}(\alpha)$&3.5336 & 4.6035 & 6.1513 & 7.7434 & 12.1898&3.8584 & 4.9070 & 6.4371 & 8.0357 & 12.6302\\
\hline
$C^{(\tau)}(\alpha)$&0.7430 & 0.8387 & 0.9162 & 0.9551 & 0.9886& 0.7824& 0.8633 & 0.9283 & 0.9613 & 0.9900 \\
$V^{(\tau)}_+(\alpha/2)$&0.1301 & 0.0808 & 0.0415 & 0.0221 & 0.0061&0.1091& 0.0680& 0.0351 & 0.0189& 0.0053 \\
$V^{(\tau)}_-(\alpha/2)$&0.1269 & 0.0805 & 0.0422 & 0.0228 & 0.0054&0.1085& 0.0687 & 0.0366 & 0.0198 & 0.0047\\
${L}^{(\tau)}(\alpha)$&3.7420 & 4.9317 & 6.6982 & 8.5677 & 13.7991&4.1045 & 5.2913& 7.0734& 8.9901& 14.4295\\
\hline
\hline
\end{tabular}
\end{table}
Inspection of the table reveals that $C^{(\ell)}(\alpha) \simeq(1-\alpha)$ and $V_+^{(\ell)}(\alpha/2) \simeq V_-^{(\ell)}(\alpha/2) \simeq \alpha/2$, which is a qualitative confirmation of the validity of our methodology (see Section~\ref{backtestSec} for more formal validation). Three remarks emerge from these results.
First, regarding the sensitivity of our procedure to capping, lower values of $\kappa_T$, in general, provide better results when $\alpha$ is higher, while larger values of~$\kappa_T$ provide better results for lower values of $\alpha$; in all cases, $\kappa_T=0.5$ yields a mostly conservative coverage frequency higher than $(1-\alpha)$. In particular, note that the choice of $\kappa_T=0$ (no capping at all), although ruled out by Assumption~(R), still provides very good results.
Second, setting $\ell=\tau$, that is, considering the entire past history to compute quantiles apparently is not the best strategy, and shorter horizons $\ell$ seem preferable. This finding is possibly related to some time variation in the distribution of the innovations of log-volatilities at horizons longer than one year. Third, for any given~$\alpha$, shorter intervals are obtained when setting $\ell=252$ or $504$ regardless of the choice of $\kappa_T$. Overall, choosing~$\kappa_T=0.1$ and~$\ell=252$ or $504$ works best for $\alpha=0.32$ and $0.2$, while~$\kappa_T=0.25$ and $\ell=126$ or $252$ works best for~$\alpha=0.1, 0.05$, and $0.01$.
\begin{figure}[t!]\caption{\small One-step-ahead 90\% conditional prediction intervals (in red; $\ell=252$): America International Group (AIG), Bank of America (BAC), Citigroup (C), Goldman Sachs (GS), JPMorgan Chase (JPM), Morgan Stanley (MS).}\label{fig:stocks}
\centering \smallskip\noindent
\setlength{\tabcolsep}{.01\textwidth}
\begin{tabular}{@{}cc}
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{AIG.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{BAC.eps} \\
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{C.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{GS.eps}\\
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{JPM.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{MS.eps}\\
\end{tabular}
\end{figure}
In Figures \ref{fig:stocks} and \ref{fig:stocks2}, we set $\kappa_T=0.25$ and $\ell=252$ and we show (in grey) $Y_{i,\tau+1}$ for some selected individual stocks, together with (in red) the estimated upper and lower bounds of the 90\% one-step-ahead prediction interval, i.e. $\widehat{\mathcal U}^{(252)}_{i,\tau+1|\tau}(0.05)$ and $\widehat{\mathcal L}^{(252)}_{i,\tau+1|\tau}(0.05)$, respectively. Figure \ref{fig:stocks} shows results for six of the most volatiles stocks in our dataset, all belonging to the financial sector: America International Group (AIG), Bank of America (BAC), Citigroup (C), Goldman Sachs (GS), JPMorgan Chase (JPM), Morgan Stanley (MS). Figure \ref{fig:stocks2} provides the same results for eight relevant non-financial stocks: Apple (AAPL), Microsoft (MSFT), Amazon (AMZN), Wallgreens (WAG), Exxon Mobil (XOM), Johnson \& Johnson (JNJ), Boeing (BA), General Electric (GE). Volatilities, in those series, which were the most seriously affected by the great financial crisis, are notoriously hard to predict.
\begin{figure}[t!]\caption{\small One-step-ahead 90\% conditional prediction intervals (in red; $\ell=252$): Apple (AAPL), Microsoft (MSFT), Amazon (AMZN), Wallgreens (WAG), Exxon Mobil (XOM), Johnson \& Johnson (JNJ), Boeing (BA), General Electric (GE).} \label{fig:stocks2}
\centering \smallskip\noindent
\setlength{\tabcolsep}{.01\textwidth}
\begin{tabular}{@{}cc}
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{AAPL.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{MSFT.eps}\\
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{AMZN.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{WAG.eps} \\
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{XOM.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{JNJ.eps}\\
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{BA.eps}&
\includegraphics[width=.45\textwidth,trim=3.5cm 1.6cm 2cm 0cm,clip]{GE.eps}\\
\end{tabular}
\end{figure}
\subsection{Coverage: comparison with GARCH}
The novelty of our prediction intervals is that they are exploiting the information contained in the available cross-section of $n=90$ stocks. This is in sharp contrast with the usual GARCH approach, which is strongly univariate, and disregards cross-sectional information by analyzing the $n$ series one by one. Moreover, estimating 90 univariate GARCH models requires much more computing time than estimating our model. GARCH nevertheless constitute the more common practice in this context, and serves as a natural benchmark.
We therefore compare our prediction intervals with those obtained by fitting, {via quasi-maximum likelihood,} univariate GARCH(1,1) models to all series in our panel. Specifically, for each series $i$, we estimate the model
\begin{align}
Y_{it} &= \mathrm E[y_{it}]+ \sigma_{it}\epsilon_{it}, \quad \epsilon_{it}\stackrel{iid}{\sim} (0,1),\qquad t=1,\ldots, \tau,\nonumber\\
\sigma_{it}^2&=\omega_i +\gamma_i Y_{it-1}^2+ \beta_i\sigma_{it-1}^2, \quad \omega_i>0,\ \gamma_i,\beta_i\ge 0,\ \gamma_i+\beta_i<1.\nonumber
\end{align}
For given $\tau=(T-M),\ldots, (T-1)$, we obtain estimated parameters $\widehat{\omega}_i, \widehat{\gamma}_i$ and $\widehat{\beta}_i$, from which we compute the estimated volatilities $\widehat{\sigma}_{it}^2$ and the innovation values~$\widehat{\epsilon}_{i t}=Y_{it}/\widehat{\sigma}_{it}$, $t=1,\ldots, \tau$.
Innovation quantiles are computed from $(\widehat \epsilon_{i,\tau-\ell+1},\ldots,\widehat \epsilon_{i,\tau})$, where as before we set $\ell\in\{126,\, 252,\, 504,\, \tau\}$. Then, for any given level $\alpha$ and window size $\ell$, and for $\tau=(T-M),\ldots, (T-1)$, given the one-step-ahead volatility pre- \linebreak
dictor~$\widehat{\sigma}_{i,\tau+1|\tau}^2 = \widehat{\omega}_i+\widehat{\gamma}_i Y_{i,\tau}^2+\widehat{\beta}_i \widehat{\sigma}_{i,\tau}^2$, we compute the the upper and lower confidence bounds
\[
\widehat{\mathcal U}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha):=\bar Y_{i} + \widehat \sigma_{i,\tau+1|\tau}\, \widehat \epsilon^{(\ell)}_{i(\lceil \ell(1-\alpha)\rceil)} \quad \text{and}\quad
\widehat{\mathcal L}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha):=\bar Y_{i} + \widehat \sigma_{i,\tau+1|\tau}\, \widehat \epsilon^{(\ell)}_{i(\lceil \ell\alpha\rceil)},\nonumber
\]
yielding the one-step-ahead prediction intervals
\[
\widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha):=\big[\widehat{\mathcal L}^{(\ell)\text{\tiny GARCH}}_{i,\tau+1|\tau}(\alpha/2),\widehat{\mathcal U}^{(\ell)\text{\tiny GARCH}}_{i,\tau+1|\tau}(\alpha/2)\big]
\]
and the indicators of correct interval prediction $\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha):=\mathbb I(Y_{i,\tau+1}\in \widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha))$. Based on these quantities, we then compute, for $\alpha\in\{0.32,\, 0.2,\, 0.1,\, 0.05,\, 0.01\}$,
the empirical coverage frequency, denoted as $C_i^{(\ell)\text{\tiny{GARCH}}}(\alpha)$, the
proportions of coverage violations in the upper and lower tail, denoted as $V_{i,+}^{(\ell)\text{\tiny{GARCH}}}(\alpha/2)$ and~$V_{i,-}^{(\ell)\text{\tiny{GARCH}}}(\alpha/2)$, respectively, and the average interval length, denoted as ${L}_i^{(\ell)\text{\tiny{GARCH}}}(\alpha)$.
Averages of these quantities over the $n$ series under study are shown in Table \ref{tab:garch}. Inspection of this table reveals that the GDFM performances are slightly better than the GARCH ones in terms of coverage frequencies, based on similar interval lengths. This, however, is mainly a descriptive and, due to cross-sectional dependence, somewhat misleading assessment, which ideally should be reinforced into a more formal testing analysis.
\begin{table}[t!]
\centering
\caption{\small Standard \& Poor's 100 Index data ($n=90$ daily returns). Empirical coverage, frequency of prediction bounds violations, and average length of prediction intervals for GARCH, averaged over the cross-section. }\label{tab:garch}
\vskip .2cm
\footnotesize
\begin{tabular}{l | ccccc }
\hline
\hline
&\multicolumn{5}{c}{$\alpha$}\\
&0.32&0.2&0.1&0.05&0.01 \\
\hline
$C^{(126)\text{\tiny{GARCH}}}(\alpha)$&0.6755 & 0.7947 & 0.8933 & 0.9429 & 0.9834\\
$V_{+}^{(126)\text{\tiny{GARCH}}}(\alpha/2)$&0.1576 & 0.0991 & 0.0507 & 0.0267 & 0.0076\\
$V_{-}^{(126)\text{\tiny{GARCH}}}(\alpha/2)$&0.1669 & 0.1062& 0.0560 & 0.0304 & 0.0090\\
$L^{(126)\text{\tiny{GARCH}}}(\alpha)$&3.4401 & 4.5562 & 6.1207 & 7.7282 & 12.3986\\
\hline
$C^{(252)\text{\tiny{GARCH}}}(\alpha)$&0.6786 & 0.7981 & 0.8968 & 0.9460 & 0.9871\\
$V_{+}^{(252)\text{\tiny{GARCH}}}(\alpha/2)$&0.1567 & 0.0978 & 0.0491 & 0.0255 & 0.0060\\
$V_{-}^{(252)\text{\tiny{GARCH}}}(\alpha/2)$& 0.1647 & 0.1041 & 0.0541 & 0.0285 & 0.0069\\
$L^{(252)\text{\tiny{GARCH}}}(\alpha)$&3.4142 & 4.5235 & 6.0755 & 7.6329 & 12.2536\\
\hline
$C^{(504)\text{\tiny{GARCH}}}(\alpha)$&0.6807 & 0.7994 & 0.8983 & 0.9479 & 0.9878\\
$V_{+}^{(504)\text{\tiny{GARCH}}}(\alpha/2)$&0.1560 & 0.0975 & 0.0488 & 0.0248 & 0.0056\\
$V_{-}^{(504)\text{\tiny{GARCH}}}(\alpha/2)$&0.1633 & 0.1031 & 0.0529 & 0.0274 & 0.0066\\
$L^{(504)\text{\tiny{GARCH}}}(\alpha)$&3.3822 & 4.4801 & 6.0220 & 7.5581 & 11.7469\\
\hline
$C^{(\tau)\text{\tiny{GARCH}}}(\alpha)$&0.6920 & 0.8077 & 0.9036 & 0.9510 & 0.9897\\
$V_{+}^{(\tau)\text{\tiny{GARCH}}}(\alpha/2)$&0.1520 & 0.0935 & 0.0458 & 0.0228 & 0.0048 \\
$V_{-}^{(\tau)\text{\tiny{GARCH}}}(\alpha/2)$& 0.1560 & 0.0988 & 0.0505 & 0.0262 & 0.0055\\
$L^{(\tau)\text{\tiny{GARCH}}}(\alpha)$&3.4139 & 4.4942 & 6.0156 & 7.5369 & 11.6268\\
\hline
\hline
\end{tabular}
\end{table}
A formal comparison between the GDFM and GARCH(1,1) coverage performances should take into account the fact that the coverage results of the two methods, for given $i$ and $\tau$, are not independent. The situation is quite similar to that of comparing paired proportions, where tests are to be carried out on the basis of the traditional \cite{McN47} test. For given $\alpha$ and $\ell$, consider, for all $i$, the events (discordant GDFM and GARCH coverage results)
\begin{align}
\mathcal A_{i,\tau+1|\tau}^{(\ell)}(\alpha)&:=\left\{Y_{i,\tau+1}\in \widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)}(\alpha) \cap Y_{i,\tau+1}\notin \widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha)\right\}\nonumber\\
\mathcal B_{i,\tau+1|\tau}^{(\ell)}(\alpha)&:=\left\{Y_{i,\tau+1}\notin \widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)}(\alpha) \cap Y_{i,\tau+1}\in \widehat{\mathcal I}_{i,\tau+1|\tau}^{(\ell)\text{\tiny GARCH}}(\alpha)\right\},\nonumber
\end{align}
and define
\[
n_{12i}^{(\ell)}(\alpha) := \sum_{\tau=T-M}^{T-1} \mathbb I\left(\mathcal A_{i,\tau+1|\tau}^{(\ell)}(\alpha) \right) \quad\text{ and }\quad
n_{21i}^{(\ell)}(\alpha) := \sum_{\tau=T-M}^{T-1} \mathbb I\left(\mathcal B_{i,\tau+1|\tau}^{(\ell)}(\alpha)\right).
\]
Consider the null hypothesis under which the indicators of a successful interval prediction in both methods are i.i.d.~Bernoulli, with identical (but otherwise unspecified) coverage probabilities.
The McNemar test of that hypothesis is conditioning on the sum $n_{\text{disc},i}^{(\ell)}(\alpha):=n_{12i}^{(\ell)}(\alpha) + n_{21i}^{(\ell)}(\alpha)$ of discordant coverage results: concordant results indeed carry no information on a difference between coverage probabilities. Conditional on $n_{\text{disc},i}^{(\ell)}(\alpha)$, the null distribution of $n_{12i}^{(\ell)}(\alpha)$ is binomial Bin$(n_{\text{disc},i}^{(\ell)}(\alpha),\, 0.5)$. At probability level~$\delta$, the test rejects in favour of a better GDFM coverage for ``large values'' of $n_{12i}^{(\ell)}(\alpha)$, in favour of a better GARCH coverage for ``small values'' of the same (equivalently, ``large values'' of~$n^{(\ell)}_{21i}$), with critical values the $(1-\delta)$ and $\delta$ binomial quantiles, respectively.
Table \ref{tab:comparetest1} reports the McNemar empirical rejection frequencies (over the $n=90$ series)---in favour of a better GDFM coverage in the left-hand panel, in favour of a better GARCH coverage in the right-hand one. We consider the cases in which $\alpha=0.1$ or $0.05$, $\ell=126$ or $252$, $\kappa_T=0.25$ (for the GDFM); testing was performed at significance levels $\delta=0.1$, $0.05$, and $0.01$. Irrespective of $\ell$ and $\alpha$, the GDFM approach appears to outperform, quite consistently and significantly, the GARCH one.
\begin{table}[t!]
\centering
\caption{\small Standard \& Poor's 100 Index data ($n=90$ daily returns). Proportions of McNemar rejections in favour of a better GDFM coverage (left-hand panel), in favour of a better GARCH coverage (right-hand panel).
}\label{tab:comparetest1}
\vskip .2cm
\footnotesize
\begin{tabular}{l | ccc | ccc }
\hline
\hline
&\multicolumn{3}{|c}{better GDFM coverage}&\multicolumn{3}{|c}{better GARCH coverage}\\
\hline
$\alpha=0.1$ & $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$ &$\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$\\
\hline
$\ell=126$&0.6000 & 0.5444 & 0.3667 & 0.1000& 0.0778 & 0.0667\\
$\ell=252$&0.5333 & 0.4556 & 0.2889 & 0.1333& 0.1222 & 0.0889\\
\hline
\hline
$\alpha=0.05$ & $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$ &$\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$\\
\hline
$\ell=126$&0.4111 & 0.2778 & 0.1556 &0.1111 & 0.0778 & 0.0556\\
$\ell=252$&0.2556 & 0.2000 & 0.0556&0.1778 & 0.1222 & 0.1111\\
\hline
\hline
\end{tabular}
\end{table}
\subsection{Coverage: backtesting}\label{backtestSec}
As explained in Section~\ref{eq:int}, a formal assessment of the validity of our approach can be based on the backtesting procedure proposed by \citet{christoffersen1998}. The idea consists in testing the null hypothesis \eqref{eq:null} under which the $\widehat{\mathcal H}_{i,\tau+1|\tau}(\alpha)$'s (the indicators of a successful interval prediction) are i.i.d.~Bernoulli$ (1-\alpha)$. Depending on the objectives, several alternatives can be considered. One can be interested (Section \ref{521}) in the validity of interval prediction or the sharpness of the nominal coverage level. Else, one may consider (Section \ref{522}) alternatives of serial dependence. Or, those two issues can be combined (Section \ref{523}) by merging the corresponding alternatives.
Irrespective of the alternative, however, it should be insisted that all those tests---one for each cross-sectional item---are intrinsically univariate. When simultaneously performing several or all of them, one should be extremely cautious with the interpretation of the results. The tables we are providing below are reporting empirical rejection frequencies (over the $n=90$ series). Those $n$ tests, however, are not functionally interrelated (as they would be if the prediction intervals were based on the quantiles of common shocks only); hence, they are not about testing the validity of {\it joint prediction intervals} with global asymptotic coverage level $(1-\alpha)$. Neither are they mildly interrelated (as they would be if the prediction intervals were exclusively based on idiosyncratic quantiles), providing joint prediction intervals with global asymptotic coverage level of the order of $(1-\alpha)^n$. High rejection frequencies across the $n$ series thus do not imply bad forecasting properties, but can result from complex cross-sectional dependencies. A standard attitude would consist in adopting a Bonferroni or a \v Sid\' ak correction; for $n=90$, and for a global testing level of $1\%$, this would lead to implementing the~$n=90$ individual tests at an overly conservative level $\delta \approx 0.0001 = 10^{-4}$---a level at which none of the null hypotheses under study is rejected.
All tests below are performed for $\kappa_T=0.25$, $\alpha=0.1$ or $0.05$, $\ell=126$ or $252$; testing significance levels are~$\delta=0.1$, $0.05$, and $0.01$.
\subsubsection{Testing for valid or sharp conditional coverage probabilities}\label{521}
If we are interested in the validity of interval prediction, the relevant testing problems are (one-sided)
\begin{equation}\label{eq:testUConesided}
H_{0i}: \mathrm E[\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha)] \geq (1-\alpha)\quad\text{ versus }\quad H_{1i}: \mathrm E[\widehat{\mathcal H}^{(\ell)}_{i,\tau+1|\tau}(\alpha)]< (1-\alpha).
\end{equation}
If instead we are interested in testing whether $(1-\alpha)$, as a nominal confidence level, is sharp, the testing problems are (still one-sided)
\begin{equation}\label{eq:testUConesided2}
H_{0i}: \mathrm E[\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha)] \leq (1-\alpha)\quad\text{ versus }\quad H_{1i}: \mathrm E[\widehat{\mathcal H}^{(\ell)}_{i,\tau+1|\tau}(\alpha)]> (1-\alpha).
\end{equation}
Both testing problems \eqref{eq:testUConesided} and \eqref{eq:testUConesided2},
admit a level-$\delta$ uniformly most powerful solution, rejecting $H_{0i}$ whenever the test statistic
$$
n_{1i}^{(\ell)}(\alpha) :=\sum_{\tau=T-M}^{T-1}\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha)
$$
falls below the binomial Bin$(M, 1-\alpha)$ quantile of order~$\delta$ when testing \eqref{eq:testUConesided}, or above the \linebreak Bin$(M, 1-\alpha)$ quantile of order~$(1-\delta)$ when testing \eqref{eq:testUConesided2}. Since $M$ is large, the same tests are well approximated by rejecting $H_{0i}$ whenever the proportion $n_{1i}^{(\ell)}(\alpha) /M$ of correct coverage is smaller than~$(1-\alpha) - z_{\delta}\sqrt{\alpha(1-\alpha)}$ when testing \eqref{eq:testUConesided}, or larger than $(1-\alpha) + z_{\delta}\sqrt{\alpha(1-\alpha)}$ when testing~\eqref{eq:testUConesided2}, where~$z_{\delta}$ stands for the $(1-\delta)$ standard normal quantile. A two-sided coverage test can also be computed
\begin{equation}\label{eq:2sided}
LR_{\text{cover},i}^{(\ell)}(\alpha):= \big(n_{1i}^{(\ell)}(\alpha) -M(1-\alpha)\big)^2/M\alpha(1-\alpha),
\end{equation}
with asymptotic $\chi^2_{(1)}$ null distribution (as $M\to\infty$).\footnote{It is easily seen that $LR_{\text{cover},i}^{(\ell)}(\alpha)$ is equivalent, up to a constant term, to the so-called ``unconditional coverage'' likelihood ratio test statistic proposed in Section 3.1 of \citet{christoffersen1998} which therefore yields the same results.}
Table \ref{tab:testUC_valid)} reports the empirical rejection frequencies (over $n=90$ series) when testing \eqref{eq:testUConesided} (left-hand panel) and~\eqref{eq:testUConesided2} (right-hand panel), respectively and using the normal approximation of the binomial. The general comments above apply when interpreting those tables: the only valid global conclusions are those resulting from Bonferroni or \v Sid\' ak corrections, which do not lead to any rejections.
\begin{table}[t!]
\centering
\caption{\small Standard \& Poor's 100 Index data ($n=90$ daily returns). Proportion of rejections when testing for valid nominal coverage~\eqref{eq:testUConesided} (left-hand panel) and for sharp nominal coverage \eqref{eq:testUConesided2} (middle panel), and when considering the two-sided test~\eqref{eq:2sided} (right-hand panel)}\label{tab:testUC_valid)}
\vskip .2cm
\footnotesize
\begin{tabular}{l | ccc | ccc | ccc}
\hline
\hline
&\multicolumn{3}{|c}{valid nominal coverage test}&\multicolumn{3}{|c}{sharp nominal coverage test}&\multicolumn{3}{|c}{two-sided coverage test}\\
\hline
$\alpha=0.1$ & $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$& $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$& $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$\\
\hline
$\ell=126$&0.1444 & 0.1222 & 0.0889&0.2444& 0.1444& 0.0333&0.2667 & 0.1889 & 0.0778\\
$\ell=252$&0.1556 & 0.1333 & 0.0778&0.3111& 0.2111& 0.0889&0.3444 & 0.2556 & 0.1333\\
\hline
\hline
$\alpha=0.05$ & $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$& $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$& $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$\\
\hline
$\ell=126$&0.2889& 0.1889& 0.1333&0.0111& 0.0000& 0.0000&0.1889& 0.1556& 0.1000\\
$\ell=252$&0.3000& 0.2000& 0.1444&0.0556& 0.0111& 0.0000&0.2111 & 0.1556 & 0.1333\\
\hline
\hline
\end{tabular}
\end{table}
\subsubsection{Testing against serial dependence}\label{522}
If the alternative of interest is serial dependence among coverage indicators, we propose considering, for each individual stock $i$, alternatives of binary first-order Markov dependence. More precisely, defining the transition probabilities
$$p_{hk,i}(\alpha)=\mathrm{P}\Big(\widehat{\mathcal H}^{(\ell)}_{i,\tau+1|\tau}(\alpha)=k\,\Big\vert\, \widehat{\mathcal H}^{(\ell)}_{i,\tau|\tau-1}(\alpha)=h\Big),\quad h,k=1,0,$$
we consider the testing problem (with unspecified unconditional probability $p_i(\alpha)$ of correct coverage)
\begin{equation}\label{eq:testID}
H_{0i}:p_{01,i}(\alpha)=p_{11,i}(\alpha)=:p_i(\alpha) \;\quad\text{ versus }\quad H_{1i}:p_{01,i}(\alpha)\neq p_{11,i}(\alpha);
\end{equation}
note that $p_{01,i}(\alpha)=p_{11,i}(\alpha)$ automatically
implies $p_{00,i}(\alpha)=p_{10,i}(\alpha)$.
Defining
\begin{align}
n_{11i}^{(\ell)}(\alpha) := \sum_{\tau=T-M+1}^{T-1}\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha)\widehat{\mathcal H}_{i,\tau|\tau-1}^{(\ell)}(\alpha),
\qquad & n_{10i}^{(\ell)}(\alpha) := n_{1i}^{(\ell)}(\alpha)-n_{11i}^{(\ell)}(\alpha),\nonumber\\
n_{01i}^{(\ell)}(\alpha) := \sum_{\tau=T-M+1}^{T-1}\widehat{\mathcal H}_{i,\tau+1|\tau}^{(\ell)}(\alpha)\big(1-\widehat{\mathcal H}_{i,\tau|\tau-1}^{(\ell)}(\alpha)\big),\qquad & n_{00i}^{(\ell)}(\alpha) := n_{0i}^{(\ell)}(\alpha)-n_{01i}^{(\ell)}(\alpha),\nonumber
\end{align}
the statistics
\[
\pi_i^{(\ell)}(\alpha):=\big( {n_{01i}^{(\ell)}(\alpha)+n_{11i}^{(\ell)}(\alpha)}\big)/
M
\]
are estimators of the $p_i(\alpha)$'s under the null, while
\begin{align}
{\pi}_{11i}^{(\ell)}(\alpha) &:= {n_{11i}^{(\ell)}(\alpha)}/{n_{1i}^{(\ell)}(\alpha)}, \qquad\qquad \qquad{\pi}_{10i}^{(\ell)}(\alpha):=1-{\pi}_{11i}^{(\ell)}(\alpha),\nonumber\\
{\pi}_{01i}^{(\ell)}(\alpha) &:=\ {n_{01i}^{(\ell)}(\alpha)}/\big( M-{n_{1i}^{(\ell)}(\alpha)}\big), \ \ \text{and}\quad {\pi}_{00i}^{(\ell)}(\alpha):=1-{\pi}_{01i}^{(\ell)}(\alpha)\nonumber
\end{align}
are estimating the transition probabilities $p_{hk,i}(\alpha)$ under the alternative. Log-likelihoods under the null and the alternative are
\begin{align}
L_{0i}(\alpha)=&\, (n_{00i}^{(\ell)}(\alpha)+n_{10i}^{(\ell)}(\alpha))\log[1-\pi_i^{(\ell)}(\alpha)]+n_{01i}^{(\ell)}(\alpha)+n_{11i}^{(\ell)}(\alpha)\log [\pi_i^{(\ell)}(\alpha)], \nonumber
\end{align}
and
\begin{align}
L_{1i}^{(\ell)}(\alpha)=& \,n_{00i}^{(\ell)}(\alpha)\log [1-{\pi}_{01i}^{(\ell)}(\alpha)]+n_{01i}^{(\ell)}(\alpha)\log[{\pi}_{01i}^{(\ell)}(\alpha)] \nonumber\\
&\hspace{33mm}+n_{10i}^{(\ell)}(\alpha)\log [1-{\pi}_{11i}^{(\ell)}(\alpha)]+n_{11i}^{(\ell)}(\alpha)\log [{\pi}_{11i}^{(\ell)}(\alpha)],\nonumber
\end{align}
respectively. For any given $i$, $\alpha$ and $\ell$, thus, we can construct a likelihood-ratio test for \eqref{eq:testID}, based on the asymptotically $ \chi^2_{(1)}$ null distribution (as $M\to\infty$) of $LR_{\text{{ind}},i}^{(\ell)}(\alpha) := 2\big[L_{1i}^{(\ell)}(\alpha)-L_{0i}(\alpha)\big]$ (see also Section 3.2 in \citealp{christoffersen1998}). More general alternatives, involving higher-order serial dependencies, could be considered as well, based on the tests proposed by \cite{DHM98}.
In Table \ref{tab:testID} (left-hand panel), we report the proportions of rejections (over the~$n$ series) when testing~\eqref{eq:testID}. The same remarks apply as in the interpretation of Table \ref{tab:testUC_valid)}.
\subsubsection{Combined test}\label{523}
Combining the above tests, a likelihood ratio test (given $i$, $\alpha$, and $\ell$) for
\begin{equation}\label{eq:testCC}
H_{0i}:p_{01,i}(\alpha)=p_{11,i}(\alpha)=(1-\alpha) \ \text{ versus }\ H_{1i}:p_{01,i}(\alpha)\neq p_{11,i}(\alpha)\ \text{ or }\ p_{01,i}(\alpha)=p_{11,i}(\alpha)\neq (1-\alpha)
\end{equation}
can be based on the asymptotically $\chi^2_{(2)}$ (as $M\to\infty$) null distribution of
\[
LR_i^{(\ell)}(\alpha) = LR_{\text{cover},i}^{(\ell)}(\alpha)+LR_{\text{{ind}},i}^{(\ell)}(\alpha)
\]
(see also Section 3.3 in \citealp{christoffersen1998}). The fraction of rejections (over $n$ series) when testing~\eqref{eq:testCC} is reported in Table \ref{tab:testID} (right-hand panel). The same remarks as in Table \ref{tab:testUC_valid)} still apply.
\begin{table}[t!]
\centering
\caption{\small Standard \& Poor's 100 Index data ($n=90$ daily returns). Proportion of rejections when testing against serial dependence \eqref{eq:testID} (left-hand panel) and in the combined problem \eqref{eq:testCC} (right-hand panel). }\label{tab:testID}
\vskip .2cm
\footnotesize
\begin{tabular}{l | ccc | ccc}
\hline
\hline
$\alpha=0.1$ & $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$& $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$\\
\hline
$\ell=126$&0.3222& 0.2222& 0.0778&0.3667& 0.2222 & 0.1222\\
$\ell=252$&0.4000& 0.3556 & 0.1889&0.4889& 0.4222 & 0.2556\\
\hline
\hline
$\alpha=0.05$ & $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$& $\delta=0.1$ & $\delta=0.05$ & $\delta=0.01$\\
\hline
$\ell=126$&0.2778& 0.2000& 0.0556&0.2778& 0.2111& 0.1222\\
$\ell=252$&0.3667 & 0.2667 & 0.1667&0.3556& 0.2667& 0.2222\\
\hline
\hline
\end{tabular}
\end{table}
\subsection{Discussion}
In Table \ref{tab:sel_ser}, we report (four panels, according to the values of $\alpha$ and $\ell$) the ten individual series for which the four tests above return the most significant rejections. Rejecting in \eqref{eq:testUConesided} the null hypothesis of a valid coverage (``small'' values of $n_{1i}^{(\ell)}/M$) means that the approximations we are making in the construction of the intervals lead to a loss of prediction accuracy for that specific series: the intervals for that series are not wide enough---equivalently, their actual coverage probability is less than the nominal $(1-\alpha)$ level. The series listed in the first column of each panel thus are ``hardest to predict''. Among them are stocks belonging to the Financial sector, as America International Group (AIG), Bank of America (BAC), and Citigroup (C). These series, in particular, were among those mostly affected by the great financial crisis. Rejecting in \eqref{eq:testUConesided2} the null hypothesis of a sharp coverage (``large'' values of~$n_{1i}^{(\ell)}/M$) also means that the approximations we are making in the construction of the intervals lead to a loss of prediction accuracy for that specific series, now in the sense that we could do better: the intervals for that series are too wide---their actual coverage probability is more than the nominal $(1-\alpha)$ level. The series listed in the second column of each panel thus are ``easiest to predict''. Among them, stocks belonging to the Energy and Consumers sectors, as Exxon Mobil (XOM), Cisco Systems (CSCO), and McDonalds (MCD).
When testing against serial dependence, rejection (``large'' values of $LR_{\text{{ind}},i}^{(\ell)}(\alpha)$) indicates that the predictive information available in past observations has not been fully exploited in the construction of the prediction intervals. This could be the case, for example, if some informative idiosyncratic cross-correlation is available: idiosyncratic cross-correlations indeed are not captured by our univariate autoregressive modelling of idiosyncratic components.
Alternative multivariate models for idiosyncratic components, such as sparse VAR, are likely to improve on this (see e.g. the approach proposed in \citealp{barigozzihallin15c}), and could be incorporated into our two-step GDFM approach. We do not explore this any further in this paper, though. Such dependencies could be related to sectoral co-movements which, being specific to some restricted sector, are not captured by the market-wide factors. This seems to be the case especially for Financial and Energy stocks.
A symptom of that phenomenon is the fact that the explained variance of the common component of the Financial stock returns is about 30\% less than the variance explained by the common component of all other stock returns. The importance of this idiosyncratic variation, which is not accounted for by our approach, may explain why combined tests of correct coverage and independence exhibit, for Financial stock returns, high rejection frequencies.
\begin{sidewaystable}[h!]
\centering
\caption{\small Standard \& Poor's 100 Index data ($n=90$ daily returns). Series tickers for which the null hypotheses considered in Section~5.3 are rejected most significantly. }\label{tab:sel_ser}
\vskip .3cm
\footnotesize
\begin{tabular}{l | llll | l | llll}
\hline
\hline
$\alpha = 0.1$ & smallest & largest & largest & largest &$\alpha = 0.05$ & smallest & largest&largest&largest\\
& $n_{1i}^{(\ell)}(\alpha)/M$ & $n_{1i}^{(\ell)}(\alpha)/M$ & $LR_{\text{{ind}},i}^{(\ell)}(\alpha)$ & $LR_{i}^{(\ell)}(\alpha)$&& $n_{1i}^{(\ell)}(\alpha)/M$ & $n_{1i}^{(\ell)}(\alpha)/M$ & $LR_{\text{{ind}},i}^{(\ell)}(\alpha)$ & $LR_{i}^{(\ell)}(\alpha)$\\
\hline
$\ell=126$ & BAC & MCD & AIG & BAC & $\ell=126$ & SPG & COST & AIG & AIG \\
& SPG & CSCO & BRK.B & SPG & & BAC & MCD & BRK.B & BAC \\
& C & CVX & AMGN & AIG & & C & TGT & AMGN & SPG \\
& AIG & GILD & SPG & C & & AIG & EMC & DVN & C \\
& WFC & MO & BAC & WFC & & WFC & WMT & MRK & WFC \\
& USB & TXN & COP & BRK.B & & BRK.B & CVX & BAC & COP \\
& JPM & WMT & SO & AMGN & & SO & GILD & XOM & BRK.B \\
& COF & XOM & AAPL & USB & & USB & T & CVS & DVN \\
& BRK.B & EMC & APC & COP & & MS & COP & COP & APC \\
& MS & SLB & JNJ & JPM & & COF & CSCO & EXC & USB \\
\hline
\hline
$\alpha = 0.1$ & smallest & largest & largest & largest &$\alpha = 0.05$ & smallest & largest&largest&largest\\
& $n_{1i}^{(\ell)}(\alpha)/M$ & $n_{1i}^{(\ell)}(\alpha)/M$ & $LR_{\text{{ind}},i}^{(\ell)}(\alpha)$ & $LR_{i}^{(\ell)}(\alpha)$&& $n_{1i}^{(\ell)}(\alpha)/M$ & $n_{1i}^{(\ell)}(\alpha)/M$ & $LR_{\text{{ind}},i}^{(\ell)}(\alpha)$ & $LR_{i}^{(\ell)}(\alpha)$\\
\hline
$\ell=252$ & BAC & GILD & AIG & BAC & $\ell=252$ & SPG & MCD & AIG & AIG \\
& SPG & MCD & SPG & SPG & & BAC & GILD & BRK.B & BAC \\
& C & XOM & COP & AIG & & C & ORCL & COP & SPG \\
& AIG & TXN & BRK.B & C & & AIG & QCOM & DVN & C \\
& WFC & CSCO & DVN & WFC & & WFC & CVX & EXC & WFC \\
& USB & CVX & APC & BRK.B & & BRK.B & CSCO & OXY & BRK.B \\
& MS & EBAY & LLY & USB & & JPM & WMT & ALL & COP \\
& BRK.B & EMC & BAC & MS & & USB & COST & BAC & SO \\
& JPM & TGT & UNH & SO & & FCX & EMC & SPG & OXY \\
& COF & MO & C & AMGN & & TWX & TGT & KO & ALL \\
\hline
\hline
\end{tabular}
\end{sidewaystable}
\section{Conclusions}\label{sec:conc}
In this paper, we consider a two-step GDFM approach for jointly modelling stock returns and their volatilities in order to build conditional prediction intervals. A careful study of the consistency properties (as the cross-sectional dimension $n$ and the sample size $T$ both tend to infinity) of the resulting estimators is conducted. Those results are the theoretical foundation of (\citealp{barigozzihallin15a,barigozzihallin15b,barigozzihallin15c}, and \citealp{BHS18}); here, we are using them in the construction of one-step-ahead prediction intervals.
We then apply our methodology to a panel of 90 daily returns of stocks listed in the S\&P100. Through a recursive exercise, we show that we are able to obtain one-step-ahead prediction intervals which are in general more accurate than univariate GARCH methods.
Many extensions of this work are possible, which are left for future research. First, our empirical results indicate that, by exploiting also the cross-sectional lagged dependencies among idiosyncratic components, we could achieve better coverage especially for those series belonging to the Financial sector, which remains strongly interconnected even after controlling for common factors. This could be achieved by computing predictions of idiosyncratic components by fitting multivariate models such as sparse VARs. Second, our methodology immediately allows us to consider bivariate or multivariate prediction intervals. Third, asymmetric prediction intervals can also be considered. In particular, Value-at-Risk indicators are readily computable; moreover, by considering many values of the coverage, we can approximate the whole conditional distribution of returns. Last, another possible application consists in the construction od prediction intervals for macroeconomic variables as GDP or inflation taking into account, in a way similar to \citet{jurado2015}, the uncertainty related to the business cycle.
\small
\bibliography{BH_biblio}
\bibliographystyle{apalike}
\setcounter{section}{0}
\setcounter{subsection}{-1}
\setcounter{equation}{0}
\setcounter{lemma}{0}