Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
65,215 characters · 18 sections · 59 citation commands
Sparse Bayesian Time-Varying Covariance Estimation in Many Dimensions
JEL classification: C32; C51; C58
Keywords: dynamic correlation, factor stochastic volatility, curse of dimensionality, shrinkage, minimum variance portfolio
The joint analysis of hundreds or even thousands of time series exhibiting a potentially time-varying variance-covariance structure has been on numerous research agendas for well over a decade. In the present paper we aim to strike the indispensable balance between the necessary flexibility and parameter parsimony by using a factor stochastic volatility (SV) model in combination with a global-local shrinkage prior. Our contribution is threefold. First, the proposed approach offers a hybrid cure to the curse of dimensionality by combining parsimony (through imposing a factor structure) with sparsity (through employing computationally efficient absolutely continuous shrinkage priors on the factor loadings). Second, the efficient construction of posterior simulators allows for conducting Bayesian inference and prediction in very high dimensions via carefully crafted Markov chain Monte Carlo (MCMC) methods made available to end-users through the \proglang{R} r:r package {\normalfont\fontseries{b}\selectfont factorstochvol} r:fac. Third, we show that the proposed method is capable of accurately predicting covariance and precision matrices which we asses via statistical and economic forecast evaluation in several simulation studies and an extensive real-world example.
Concerning factor SV modeling, early key references include har-etal:mul, pit-she:tim, agu-wes:bay which were later picked up and extended by e.g. phi-gli:fac, chi-etal:ana, han:ass, lop-car:fac, nak-wes:dynJFE, zho-etal:bay, ish-omo:por. While reducing the dimensionality of the problem at hand, models with many factors are still rather rich in parameters. Thus, we further shrink unimportant elements of the factor loadings matrix to zero in an automatic way within a Bayesian framework. This approach is inspired by high-dimensional regression problems where the number of parameters frequently exceeds the size of the data. In particular, we adopt the approach brought forward by car-dou:spa, gri-bro:inf who suggest to use a special continuous prior structure -- the Normal-Gamma prior -- on the regression parameters (in our case the factor loadings matrix). This shrinkage prior is a generalization of the Bayesian Lasso par-cas:bay and has recently received attention in the econometrics literature bit-fru:ach, hub-fel:ada.
Another major issue for such high-dimensional problems is the computational burden that goes along with statistical inference, in particular when joint modeling is attempted instead of multi-step approaches or rolling-window-like estimates. Suggested solutions include eng-kel:dyn who propose an estimator assuming that pairwise correlations are equal at every point in time, pal-etal:fit who consider composite likelihood estimation, gru-wes:gpu who use a decoupling-recoupling strategy to parallelize estimation (executed on graphical processors), lop-etal:par who treat the Cholesky-decomposed covariance matrix within the framework of Bayesian time-varying parameter models, and oh-pat:mod who choose a copula-based approach to link separately estimated univariate models. We propose to use a Gibbs-type sampler which allows to jointly take into account both parameter as well as sampling uncertainty in a finite-sample setup through fully Bayesian inference, thereby enabling inherent uncertainty quantification. Additionally, this approach allows for fully probabilistic in- and out-of-sample density predictions.
For related work on sparse Bayesian prior distributions in high dimensions, see e.g. kau-sch:bay who use a point mass prior specification for factor loadings in dynamic factor models or ahe-etal:bay who use a graphical representation of vector autoregressive models to select sparse graphs. From a mathematical point of view, pat-etal:pos investigate posterior contraction rates for a related class of continuous shrinkage priors for static factor models and show excellent performance in terms of posterior rates of convergence with respect to the minimax rate. All of these works, however, assume homoskedasticity and are thus potentially misspecified when applied to financial or economic data. For related methods that take into account heteroskedasticity, see e.g. nak-wes:dynJFE, nak-wes:dynBJPS who employ a latent thresholding process to enforce time-varying sparsity. Moreover, zha-etal:dyn approach this issue via dependence networks, lod-etal:sel use stochastic search for model selection, and bas-etal:tim use time-varying combinations of dynamic models and equity momentum strategies. These methods are typically very flexible in terms of the dynamics they can capture but are applied to moderate dimensional data only.
We illustrate the merits of our approach through extensive simulation studies and an in-depth financial application using 300 S&P 500 members. In simulations, we find considerable evidence that the Normal-Gamma shrinkage prior leads to substantially sparser factor loadings matrices which in turn translate into more precise correlation estimates when compared to the usual Gaussian prior on the loadings.\footnote{Note that in contrast to e.g.\ fru-tue:bay, we do not attempt to identify exact zeros in a covariance matrix. Rather, we aim to find a parsimonious factor representation of the underlying heteroskedastic data which may (or may not) imply covariances that are close to zero at times.} In the real-world application, we evaluate our model against a wide range of alternative specifications via log predictive scores and minimum variance portfolio returns. Factor SV models with sufficiently many factors turn out to imply extremely competitive portfolios in relation to well-established methods which typically have been specifically tailored for such applications. Concerning density forecasts, we find that our approach outperforms all included competitors by a large margin.
The remainder of this paper is structured as follows. In Section (ref), the factor SV model is specified and the choice of prior distributions is discussed. Section (ref) treats statistical inference via MCMC methods and sheds light on computational aspects concerning out-of-sample density predictions for this model class. Extensive simulation studies are presented in Section (ref), where the effect of the Normal-Gamma prior on correlation estimates is investigated in detail. In Section (ref), the model is applied to 300 S&P 500 members. Section (ref) wraps up and points out possible directions for future research.
Consider an $m$-variate zero-mean return vector ${\bm{ y}}_t = (y_{1t}, \ldots, y_{m t}) ^{'}$ for time $t = 1, \dots, T$ whose conditional distribution is Gaussian, i.e. \[ {\bm{ y}}_t|{\bm{\Sigma}}_t \sim \mathcal{N}_{m}\!\left({\bm{{0}}}, {\bm{\Sigma}}_t\right). \]
To reduce dimensionality, factor SV models utilize a decomposition of the $m \times m$ covariance matrix ${\bm{\Sigma}}_t$ with $m(m+1)/2$ free elements into a factor loadings matrix $\bm{\Lambda}$ of size $m \times r$, an $r$-dimensional diagonal matrix ${\bm{ V}}_t$ and an $m$-dimensional diagonal matrix ${\bm{{{U}}}}_t$ in the following fashion:
This reduces the number of free elements to $m r + m + r$. Because $r$ is typically chosen to be much smaller than $m$, this specification constrains the parameter space substantially, thereby inducing parameter parsimony. For the paper at hand, $\bm{\Lambda}$ is considered to be time invariant whereas the elements of both ${\bm{ V}}_t$ and ${\bm{{{U}}}}_t$ are allowed to evolve over time through parametric stochastic volatility models, i.e.\ ${\bm{{{U}}}}_t = \mbox{\text diag}\!\left(\exp(h_{1t}),\ldots,\exp(h_{m t})\right)$ and ${\bm{ V}}_t = \mbox{\text diag}\!\left(\exp(h_{m+1,t}),\ldots,\exp(h_{m+r,t})\right)$ with
More specifically, ${\bm{{{U}}}}_t$ describes the idio\-syncratic (series-specific) variances while ${\bm{ V}}_t$ contains the variances of underlying orthogonal factors ${\bm{ f}}_t \sim \mathcal{N}_{r}\!\left({\bm{{0}}}, {\bm{ V}}_t\right)$ that govern the contemporaneous dependence. The autoregressive process in ((ref)) is assumed to have mean zero to identify the unconditional scaling of the factors.
This setup is commonly written in the following hierarchical form chi-etal:ana:
where the distributions are assumed to be conditionally independent for all points in time. To make further exposition clearer, let ${\bm{ y}} = ({\bm{ y}}_1 \cdots {\bm{ y}}_T)$ denote the $m \times T$ matrix of all observations, ${\bm{ f}}=({\bm{ f}}_1 \cdots {\bm{ f}}_T)$ the $r \times T$ matrix of all latent factors and ${\bm{h}} = ({\bm{h}}_1'\cdots{\bm{h}}_{m+r}')'$ the $(T+1)\times(m+r)$ matrix of all $m+r$ log-variance processes ${\bm{h}}_i=(h_{i0}, h_{i1}, \ldots,h_{iT})$, $i = 1,\dots,m+r$. The vector $\bm{\theta}_i = (\mu_i, \phi_i, \sigma_i)'$ is referred to as the vector of parameters where $\mu_i$ is the level, $\phi_i$ the persistence, and $\sigma_i^2$ the innovation variance of ${\bm{h}}_i$. To denote specific rows and columns of matrices, we use the “dot” notation, i.e. $\bm{X}_{i\cdot}$ refers to the $i$th row and $\bm{X}_{\cdot j}$ to the $j$th column of $\bm{X}$. The proportions of variances explained through the common factors for each component series, $C_{it} = 1-{{{U}}_{ii,t}}/{\Sigma_{ii,t}}$ for $i = 1,\dots,m,$ are referred to as the communalities. Here, ${{U}}_{ii,t}$ and $\Sigma_{ii,t}$ denote the $i$th diagonal element of ${\bm{{{U}}}}_t$ and ${\bm{\Sigma}}_t$, respectively. As by construction $0 \leq {{U}}_{ii,t} \leq \Sigma_{ii,t}$, the communality for each component series and for all points in time lies between zero and one. The joint (overall) communality $C_t = m^{-1}\sum_{i=1}^m C_{it}$ is simply defined as the arithmetic mean over all series.
Three comments are in order. First, the variance-covariance decomposition in ((ref)) can be rewritten as ${\bm{\Sigma}}_t = \bm{\Lambda}_t \bm{\Lambda}'_t + {\bm{{{U}}}}_t$ with $\bm{\Lambda}_t := \bm{\Lambda} {\bm{ V}}_t^{1/2}$. An essential assumption within the factor framework is that both ${\bm{ V}}_t$ as well as ${\bm{{{U}}}}_t$ are diagonal matrices. This implies that the factor loadings $\bm{\Lambda}_t$ are dynamic but can only vary column-wise over time. Consequently, the time-variability of ${\bm{\Sigma}}_t$'s off-diagonal elements are cross-sectionally restricted while its diagonal elements are allowed to move independently across series. Hence, the “strength” of a factor, i.e.\ its cross-sectional explanatory power, varies jointly for all series loading on it. Consequently, it is likely that more factors are needed to properly explain the co-volatility dynamics of a multivariate time series than in models which allow for completely unrestricted time-varying factor loadings lop-car:fac, correlated factors zho-etal:bay, or approximate factor models bai-ng:det. Our specification, however, is less prone to overfitting and has the significant advantage of vastly simplified computations.
Second, identifying loadings for latent factor models is a long-standing issue that goes back to at least and-rub:sta who discuss identification of factor loadings. Even though this problem is alleviated somewhat when factors are allowed to exhibit conditional heteroskedasticity sen-fio:ide, rig:ide, most authors have chosen an upper triangular constraint of the loadings matrix with unit diagonal elements, thereby introducing dependence on the ordering of the data fru-lop:par. However, when estimation of the actual factor loadings is not the primary concern (but rather a means to estimate and predict the covariance structure), this issue is less striking because a unique identification of the loadings matrix is not necessary.\footnote{The conditional covariance matrix ${\bm{\Sigma}}_t = \bm{\Lambda} {\bm{ V}}_t \bm{\Lambda}' + {\bm{{{U}}}}_t$ involves a rotation-invariant transformation of $\bm{\Lambda}$.} This allows leaving the factor loadings matrix completely unrestricted, thus rendering the method invariant with respect to the ordering of the series.
Third, note that even though the joint distribution of the data is conditionally Gaussian, its stationary distribution has thicker tails. Nevertheless, generalizations of the univariate SV model to cater for even more leptokurtic distributions lie-jun:sto or asymmetry yu:on can straightforwardly be incorporated in the current framework. All of these extensions, however, tend to increase both sampling inefficiency as well as running time considerably and could thus preclude inference in very high dimensions.
The usual prior for each (unrestricted) element of the factor loadings matrix is a zero-mean Gaussian distribution, i.e.\ $\Lambda_{ij} \sim \mathcal{N}\!\left(0, \tau^2_{ij}\right)$ independently for each $i$ and $j$, where $\tau^2_{ij} \equiv \tau^2$ is a constant specified a priori pit-she:tim,agu-wes:bay,chi-etal:ana,ish-omo:por,kas-etal:eff. To achieve more shrinkage, we model this variance hierarchically by placing a hyperprior on $\tau^2_{ij}$. This approach is related to bha-dun:spa,pat-etal:pos who investigate a similar class of priors for homoskedastic factor models. More specifically, let
Intuitively, each prior variance $\tau_{ij}^2$ provides element-wise shrinkage governed independently for each row by $\lambda_i^2$. Integrating out $\tau^2_{ij}$ yields a density for $\Lambda_{ij}|\lambda_i^2$ of the form $ p(\Lambda_{ij}|\lambda_i^2) \propto |\Lambda_{ij}|^{a_i-1/2}K_{a_i-1/2}(\sqrt a_i \lambda_i|\Lambda_{ij}|), $ where $K$ is the modified Bessel function of the second kind. This implies that the conditional variance of $\Lambda_{ij}|\lambda_i^2$ is $2/\lambda_i^2$ and the excess kurtosis of $\Lambda_{ij}$ is $3/a_i$. The hyperparameters $a_i$, $c_i$, and $d_i$ are fixed a priori, whereas $a_i$ in particular plays a crucial role for the amount of shrinkage this prior implies. Choosing $a_i$ small enforces strong shrinkage towards zero, while choosing $a_i$ large imposes little shrinkage. For more elaborate discussions on Bayesian shrinkage in general and the effect of $a_i$ specifically, see gri-bro:inf and pol-sco:shr. Note that the Bayesian Lasso prior par-cas:bay arises as a special case when $a_i = 1$.
One can see prior ((ref)) as row-wise shrinkage with element-wise adaption in the sense that all variances in row $i$ can be thought of as “random effects” from the same underlying distribution. In other words, each series has high and a priori independent mass not to load on any factors and thus can be thought of as series-specific shrinkage. For further aspects on introducing hierarchical prior structure via the Normal-Gamma distribution, see gri-bro:hie, hub-fel:ada. Analogously, it turns out to be fruitful to also consider column-wise shrinkage with element-wise adaption, i.e.
This means that each factor has high and a priori independent mass not to be loaded on by any series and thus can be thought of as factor-specific shrinkage.
Concerning the univariate SV priors, we follow kas-fru:anc. For the $m$ idiosyncratic and $r$ factor volatilities, the initial states $h_{i0}$ are distributed according to the stationary distributions of the AR($1$) processes ((ref)) and ((ref)), respectively. Furthermore, $p(\mu_i,\phi_i,\sigma_i)$ = $p(\mu_i)p(\phi_i)p(\sigma_i)$, where the level $\mu_i \in \mathbb{R}$ is equipped with the usual Gaussian prior $\mu_i \sim \mathcal{N}\!\left(b_{\mu}, B_{\mu}\right)$, the persistence parameter $\phi_i \in (-1,1)$ is implied by $(\phi_i+1)/2 \sim \mathcal{B}\!\left(a_0, b_0\right)$ and the volatility of volatility parameter $\sigma_i \in \mathbb{R}^+$ is chosen according to $\sigma_i^2 \sim B_\sigma \chi^2_1= \mathcal{G}\!\left(1/2,1/(2B_\sigma)\right)$.
There are a number of methods to estimate factor SV models such as quasi-maximum likelihood har-etal:mul, simulated maximum likelihood lie-ric:cla,jun-koo:msv, and Bayesian MCMC simulation pit-she:tim, agu-wes:bay, chi-etal:ana, han:ass. For high dimensional problems of this kind, Bayesian MCMC estimation proves to be a very efficient estimation method because it allows simulating from the high dimensional joint posterior by drawing from lower dimensional conditional posteriors.
One substantial advantage of MCMC methods over other ways of learning about the posterior distribution is that it constitutes a modular approach due to the conditional nature of the sampling steps. Consequently, conditionally on the matrix of variances $\bm{\tau} = (\tau_{ij})_{1\leq i\leq m;\, 1\leq j\leq r}$, we can adapt the sampling steps of kas-etal:eff. For obtaining draws for $\bm{\tau}$, we follow gri-bro:inf. The MCMC sampling steps for the factor SV model are:
The above sampling steps are implemented in an efficient way within the \proglang{R} package {\normalfont\fontseries{b}\selectfont factorstochvol} r:fac. Table (ref) displays the empirical run time in milliseconds per MCMC iteration. Note that using more efficient linear algebra routines such as Intel MKL leads to substantial speed gains only for models with many factors. To a certain extent, computation can further be sped up by computing the individual steps of the posterior sampler in parallel. In practice, however, doing so is only useful in shared memory environments (e.g.\ through multithreading/multiprocessing) as the increased communication overhead in distributed memory environments easily outweighs the speed gains.
Given draws of the joint posterior distribution of parameters and latent variables, it is in principle straightforward to predict future covariances and consequently also future observations. This gives rise to the predictive density gew-ami:com, defined as
where $\bm{\kappa}$ denotes the vector of all unobservables, i.e.\ parameters and latent variables. The superscript $o$ in $\bm{y}^o_{[1:t]}$ denotes ex post realizations (observations) for the set of points in time $\{1,\dots ,t\}$ of the ex ante random values $\bm{y}_{[1:t]} = (\bm{y}_1 \cdots \bm{y}_t)$. The integration space ${\bm K}$ simply stands for the space of the possible values for $\bm{\kappa}$. Because ((ref)) is the integral of the likelihood function where the values of $\bm{\kappa}$ are weighted according to their posterior distribution, it can be seen as the forecast density for an unknown value $\bm{y}_{t+1}$ after accounting for the uncertainty about $\bm{\kappa}$, given the history $\bm{y}^o_{[1:t]}$.
As with most quantities of interest in Bayesian analysis, computing the predictive density can be challenging because it constitutes an extremely high-dimensional integral which cannot be solved analytically. However, it may be approximated at a given “future” point $\bm{y}^f$ through Monte Carlo integration,
where $\bm{\kappa}^{(k)}_{[1:t]}$ denotes the $k${th} draw from the posterior distribution up to time $t$. If ((ref)) is evaluated at $\bm{y}^f=\bm{y}_{t+1}^o$, it is commonly referred to as the (one-step-ahead) predictive likelihood at time $t+1$, denoted $P\!L_{t+1}$. Also, draws from ((ref)) can straightforwardly be obtained by generating values $\bm{y}_{t+1}^{(k)}$ from the distribution given through the (in our case multivariate Gaussian) density $p(\bm{y}_{t+1}|\bm{y}^o_{[1:t]}, \bm{\kappa}^{(k)}_{[1:t]})$.
For the model at hand, two ways of evaluating the predictive likelihood particularly stand out. First, one could average over $k=1,\dots, K$ densities of $$\mathcal{N}_{m}\!\left(\bm{\Lambda}^{(k)}_{[1:t]}{\bm{ f}}^{(k)}_{t+1,[1:t]}, {\bm{{{U}}}}^{(k)}_{t+1,[1:t]}\right),$$ evaluated at $\bm{y}^o_{t+1}$, where the subscript $t+1$ denotes the corresponding one-step ahead predictive draws and ${\bm{{{U}}}}^{(k)}_{t+1,[1:t]} = \mbox{\text diag}\!\left(\exp{h_{1,t+1,[1:t]}^{(k)}}, \dots, \exp{h_{m,t+1,[1:t]}^{(k)}}\right)$. Note that because ${\bm{{{U}}}}^{(k)}_{t+1,[1:t]}$ is by construction diagonal, this method only requires univariate Gaussian evaluations and is thus computationally efficient. Nevertheless, because evaluation is done conditionally on realized values of ${\bm{ f}}_{t+1,[1:t]}$, it is extremely unstable in many dimensions. Moreover, since the numerical inaccuracy increases with an increasing number of factors $r$, this approach can lead to systematic undervaluation of $P\!L_{t+1}$ for larger $r$. Thus, in what follows, we recommended an alternative approach.
To obtain $P\!L_{t+1}$, we suggest to average over $k=1,\dots, K$ densities of $$\mathcal{N}_{m}\!\left({\bm{{0}}}, \bm{\Lambda}^{(k)}_{[1:t]} {\bm{ V}}^{(k)}_{t+1,[1:t]} (\bm{\Lambda}^{(k)}_{[1:t]})' + {\bm{{{U}}}}^{(k)}_{t+1,[1:t]}\right),$$ evaluated at $\bm{y}^o_{t+1}$, where ${\bm{ V}}^{(k)}_{t+1,[1:t]} = \mbox{\text diag}\!\left(\exp{h_{m+1,t+1,[1:t]}^{(k)}}, \dots, \exp{h_{m+r,t+1,[1:t]}^{(k)}}\right)$. This form of the predictive likelihood is obtained by analytically performing integration in ((ref)) with respect to ${\bm{ f}}_{t+1,[1:t]}$. Consequently, it is numerically more stable, irrespectively of the number of factors $r$. However, it requires a full $m$-variate Gaussian density evaluation for each $k$ and is thus computationally much more expensive. To a certain extent, the computational burden can be mitigated by using the Woodbury matrix identity, ${\bm{\Sigma}}_t^{-1} = {\bm{{{U}}}}_t^{-1} - {\bm{{{U}}}}_t^{-1}\bm{\Lambda} \left({\bm{ V}}_t^{-1}+\bm{\Lambda}'{\bm{{{U}}}}_t^{-1}\bm{\Lambda} \right)^{-1} \bm{\Lambda}'{\bm{{{U}}}}_t^{-1}$, along with the matrix determinant lemma, $\det({\bm{\Sigma}}_t) = \det({\bm{ V}}_t^{-1} + \bm{\Lambda}'{\bm{{{U}}}}_t^{-1}\bm{\Lambda}) \det({\bm{ V}}_t)\det({\bm{{{U}}}}_t)$. This substantially speeds up the repetitive evaluation of the multivariate Gaussian distribution if $r \ll m$.
We apply these results for comparing competing models $A$ and $B$ between time points $t_1$ and $t_2$ and consider cumulative log predictive Bayes factors defined through $\log B\!F_{t_1,t_2}(A, B) = \sum_{t=t_1+1}^{t_2} \log P\!L_t(A) - \log P\!L_t(B)$, where $P\!L_t(A)$ and $P\!L_t(B)$ denote the predictive likelihood of model $A$ and $B$ at time $t$, respectively. When the cumulative log predictive Bayes factor is greater than 0 at a given point in time, there is evidence in favor of model $A$, and vice versa. Thereby, data up to time $t_1$ is regarded as prior information and out-of-sample evaluation starts at time $t_1+1$.
The aim of this section is to apply the model to a simulated data set in order to illustrate the shrinkage properties of the Normal-Gamma prior for the factor loadings matrix elements. For this purpose, we first illustrate several scenarios on a single ten dimensional data set. Second, we investigate the performance of our model in a full Monte Carlo simulation based on $100$ simulated data sets. Third, and finally, we investigate to what extend these results carry over to higher dimensions.
In what follows, we compare five specific prior settings. Setting 1 refers to the usual standard Gaussian prior with variance $\tau_{ij}^2 \equiv \tau^2 = 1$ and constitutes the benchmark. Setting 2 is the row-wise Bayesian Lasso where $a_i = 1$ for all $i$. Setting 3 is the column-wise Bayesian Lasso where $a_j=1$ for all $j$. Setting 4 is the Normal-Gamma prior with row-wise shrinkage where $a_i = 0.1$ for all $i$. Setting 5 is the Normal-Gamma prior with column-wise shrinkage where $a_j = 0.1$ for all $j$. Throughout this section, prior hyperparameters are chosen as follows: $b_\mu = 0$, $B_\mu = 1000$, $B_\sigma = 1$. The prior hyperparameters for the persistence of the latent log variances are fixed at $a_0 = 10$, $b_0 = 2.5$ for the idiosyncratic volatilities and $a_0 = 2.5$, $b_0 = 2.5$ for the factor volatilities; note that the parameters of the superfluous factor are only identified through the prior. The shrinkage hyperparameters are set as in bel-etal:hie, i.e.\ $c_i=c_j=d_i=d_j=0.001$ for all applicable $i$ and $j$. For each setting, the algorithm is run for $110\,000$ iterations of which the first $10\,000$ draws are discarded as burn-in.
To investigate the effects of different priors on the posteriors of interest, we simulate a single data set from a two factor model for $m=10$ time series of length $T=1000$. For estimation, an overfitting model with three latent factors is employed. The nonzero parameter values used for simulation are picked randomly and are indicated as black circles in Figure (ref); some loadings are set to zero, indicated by black dots. We set $\Lambda_{ij}$ to zero if $j>i$ for simulation and estimation.
Figure (ref) shows smoothed kernel density estimates of posterior loadings under the different prior assumptions. The signs of the loadings have not been identified so that a multimodal posterior distribution hints at a “significant” loading whereas a unimodal posterior hints at a zero loading, see also fru-wag:sto. It stands out that only very little shrinkage is induced by the standard Gaussian prior. The other priors, however, impose considerably tighter posteriors. For the nonzero loadings on factor one, e.g., the row-wise Bayesian Lasso exhibits the strongest degree of shrinkage. Little difference between the various shrinkage priors can be spotted for the nonzero loadings on factor two.
Turning towards the zero loadings, the strongest shrinkage is introduced by both variants of the Normal-Gamma prior, followed by the different variants of the Bayesian Lasso and the standard Gaussian prior. This is particularly striking for the loadings on the superfluous third factor. The difference between row- and column-wise shrinkage for the Lasso variants can most clearly be seen in row 9 and column 3, respectively. The row-wise Lasso captures the “zero-row” 9 better, while the column-wise Lasso captures the “zero-column” 3 better. Because of the increased element-wise shrinkage of the Normal-Gamma prior, the difference between the row-wise and the column-wise variant are minimal.
In the context of covariance modeling, however, factor loadings can be viewed upon as a mere means to parsimony, not the actual quantity of interest. Thus, Figure (ref) displays selected time-varying correlations. The top panel shows a posterior interval estimate (mean plus/minus two standard deviations) for the correlation of series 1 and series 2 (which is nonzero) under all five prior settings; the bottom panel depicts the interval estimate for the correlation of series 9 and 10 (which is zero). While the relative differences between the settings in the nonzero correlation case are relatively small, the zero correlation case is picked up substantially better when shrinkage priors are used. Posterior means are closer to zero and the posterior credible intervals are tighter.
To conclude, we briefly examine predictive performance by investigating cumulative log predictive Bayes factors. Thereby, the first $1000$ points in time are treated as prior information, then 1-day- and 10-days-ahead predictive likelihoods are recursively evaluated until $t=1500$. Table (ref) displays the sum of these values for the respective models in relation to the 2-factor model with the standard Gaussian prior. This way, numbers greater than zero can be interpreted as evidence in favor of the respective model. Not very surprisingly, log Bayes factors are highest for the 2-factor model; within this class, models imposing stronger shrinkage perform slightly better, in particular when considering the longer 10-day horizon. Underfitting models predict very poorly both on the short and the longer run, while overfitting models appear almost en par with the baseline model when shrinkage priors are used. This suggests that shrinkage safeguards against overfitting, at least to a certain extent.
For a more comprehensive understanding of the shrinkage effect, the above study is repeated for $100$ different data sets where all latent variables are generated randomly for each realization. In Table (ref), the medians of the respective relative RMSEs (root mean squared errors, averaged over time) between the true and the estimated pairwise correlations are depicted. The part above the diagonal represents the relative performance of the row-wise Lasso prior (setting 2) with respect to the baseline prior (setting 1), the part below the diagonal represents the relative performance of the row-wise Normal-Gamma prior (setting 4) with respect to the row-wise Lasso prior (setting 2). Clearly, gains are highest for series 9 which is by construction completely uncorrelated to the other series. Additionally, geometric averages of these performance indicators are displayed in the first row (setting 2 vs.\ baseline) and in the last row (setting 4 vs.\ baseline). They can be seen as the average relative performance of one specific series' correlation estimates with all other series.
To illustrate the fact that extreme choices of $c_i$ and $d_i$ are crucial for the shrinkage effect of the Bayesian Lasso, Table (ref) displays relative RMSEs for moderate hyperparameter choices $c_i = d_i = 1$. Note that the performance of the Bayesian Lasso deteriorates substantially while performance of the Normal-Gamma prior is relatively robust with regard to these choices. This indicates that the shrinkage effect of the Bayesian Lasso is strongly dependent on the particular choice of these hyperparameters (governing row-wise shrinkage) while the Normal-Gamma can adapt better through increased element-wise shrinkage.
An overall comparison of the errors under different priors is provided in Table (ref) which lists RMSEs and MAEs for all prior settings, averaged over the non-trivial correlation matrix entries as well as time. Note again that results under the Lasso prior are sensitive to the particular choices of the global shrinkage hyperparameters as well as the choice of row- or column-wise shrinkage, which is hardly the case for the Norma-Gamma prior. Interestingly, the performance gains achieved through shrinkage prior usage are higher when absolute errors are considered. This is coherent with the extremely high kurtosis of Normal-Gamma-type priors which, while placing most mass around zero, allow for large values.
The findings are similar if dimensionality is increased; in analogy to above, we report overall RMSEs and MAEs for $495\,000$ pairwise correlations, resulting from $m=100$ component series at $T=1000$ points in time. The factor loadings for the $r = 10$ factors are again randomly sampled with 43.8% of the loadings being equal to zero, resulting in about 2.6% of the pairwise correlations being zero. Using this setting, $100$ data sets are generated; for each of these, a separate (overfitting) factor SV model using $r = 11$ factors without any prior restrictions on the factor loadings matrix is fit. The error measures are computed and aggregated. Table (ref) reports the medians thereof. In this setting, the shrinkage priors outperform the standard Gaussian prior by a relatively large margin; the effect of the specific choice of the global shrinkage hyperparameters is less pronounced.
In this section we apply the SV factor model to stock prices listed in the Standard & Poor's 500 index. We only consider firms which have been continuously included in the index from November 1994 until December 2013, resulting in $m=300$ stock prices on $5001$ days, ranging from 11/1/1994 to 12/31/2013. The data was obtained from Bloomberg Terminal in January 2014. Instead of considering raw prices we investigate percentage log-returns which we demean a priori.
The presentation consists of two parts. First, we exemplify inference using a multivariate stochastic volatility model and discuss the outcome. Second, we perform out-of-sample predictive evaluation and compare different models. To facilitate interpretation of the results discussed in this section, we consider the GICS\footnote{Global Industry Classification Standard, retrieved from \url{https://en.wikipedia.org/w/index.php?title=List_of_S%26P_500_companies&oldid=589980759} on April 11, 2016.} classification into 10 sectors listed in Table (ref).
To keep graphical representation feasible, we only focus on the latest 2000 returns of our data set, i.e.\ 5/3/2006 to 12/31/2013. This time frame is chosen to include both the 2008 financial crisis as well as the period before and thereafter. Furthermore, we restrict our discussion to a four-factor model. This choice is somewhat arbitrary but allows for a direct comparison to a popular model based on four observed (Fama-French plus Momentum) factors. A comparison of predictive performance for varying number of factors is discussed in Section (ref); the Fama-French plus Momentum model is introduced in Section (ref).
We run our sampler employing the Normal-Gamma prior with row-wise shrinkage for $110\,000$ draws and discard the first $10\,000$ draws as burn-in.\footnote{To keep presentation at a reasonable length and because qualitative as well as quantitative results are very similar, we omit details about the Normal-Gamma prior with column-wise shrinkage.} Of the remaining $100\,000$ draws every $10$th draw is kept, resulting in $10\,000$ draws used for posterior inference. Hyperparameters are set as follows: $a_i \equiv a = 0.1$, $c_i \equiv c=1$, $d_i \equiv d = 1$, $b_\mu=0$, $B_\mu=100$, $a_0=20$, $b_0=1.5$, $B_\sigma=1, B_{m+j} = 1$ for $j = 1,\dots,r$. To prevent factor switching, we set all elements above the diagonal to zero. The leading series are chosen manually after a preliminary unidentified run such that series with high loadings on that particular factor (but low loadings on the other factors) become leaders. Note that this intervention (which introduces an order dependency) is only necessary for interpreting the factor loadings matrix but not for covariance estimation or prediction. Concerning MCMC convergence, we observe excellent mixing for both the covariance as well as the correlation matrix draws. To exemplify, trace plots of the first $1000$ draws after burn-in and thinning for posterior draws of the log determinant distribution of the covariance and correlation matrices at $t = T$ are displayed in Figure (ref).
To illustrate the substantial degree of volatility co-movement, mean posterior variances are displayed in the top panel of Figure (ref). This depiction resembles one where all series are modeled with independent univariate stochastic volatility models. Clear spikes can be spotted during the financial crisis in late 2008 but also in early 2010 and late 2011. This picture is mirrored (to a certain extent) in the bottom panel which displays the posterior distribution of the joint communality $C_{t}$. In particular during the financial crisis, the first half of 2010 and late 2011 the joint communality reaches high values of 0.7 and more.
Median posterior factor loadings are visualized in Figure (ref). In the top panel it can be seen that all series significantly load on the first factor which consequently could be interpreted to represent the joint dynamics of the wider US equity market. Highly loading elements include United States Steel Corp.\ (X) and Cliffs Natural Resources Inc.\ (CLF), both of which belong to the sector Materials and both of which have been dropped from the S&P 500 index in 2014 due to market capitalization changes. Cummins Inc.\ (CMI, Industrials) and PulteGroup, Inc.\ (PHM, Consumer Discretionary) rank third and fourth. Companies in sectors Consumer Staples, Utilities and Health Care tend to load comparably low on this factor.
Investigating the second factor, it stands out that due to the use of the Normal-Gamma prior a considerable amount of loadings are shrunk towards zero. Main drivers are all in the sector Utilities. Also, companies in sectors Consumer Staples, Health Care and (to a certain extent) Financials load positively here. Both the loadings on factor 3 as well as the loadings on factor 4 are substantially shrunk towards zero. Notable exceptions are Energy and Materials companies for factor 3 and \emph{Financials} for factor 4.
The corresponding factor log variances are displayed in Figure (ref). Apart from featuring similar low- to medium frequency properties, each process exhibits specific characteristics. First, notice the sharp increase of volatility in early 2010 which is mainly visible for the “overall” factor 1. The second factor (Utilities) displays a pre-crisis volatility peak during early 2008. The third factor, driven by Energy and Materials, shows relatively smooth volatility behavior while the fourth factor, governed by the Financial, exhibits a comparably “nervous” volatility evolution.
Finally, we show three examples of the posterior mean of the correlation matrix ${\bm{\Sigma}}_t$ in Figure (ref). The series are grouped according to the alphabetically ordered industry sectors (and simply sorted according to their ticker symbol therein). An animation displaying the mean correlation matrix for all points in time is available at \url{https://vimeo.com/217021226}.
Considering the last trading day in 2006, highly correlated clusters appear within Energy and Utilities, to a certain extent also within Financials, Industrials and Materials. Not very surprisingly, there exists only low correlation between companies in the sectors Consumer Discretionary/Staples and \emph{Energy} but higher correlation between \emph{Energy}, \emph{Industrials} and \emph{Materials}. Looking at the last trading day of 2008, the overall picture changes radically. Higher correlation can be spotted throughout, both within sectors but also between sectors. There are only few companies that show little and virtually no companies that show no correlation with others. Another two years later, we again see a different overall picture. Lower correlations throughout become apparent with moderate correlations remaining within the sectors \emph{Energy}, \emph{Utilities}, and in particular \emph{Financials}.
Even for univariate volatility models evaluating in- or out-of-sample fit is not straightforward because the quantity of interest (the conditional standard deviation) is not directly observable. While in lower dimensions this issue can be circumvented to a certain extent using intraday data and computing realized measures of volatility, the difficulty becomes more striking when the dimension increases. Thus, we focus on iteratively predicting the observation density out-of-sample which is then evaluated at the actually observed values. Because this approach involves re-estimating the model for each point in time, it is computationally costly but can be parallelized in a trivial fashion on multi-core computers.
For the S&P 500 data set, we begin by using the first $3000$ data points (until 5/2/2006) to estimate the one-day-ahead predictive likelihood for day $3001$ as well as the ten-day-ahead predictive likelihood for day $3010$. In a separate estimation procedure, the first $1001$ data points (until 5/3/2006) are used to estimate the one-day-ahead predictive likelihood for day $3002$ and the corresponding ten-day-ahead predictive likelihood for day $3011$, etc. This procedure is repeated for $1990$ days until the end of the sample is reached.
We use a no-factor model as the baseline which corresponds to $300$ individual stochastic volatility models fitted to each component series separately. For each date, values greater than zero mean that the model outperforms the baseline model up to that point in time. Competitors of the no-factor SV model are $r$-factor SV models with $r = 1,\dots,20$ under the usual standard Gaussian prior and under the Normal-Gamma prior with $a_i \equiv 0.1$ and $c_i \equiv d_i \equiv 1$ employing row-wise-shrinkage. All other parameters are kept identical, i.e. $b_\mu=0$, $B_\mu=100$, $a_0=20$, $b_0=2.5$ (idiosyncratic persistences), $a_0=2.5$, $b_0=2.5$ (factor persistences), $B_\sigma=0.1, B_{m+j} = 0.1$ for $j = 1,\dots,r$. Note that because the object of interest in this exercise does not require the factor loadings matrix to be identified, no a priori restrictions are placed on $\bm{\Lambda}$. This alleviates the problem of arranging the data in any particular order before running the sampler. Other competing models are discussed in the following sections.
Accumulated log predictive likelihoods for the entire period are displayed in Figure (ref). Gains in predictive power are substantial up to around 8 factors with little difference for the two priors. After this point, the benefit of adding even more factors turns out to be less pronounced. On the contrary, the effect of the priors becomes more pronounced. Again, while differences in scores tend to be muted for models with fewer factors, the benefit of shrinkage grows when $r$ gets larger.
While joint models with $r > 0$ outperform the marginal (no-factor) model for all points in time, days of particular turbulence particularly stand out. To illustrate this, we display average and top three log predictive gains over the no-factor model in Table (ref). The biggest gains can be seen on “Black Monday 2011” (August 8), when US stock markets tumbled after a credit rating downgrade of US sovereign debt by Standard and Poor's. The trading day before this, August 4, displays the third highest gain. The 27th of February in 2007 also proves to be an interesting date to consider. This day corresponds to the burst of the Chinese stock bubble that led to a major crash in Chinese stock markets, causing a severe decline in equity markets worldwide. It appears that joint modeling of stock prices is particularly important on days of extreme events when conditional correlations are often higher.
An alternative to estimating latent factors from data is to use observed factors instead wan-etal:dyn. To explore this route, we investigate an alternative model with four observed factors, the three Fama-French plus the Momentum factor.\footnote{The Fama-French+Momentum factors are available at a daily frequency at Kenneth French's web page at \url{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}. The data was downloaded on February 21, 2017; missing values were replaced with zeroes and the data was standardized to have unconditional mean zero and variance one.} The Fama-French factors consist of the excess return on the market, a size factor SMB, and a book-to-market factor HML fam-fre:com; the momentum factor MOM car:on captures the empirically observed tendency for falling asset prices to fall further, and rising prices to keep rising. For estimation of this model we proceed exactly as before, except that we omit the last step of our posterior sampler and keep $\bm{f}$ fixed at the observed values.
Without presenting qualitative results in detail due to space constraints, we note that both the loadings on and the volatilities of the market excess returns show a remarkably close resemblance to those corresponding to the first latent factor displayed in Figures (ref) and (ref). To a certain extent (although much less pronounced), this is also true for SMB and the second latent factor as well as HML and the fourth latent factor. However, most of the loadings on MOM are shrunk towards zero and there is no recognizable similarity to the remaining latent factor from the original model. Log predictive scores for this model are very close to those for the SV model with two latent factors. In what follows, we term this approach FF+MOM.
We now turn to investigating the statistical performance of the factor SV model via out-of-sample predictive measures as well as its suitability for optimal asset allocation. The competitors are: Moving averages (MAs) of sample covariance matrices over a window of 500 trading days, exponentially weighted moving averages (EWMAs) of sample covariances defined by $ \bm{\Sigma}_{t+1} = (1-\alpha)\bm{y}_t' \bm{y}_t + \alpha \bm{\Sigma}_{t}, $ the Ledoit-Wolf shrinkage estimator led-wol:hon and FF+MOM described above. While this choice is certainly not exhaustive, it includes many of the approaches most widely used in practice.
For comparison, we use two benchmarking methods. First, we consider the minimum variance portfolio implied by $\bm{\hat\Sigma}_{t+1}$ (the point estimate or posterior mean estimate, respectively) which uniquely defines the optimal portfolio weights $$ \bm{\omega}_{t+1} = \frac{\bm{\hat\Sigma}_{t+1}^{-1}\bm{\iota}}{\bm{\iota}'\bm{\hat\Sigma}_{t+1}^{-1}\bm{\iota}}, $$ where $\bm{\iota}$ denotes an $m$-variate vector of ones. Using these weights, we compute the corresponding realized portfolio returns $r_{t+1}$ for $t = 3000,3001,\dots,3999$, effectively covering an evaluation period from 5/3/2006 to 3/1/2010. In the first three columns of Table (ref), we report annualized empirical standard deviations, annualized average excess returns over those obtained from the equal weight portfolio, and the quotient of these two measures, the Sharpe ratio sha:mut.
Considering the portfolio standard deviation presented in the first column of Table (ref), it turns out that the Ledoit-Wolf shrinkage estimator implies an annualized standard deviation of about $12.5$ which is only matched by factor SV models with many factors. Lower-dimensional factor SV models, including FF+MOM, as well as simple MAs and highly persistent EWMAs do not perform quite as well but are typically well below 20. Less persistent EWMAs, the no-factor SV model and the na\"{i}ve equal weight portfolio exhibit standard deviations higher than 20. Considering average returns, FF+MOM and factor SV models with around 10 to 20 factors tend to do well for the given time span. In the third column, we list Sharpe ratios, where factor SV models with 10 to 20 factors show superior performance, in particular when the Normal-Gamma prior is employed. Note however that column two and three have to be interpreted with some care, as average asset and portfolio returns generally have a high standard error.
Second, we use what we coin pseudo log predictive scores (PLPSs), i.e.\ Gaussian approximations to the actual log predictive scores. This simplification is necessary because most of the above-mentioned methods only deliver point estimates of the forecast covariance matrix and it is not clear how to properly account for estimation uncertainty. Moreover, the PLPS is simpler to evaluate as there is no need to numerically solve a high-dimensional integral. Consequently, it is frequently used instead of the actual LPS in high dimensions while still allowing for evaluation of the covariance accuracy ado-etal:for,car-etal:com,hub:den. More specifically, we use data up to time $t$ to determine a point estimate $\bm{\hat\Sigma}_{t+1}$ for $\bm{\Sigma}_{t+1}$ and compute the logarithm of the multivariate Gaussian density $\mathcal{N}_{m}\!\left({\bm{{0}}}, \bm{\hat\Sigma}_{t+1}\right)$ evaluated at the actually observed value $\bm{y}^o_{t+1}$ to obtain the one-day-ahead PLPS for time $t + 1$.
In terms of average PLPSs (the last column in Table (ref)), factor SV models clearly outperform all other models under consideration. In particular, even if $r$ is chosen as small as $r=4$ to match the number of factors, the model with latent factors outperforms FF+MOM. Using latent factors is generally preferable; note however that the 4-factor FF+MOM does better than the single- and no-factor SV models. Generally speaking, many factors appear to be needed for accurately representing the underlying data structure, irrespectively of the prior choice. Considering the computational simplicity of the Ledoit-Wolf estimator, its prediction accuracy is quite remarkable. It clearly outperforms the no-factor SV model which, in turn, beats simple MAs and EWMAs.
The aim of this paper was to present an efficient and parsimonious method of estimating high-dimensional time-varying covariance matrices through factor stochastic volatility models. We did so by proposing an efficient Bayesian MCMC algorithm that incorporates parsimony by modeling the covariance structure through common latent factors which themselves follow univariate SV processes. Moreover, we added additional sparsity by utilizing a hierarchical shrinkage prior, the Normal-Gamma prior, on the factor loadings. We showed the effectiveness of our approach through simulation studies and illustrated the effect of different shrinkage specifications. We applied the algorithm to a high-dimensional data set consisting of stock returns of $300$ S&P 500 members and conducted an out-of-sample predictive study to compare different prior settings and investigate the choice of the number of factors. Moreover, we discussed the out-of-sample performance of a minimum variance portfolio constructed from the model-implied weights and related it to a number of competitors often used in practice.
Because the algorithm scales linearly in both the series length $T$ as well as the number of component series $m$, applying it to even higher dimensions is straightforward. We have experimented with simulated data in thousands of dimensions for thousands of points in time and successfully recaptured the time-varying covariance matrix.
Further research could be directed towards incorporating prior knowledge into building the hierarchical structure of the Normal-Gamma prior, e.g. by choosing the global shrinkage parameters according to industry sectors. Alternatively, vil-etal:reg propose a mixture of experts model to cater for smoothly changing regression densities. It might be fruitful to adopt this idea in the context of covariance matrix estimation by including either observed (Fama-French) or latent factors as predictors and allowing for other mixture types than the ones discussed there. While not being the focus of this work, it is easy to extend the proposed method by exploiting the modular nature of Markov chain Monte Carlo methods. In particular, it is straightforward to combine it with mean models such as (sparse) vector autoregressions ban-etal:lar, kas-hub:spa, dynamic regressions kor:hie, or time-varying parameter models koo-kor:lar, hub-etal:new.
Variants of this paper were presented at the 6th European Seminar on Bayesian Econometrics (ESOBE 2015), the 2015 NBER-NSF Time Series Conference, the 2nd Vienna Workshop on High-Dimensional Time Series in Macroeconomics and Finance 2015, the 2015 International Work-Conference on Time Series Analysis, the 2016 ISBA World Meeting, the 10th International Conference on Computational and Financial Econometrics 2016, the CORE Econometrics & Finance Seminar 2017 and the 61st World Statistics Congress 2017. The author thanks all participants, in particular Luc Bauwens, Manfred Deistler, John Geweke, Florian Huber, Sylvia Kaufmann, Hedibert Lopes, Ruey Tsay, Stefan Voigt, Mike West, as well as the handling editor Herman van Dijk and two anonymous referees for crucially valuable comments and suggestions. Special thanks go to Mark Jensen who discussed this paper at the ESOBE 2015 and Sylvia Fr\"uhwirth-Schnatter who continuously supported the author throughout the development of this work.
\singlespacing \printbibliography[heading = bibintoc]