Generalized Dynamic Factor Models and Volatilities: Consistency, rates, and prediction intervals

Generalized Dynamic Factor Models and Volatilities: Consistency, Rates, and Prediction Intervals

\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. }

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, tiao1972asymptotic, tiao1978some, tiao1980forecasting, tiao1981modeling, tsay1985use, pena1987identifying, and tiao1989model.

The type of factor model we are considering here is the General or Generalized Dynamic Factor Model (GDFM) introduced by 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 baing02, stockwatson2002, and FLM13. Moreover, as emphasised in fornilippi01 and 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 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., stockwatson2002, baing08JoE, 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. chamberlainrotshild83, CK93, or 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 BLR06 and 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. DN89, ENR92, HRS92, and 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. CKL06 and fan15.} For these reasons, 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 barigozzihallin15b, that two-step GDFM is combined with a GARCH strategy in order to produce point-forecasts for volatilities (see also Trucios19 for a recent example), while barigozzihallin15c and 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 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 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), 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) describes the estimation of the model, and Section (ref) establishes the consistency properties (with rates) of the proposed estimators. In Section (ref), we define the one-step-ahead conditional prediction confidence limits and intervals. In Section (ref), we study the finite-sample properties of our estimators via simulations. Section (ref) 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), we conclude. Proofs are postponed to an Appendix.

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 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$.

A General Dynamic Factor Model for levels and volatilities

We throughout assume that all stochastic variables in this paper belong to the Hilbert space $L_2(\Omega, \mathcal F , \mathrm P)$, where $(\Omega, \mathcal F , \mathrm P)$ is some common probability space. We study double-indexed stochastic processes of the 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(

\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) are mainly taken from FHLZ17, with some modifications, mostly concerning the idiosyncratic components. On the other hand, the assumptions in Section (ref) are new and are related to the log-volatility proxies originally introduced in barigozzihallin15a,barigozzihallin15b.

Model and assumptions for levels

The Generalized Dynamic Factor Model (GDFM) for the levels process $\mbf Y$ is a decomposition of $Y_{it}$ into

with

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}\}$ 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 (ref) in vector notation takes the form

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))$.

More precisely, we assume that (ref)-(ref) hold and satisfy the following assumptions:

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 FLM13 and is empirically verified by boivinng06 and baing08JoE for US macroeconomic data, and by 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) 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. baing02, and 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 hannan1970, for the autocovariances, and the results in Section 6.2 in priestley01, and Theorem 5A in 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 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 fornilippi01 and hallinlippi13, (ref)-(ref) 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:

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) 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 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 (ref) as

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 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.}

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.

As a consequence of Lemma (ref), 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 andersondeistler08 for singular vector processes with a rational spectrum, fornilippi11 and 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

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 (ref), and $\mbf H^\ddag$ an appropriate $(q+1)\times q$ matrix. On that representation, we make the following assumptions.

This assumption allows us to derive an alternative representation of the GDFM (ref) which is particularly useful for estimation and for the construction, in Section (ref) 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 (ref) 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

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 (ref), we have $\left[\mbf A_n(L)\right]^{-1} \mbf H_n= \mbf B_n(L)$ (see Proposition 3 in FHLZ17). Then, the following alternative and equivalent representation of the GDFM holds:

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 (ref).

To conclude with, note that the Yule-Walker equations

characterising the $S$ matrix coefficients of $\mbf A^\ddag(L)$ in (ref) are well defined in view of part {\it (v)} of Assumption (L4); the same conclusion holds, blockwise, for the $n$ -dimensional VAR (ref).

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 (ref) 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 FGLR09 or Assumption 6 in 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$.

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.

Model and assumptions for volatilities

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

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 EM06 and our previous work (barigozzihallin15a,barigozzihallin15b,barigozzihallin15c, and BHS18). In order for such processes to be well defined we make the following 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

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 (ref) in vector notation take the form

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))$.

The following assumptions then are the analogues, for log-volatilities and (ref)-(ref) , of Assumption (L1).

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 UZ11). Pursuing with assumptions, the following one is the log-volatility counterpart of (L2).

Assumptions (V2.iv) implies that we can rewrite (ref) also as

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.

Finally, the analogue (V4) of (L4) again is based on the representation results in FHLZ15: any $(Q+1)$- dimensional subvector $\bm\chi^\ddag_t$ of $\bm\chi_{nt}$ admits an autoregressive representation of the form

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 (ref), and $\mbf R^\ddag$ an appropriate $(Q+1)\times Q$ matrix. On that representation, we make the following assumptions:

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

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) and (ref) for the log-volatility panels.

We then have the following.

\setcounter{equation}{0}

Estimation, consistency, and rates

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”.

Summary of estimation

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 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).

\vskip .2cm

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 (L.iv) and (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 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 FHLZ17 and verified empirically also in 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) 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) on simulated data (see also the results in Appendix (ref)) and in Section (ref) 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 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.

Consistency and rates

Consistency of the estimators of the GDFM model for levels is proved in FHLZ17. Some differences exist, though, between their approach and ours. First, FHLZ17 make slightly weaker assumptions on idiosyncratic serial dependence and, by exploiting results in 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.

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) for more details on the choice of the bandwidths).

Second, the results in 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.

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 vershynin12, and Appendix (ref) 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 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) and (ref) 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 FGLR09 for similar conditions).

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. 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 (L.i)-(L.vi) are as follows.

The proof of parts {\it (a)} and {\it (b)} of Proposition (ref) follows directly from 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) 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 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 (R).

The intuition behind this assumption is as follows. As shown in Appendix (ref), an immediate consequence of Proposition (ref) 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 (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)) 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 (R) and (V.i)-(V.vi) then follows.

This result, which is new, provides the theoretical foundation for the consistency of the estimators used in 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), are proved under Assumption (V2{\it iv}) of a finite-order autoregressive representation for the idiosyncratic components; the same comments as for Proposition (ref) 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.

}

Conditional prediction intervals

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,

where $d_{i0}=1$ because of (ref) and, for the log-volatilities,

where $g_{i0}=1$ because of (ref).

The optimal one-step-ahead linear predictors of level $Y_{it}$ and log-volatility $h_{it}$ are thus

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

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)$ for $\alpha <1-{\rm P}[w_{it}\leq 0]$. Prediction intervals with coverage probability $(1-\alpha )$ can be constructed as

with $\alpha^\pm<1/2$ \special{color cmyk 0 0 0 1.} and $\alpha^- + \alpha^+ =\alpha $, covering $Y_{it|t-1}$ (see (ref)) provided that

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 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 FZ11, Definition 12.1).}

The advantage of quantile-based prediction intervals of the form (ref) 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 (ref). This quantile-based approach moreover allows for unequal tails ($\alpha^- \neq \alpha^+ $ in (ref)---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 \[ \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 \[\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 (ref) and (ref) are

and

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 (ref) to (ref) 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 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)$,

For a given $i$, we say that $ \widehat{\mathcal I}_{i,T+1|T}(\alpha)$ provides the correct coverage if

which is equivalent (see e.g. Lemma 1 in christoffersen1998) to the hypothesis that

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) for details and implementation.

\setcounter{equation}{0}

Simulation study

Setup

To study the performance of our estimator on finite samples, we simulate data ($\mathcal M$ replications) according to the model described in (ref)-(ref).

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$, 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) 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

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

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

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

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 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) 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), 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): 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.

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) 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$.

Results

For each replication, we estimate the model as described in Section (ref). 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) 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

and the maximal errors over all realizations:

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 FHLZ17 and 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). 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) 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), 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.

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). 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,

and

respectively. Results are shown in Table (ref). 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)).

\setcounter{equation}{0}

Interval prediction for S&P100 returns

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 (ref). 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) 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 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

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:

The estimation of the GDFM is based on 10 cross-sectional permutations, as explained at the end of Section (ref). 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 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

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

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) 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)$.

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) 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$.

In Figures (ref) and (ref), 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) 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) 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.

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

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- 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). 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.

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 McN47 test. For given $\alpha$ and $\ell$, consider, for all $i$, the events (discordant GDFM and GARCH coverage results)

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) 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.

Coverage: backtesting

As explained in Section (ref), a formal assessment of the validity of our approach can be based on the backtesting procedure proposed by christoffersen1998. The idea consists in testing the null hypothesis (ref) 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)) in the validity of interval prediction or the sharpness of the nominal coverage level. Else, one may consider (Section (ref)) alternatives of serial dependence. Or, those two issues can be combined (Section (ref)) 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$.

Testing for valid or sharp conditional coverage probabilities

If we are interested in the validity of interval prediction, the relevant testing problems are (one-sided)

If instead we are interested in testing whether $(1-\alpha)$, as a nominal confidence level, is sharp, the testing problems are (still one-sided)

Both testing problems (ref) and (ref), 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 (ref), or above the Bin$(M, 1-\alpha)$ quantile of order $(1-\delta)$ when testing (ref). 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 (ref), or larger than $(1-\alpha) + z_{\delta}\sqrt{\alpha(1-\alpha)}$ when testing (ref), where $z_{\delta}$ stands for the $(1-\delta)$ standard normal quantile. A two-sided coverage test can also be computed

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 christoffersen1998 which therefore yields the same results.}

Table (ref) reports the empirical rejection frequencies (over $n=90$ series) when testing (ref) (left-hand panel) and (ref) (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.

Testing against serial dependence

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)

note that $p_{01,i}(\alpha)=p_{11,i}(\alpha)$ automatically implies $p_{00,i}(\alpha)=p_{10,i}(\alpha)$. Defining

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

are estimating the transition probabilities $p_{hk,i}(\alpha)$ under the alternative. Log-likelihoods under the null and the alternative are

and

respectively. For any given $i$, $\alpha$ and $\ell$, thus, we can construct a likelihood-ratio test for (ref), 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 christoffersen1998). More general alternatives, involving higher-order serial dependencies, could be considered as well, based on the tests proposed by DHM98.

In Table (ref) (left-hand panel), we report the proportions of rejections (over the $n$ series) when testing (ref). The same remarks apply as in the interpretation of Table (ref).

Combined test

Combining the above tests, a likelihood ratio test (given $i$, $\alpha$, and $\ell$) for

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 christoffersen1998). The fraction of rejections (over $n$ series) when testing (ref) is reported in Table (ref) (right-hand panel). The same remarks as in Table (ref) still apply.

Discussion

In Table (ref), 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 (ref) 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 (ref) 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 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.

Conclusions

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 (barigozzihallin15a,barigozzihallin15b,barigozzihallin15c, and 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 jurado2015, the uncertainty related to the business cycle.

\setcounter{section}{0} \setcounter{subsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0}