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.
69,328 characters · 16 sections · 49 citation commands
Efficient Bayesian Inference for Multivariate Factor Stochastic Volatility Models
\def\spacingset#1{ {#1}} \spacingset{1}
{\it Keywords:} Ancillarity-sufficiency interweaving strategy (ASIS), Curse of dimensionality, Data augmentation, Dynamic correlation, Dynamic covariance, Exchange rate data, Markov chain Monte Carlo (MCMC)
\thispagestyle{empty}
\spacingset{1.45}
The analysis of multivariate time series has become a vivid research area over the last decades, where both methodological as well as computational advances have made it possible to estimate more and more complex models. In parallel, real-world applications with an ever-increasing amount of data call for the joint modeling of many simultaneous and often co-varying observations over time. However, already the number of pair-wise co-movements increases quadratically with the number of time series, let alone higher-dimensional dependency structures. This property, often referred to as the curse of dimensionality, can often be mitigated in various ways by imposing a lower-dimensional latent factor structure, thereby effectively reducing the number of parameters to a feasible amount. In the paper at hand, we particularly focus on the case where these factors are allowed to have time-varying variances which in turn drive the multivariate dynamics. To the best of our knowledge, models of this type have first been discussed by jac-etal:bayJBES, she:sta, and kim-etal:sto. We particularly focus on the model formulation brought forward by chi-etal:ana.
Applications of multivariate factor stochastic volatility models typically reside in the field of financial econometrics, most prominently in areas that involve accurate quantification of uncertainty and risk. Examples thereof are asset allocation agu-wes:bay, han:ass, zho-etal:bay and asset pricing nar-scr:bay. These models extend standard factor pricing models such as the arbitrage pricing theory ros:arb and the capital asset pricing model sha:cap,lin:val by relaxing the assumption that the multivariate volatility dynamics are constant over time.
Statistical estimation of these models can be challenging, and a variety of solutions such as quasi-maximum likelihood har-etal:mul or simulated maximum likelihood lie-ric:cla,jun-koo:msv have been proposed. For medium to high dimensional problems, Bayesian MCMC estimation pit-she:tim, agu-wes:bay, chi-etal:ana, han:ass, omo-etal:sto is probably the most efficient estimation method, however, it is associated with a considerable computational burden when the number of assets is moderate to large.
The aim of this work is to outline a reliable method for Bayesian inference that performs well for a wide range of data sets while at the same time being easy to implement and convenient to extend. Therefore, we combine an efficient method for estimating univariate stochastic volatility models introduced by kas-fru:anc with a standard Gibbs sampler for regression problems. To ensure fast convergence and proper mixing of the MCMC chains we augment this simple procedure with interweaving strategies introduced by yu-men:cen. Through extensive simulation studies and a real-world example, we demonstrate the effectiveness of our procedure which can boost sampling efficiency by a factor of 100 and more.
The remainder of this paper is structured as follows. Section (ref) establishes notation for the factor stochastic volatility model framework and discusses questions about model specification and identification. Section (ref) gives an in-depth exposure to the estimation algorithm and its implementation, whereby the focus is placed on the novel interweaving strategies employed. Section (ref) presents measures of sampling efficiency for simulated data sets and compares the algorithms presented. Section (ref) discusses a case study with $26$ daily EUR exchange rates. Section (ref) concludes.
In a multivariate framework, the quadratic growth of the number of covariances alongside their inherent time-variability calls for a model which is sufficiently parsimoniously specified. At the same time, the model needs to be flexible enough to have the potential to capture typical features of financial and economic time series such as volatility clustering and volatility co-movement. On top of that, common irregularities in the data require the model to be robust with respect to idiosyncratic shocks.
The multivariate factor stochastic volatility (SV) model chi-etal:ana aims at uniting simplicity with flexibility and robustness. It is simple in the sense that the potentially high-dimensional observation space is reduced to a lower-dimensional orthogonal latent factor space, just like in the case of the classic factor model. It is flexible in the sense that these factors are allowed to exhibit volatility clustering, and it is robust in the sense that idiosyncratic deviations are themselves stochastic volatility processes, thereby allowing for the degree of volatility co-movement to be time-varying.
For each point in time $t=1,\ldots,T$, let ${\bm{ y}}_t=(y_{1t}, \ldots, y_{m t}) ^{'}$ be a zero-mean vector of $m$ observed returns and let ${\bm{ f}}_t=(f_{1t}, \ldots, f_{r t}) ^{'}$ be a vector of $r$ unobserved latent factors. In analogy to the static factor model, the observations are assumed to be driven by the latent factors and the idiosyncratic innovations. In the case of the factor stochastic volatility model, however, both the idiosyncratic innovations as well as the latent factors are allowed to have time-varying variances, depending on $m+r$ latent volatilities ${\bm{h}}_t=({\bm{h}}_t^U,{\bm{h}}_t^V)$, where ${\bm{h}}_t^U= (h_{1t}, \ldots, h_{m t}) ^{'}$ and ${\bm{h}}_t^V= (h_{m+1,t}, \ldots, h_{m+r, t}) ^{'}$. In short, we have
where $\bm{\Lambda}$ is an unknown $m\times r$ factor loadings matrix, ${\bm{{{U}}}}_t ({\bm{h}}^U_t)=\mbox{\text diag}\!\left(\exp(h_{1t}),\ldots,\exp(h_{m t})\right)$ is a diagonal $m \times m$ matrix containing the idiosyncratic (series-specific) variances, and ${\bm{ V}}_t ({\bm{h}}^V_t)=\mbox{\text diag}\!\left(\exp(h_{m+1,t}),\ldots,\exp(h_{m+r, t})\right)$ is a diagonal $r \times r$ matrix containing the factor variances. These variances are themselves modeled as latent variables whose logarithms follow independent autoregressive processes of order one, i.e.\ for $i=1,\ldots,m+r$:
with unknown initial value $h_{i0}$.
All innovations are assumed to follow independent standard normal distributions, i.e.\ $\bm{\epsilon}_t \sim \mathcal{N}_{m}\!\left({\bm{{0}}},{{\bm{I}}}_{m}\right)$, $\bm{\zeta}_t \sim \mathcal{N}_{r}\!\left({\bm{{0}}},{{\bm{I}}}_{r}\right)$, and $\bm{\eta}_{t} \sim \mathcal{N}_{m + r}\!\left({\bm{{0}}},{{\bm{I}}}_{m+r}\right)$, where $\bm{\eta}_{t} = (\eta_{1t}, \dots, \eta_{m+r,t})'$. This implies following structure:
with $\bm{\varepsilon}_t|{\bm{h}}_t \sim \mathcal{N}_{m}\!\left({\bm{{0}}},{\bm{{{U}}}}_t ({\bm{h}}_t)\right)$. One of the main reasons for estimating a factor SV model is to reliably estimate the potentially time-varying conditional covariance matrix of ${\bm{ y}}_t$ which, for the model at hand, is given by $\text{cov}({\bm{ y}}_t|{\bm{h}}_t) = {\bm{\Sigma}}_t({\bm{h}}_t) = \bm{\Lambda} {\bm{ V}}_t({\bm{h}}^V_t) \bm{\Lambda}' + {\bm{{{U}}}}_t({\bm{h}}^U_t)$. Note that because ${\bm{{{U}}}}_t ({\bm{h}}^U_t)$ is diagonal, all covariances between the component series are governed by the latent factors. Marginally with respect to ${\bm{h}}_t$, ${\bm{ y}}_t$ is a process with non-Gaussian stationary distribution.
Whenever certain combinations of parameter values result in (almost) identical maxima in the likelihood function, estimation of the corresponding parameter values from data can become impossible. Consequently, observationally equivalent parameter constellations must be ruled out for reliable statistical inference and a large body of literature dealing with this issue has arisen. In particular, fru-lop:par give an overview of recent advances in the context of static Bayesian factor models, and sen-fio:ide specifically discuss identification for models where the factors exhibit conditional heteroscedasticity (but innovations are assumed to be homoscedastic). If not dealt with properly, usually through certain restrictions on the parameter space, sensible interpretation of the posterior distribution is not possible (“nonidentifiability”). In less severe cases (“near-nonidentifiability”), MCMC algorithms and other estimation procedures often lack convergence and thus provide unreliable results. For the model at hand, we face several issues related to this problem.
First, to prevent factor rotation and column switching, one option is to follow the usual convention and set the upper triangular part of $\bm{\Lambda}$ to zero and $\mbox{\text diag}\!\left(\bm{\Lambda}\right)$ nonzero gew-zho:mea. Doing so, however, imposes an -- often unwanted -- order dependence. We therefore also discuss the possibility to leave the factor loadings matrix unrestricted and deal with column switching through post-processing of the MCMC draws.
Second, without identifying the scaling of either the $j$th column of $\bm{\Lambda}$ or the variance of $f_{jt}$, the model is not identified. The usual remedy agu-wes:bay, chi-etal:ana, han:ass, lop-car:fac, nak-wes:dyn, zho-etal:bay is that the diagonal loading elements in model ((ref)) are fixed to one, i.e.\ $\Lambda_{jj}=1$, for $j=1, \ldots, r$, while the level $\mu_{m+j}$ of the factor volatilities $h_{m+j,t}$ in model ((ref)) (which corresponds to the scaling of $f_{jt}$) is modeled to be unknown. This approach implies that the first $r$ variables are leading the factors and thus makes variable ordering an even more important modeling decision. To alleviate this issue, we leave the diagonal elements $\Lambda_{jj}$ in model ((ref)) unrestricted, an intuitive interpretation being that “leadership” of a factor can be shared by several series. Instead, we fix the level $\mu_{m+j}$ of the factor volatilities $h_{m+j,t}$ at zero:
This assumption, alongside the prior distribution on the loadings introduced in Section (ref), identifies the factor variance.
Finally, each column of $\bm{\Lambda}$ is only identified up to a possible sign switch. We deal with this (lightweight) identification issue a posteriori, meaning that we run our MCMC sampler in the unrestricted model and identify signs afterwards, see Section (ref) in Appendix (ref) for details.
Factor model ((ref)) together with the $m+r$ SV models ((ref)) defines our baseline parameterization, however alternative parameterizations will be exploited in Section (ref) in the context of efficient MCMC estimation of the factor SV model.
We perform Bayesian inference based on a set of carefully selected proper priors which are introduced in Section (ref) and develop efficient schemes for full conditional MCMC sampling in the remaining subsections.
Independently for each $i \in \{1, \dots, m+r\}$, priors for the univariate SV processes are chosen as in kas-fru:anc: $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 normal prior $\mu_i \sim \mathcal{N}\!\left(b_{\mu}, B_{\mu}\right)$, the persistence parameter $\phi_i \in (-1,1)$ is chosen according to $(\phi_i+1)/2 \sim \mathcal{B}\!\left(a_0, b_0\right)$ as in kim-etal:sto, and the volatility of log variance $\sigma_i \in \mathbb{R}^+$ is implied by $\sigma_i^2 \sim B_\sigma\times \chi^2_1= \mathcal{G}\!\left(\frac{1}{2},\frac{1}{2B_\sigma}\right)$. The initial state $h_{i0}$ is distributed according to the stationary distribution of the AR($1$) process ((ref)), i.e.\ $h_{i0}|\mu_i,\phi_i,\sigma_i \sim \mathcal{N}\!\left(\mu_i,\sigma_i^2/(1-\phi_i^2)\right)$. For every unrestricted element of the factor loadings matrix we choose independent zero-mean Gaussian distributions, i.e.\ $\Lambda_{ij} \sim \mathcal{N}\!\left(0, B_\Lambda\right)$.
Bayesian inference operates directly in the latent variable model ((ref)) and ((ref)) and relies on data augmentation by introducing the latent volatilities ${\bm{h}}=\{{\bm{h}}_{i,\bullet}\}, i=1,\ldots, m+r$, where ${\bm{h}}_{i,\bullet}=(h_{i0}, h_{i1},\dots, h_{iT})'$, and the latent factors ${\bm{ f}}=\{{\bm{f}}_{j,\bullet}\}, j=1,\ldots, r$, where ${\bm{f}}_{j,\bullet}=(f_{j1},\dots, f_{jT})'$, as latent data. This allows to set up a simple scheme for full conditional MCMC sampling which is outlined in Algorithm (ref) and discussed in detail thereafter.
For Step (a), observe that conditional on knowing the latent factors ${\bm{ f}}$ and the loadings $\bm{\Lambda}$, we are dealing with $ m+r$ independent, univariate SV models where the latent state equations ((ref)) are combined with following observation equations:
Hence, sampling the latent volatilities ${\bm{h}}_{i,\bullet}$ as well as the parameters $(\mu_i, \phi_i,\sigma_i)$ for $ i=1,\ldots,m+r$ (with $\mu_i=0$ for $i>m$) in Step (a) amounts to $m+r$ univariate SV updates. Consequently, the substantial amount of research on this matter which has emerged in the last two decades can directly be applied. In particular, we follow recent findings in kas-fru:anc, where an efficient sampling scheme is proposed and evaluated, and simply use the implementation in the \proglang{R} package {\normalfont\fontseries{b}\selectfont stochvol} kas:dea as a “plug-in” for Step (a) of the factor SV sampler presented in Algorithm (ref); see Appendix (ref) for more details and additional references on MCMC estimation for univariate SV models.
On the other hand, conditional on knowing the latent volatilities ${\bm{h}}$, we are dealing in ((ref)) with a factor model with heteroscedastic errors. Nevertheless, given ${\bm{h}}$, ${\bm{ f}}$ and $\bm{\Lambda} $ may be sampled conditionally on each other from the respective multivariate normal distributions in a similar manner as for a standard factor model lop-wes:bay. This approach is conceptually straightforward, see Appendix (ref) for details how to sample in Step (b) each row $\bm{\Lambda}_{i,\bullet}$ of the factor loading matrix from $\bm{\Lambda}_{i,\bullet}|{\bm{ f}},{\bm{ y}}_{i,\bullet},{\bm{h}}_{i,\bullet}$, where ${\bm{ y}}_{i,\bullet}=(y_{i1},\ldots, y_{iT})'$, and Appendix (ref) for details how to sample in Step (c) the factor ${\bm{ f}}_t$ from ${\bm{ f}}_t|\bm{\Lambda},{\bm{ y}}_{ t},{\bm{h}}_{ t}$ for $t=1,\dots,T$.
After discarding a certain amount of initial draws (the burn-in), the {\it standard full conditional sampler} iterates steps (a), (b) and (c) of Algorithm (ref), but not (b*), and should, in principle, yield draws from the joint posterior distribution. However, when estimating factor SV models through such an MCMC scheme, slow convergence and poor mixing (i.e.\ high correlation of posterior draws) can become a potentially prohibitive issue. This phenomenon substantiates in enormous autocorrelation of posterior draws -- even after thinning -- and can render MCMC output practically useless. For certain data sets, the burn-in phase may take extremely long and a huge amount of samples has to be discarded before the draws can be considered to emerge from the posterior distribution. Additionally, even after burn-in, these draws often show extraordinarily high autocorrelation and thus only explore the target distribution painstakingly slowly. These so-called badly mixing samplers do not only prolong computation time, they also frequently lead to unreliable estimates and misleading results. The simulation study in Section (ref) illustrates that this can happen for the standard full conditional sampler even with data simulated from the true model, see e.g.\ the top of the two panels in Figure (ref). Consequently, a carefully crafted posterior simulator is of utmost importance.
To overcome this problem, chi-etal:ana propose to sample the factor loading matrix $\bm{\Lambda}$ from the marginalized conditional posterior $p(\bm{\Lambda}|{\bm{ y}},{\bm{h}})$, without conditioning on the factors ${\bm{ f}}$. This distribution, however, is not available in closed form, and to sample from it requires a rather involved Metropolis-Hastings update where the proposal distribution is based on numerically maximizing the often high-dimensional conditional likelihood function and approximating its Hessian matrix at every MCMC iteration. To avoid this potential bottleneck, we employ the simpler full conditional procedure outlined in Algorithm (ref) but enhance it in Step (b*) by employing two variants of an ancillarity-sufficiency interweaving strategy (ASIS) yu-men:cen, called {\it shallow interweaving} and {\it deep interweaving}, which are explained in detail in Section (ref).
Applications to simulated data in Section (ref) as well as to exchange rate data in Section (ref) illustrate how adding Step (b*) boosts MCMC dramatically, in particular for deep interweaving; compare e.g.\ the top panel in Figure (ref) to the remaining panels. Section (ref) in Appendix (ref) provides comments on practical implementation of the boosted Algorithm (ref) using the \proglang{R} package {\normalfont\fontseries{b}\selectfont factorstochvol} kas:fac.
As discussed in Section (ref), the {\it standard full conditional sampler} outlined in Algorithm (ref) is based on data augmentation in the parameterization ((ref)) and ((ref)) of the factor SV model and suffers from slow convergence like so many other MCMC schemes which alternate between sampling from the full conditionals of the latent states and the model parameters. A large literature has emerged discussing various techniques to improve such algorithms, in particular reparameterization pap-etal:gen, marginal data augmentation van-men:art, and interweaving strategies yu-men:cen.
Reparameterization relies on data augmentation in a different parameterization of the model with alternative latent variables. In particular, so-called non-centered parameterizations where unknown model parameters are moved from the latent state equation to the observation equation proved to be useful, see e.g. fru-wag:sto in the context of state space modeling of time series. However, MCMC estimation based on different data augmentation schemes will often be efficient in separate regions of the parameter space, as demonstrated e.g. by kas-fru:anc in the context of univariate SV models. This suggests to combine different data augmentation schemes to obtain an improved sampler.
Marginal data augmentation employs a randomly sampled “working parameter” to transform the baseline parameterization to an expanded, unidentified latent variable model in which the model parameters are updated conditional on the (randomly) transformed latent variables. This technique has been applied to the basic factor model, using the undefined scaling of the factors as a working parameter gho-dun:def,fru-lop:par, however, it is not easily extended to factor SV models, in particular if the latent volatilities should be part of the acceleration scheme.
The ancillarity-sufficiency interweaving strategy (ASIS), introduced by yu-men:cen, provides another principled way to interweave two different data augmentation schemes by re-sampling certain parameters conditional on the latent variables in an alternative parameterization of the model, thereby combining “best of different worlds”. ASIS has been successfully employed in a variety of contexts such as univariate SV models kas-fru:anc and dynamic linear state space models sim:app,sim-etal:int. To boost Algorithm (ref), we apply ASIS to the factor SV model in the present paper. Two interweaving strategies -- called shallow interweaving and deep interweaving -- are derived in Section (ref), where the diagonal elements $\Lambda_{11}, \ldots, \Lambda_{r r}$ of the factor loadings matrix are resampled in Step (b*) in two alternative parameterizations of the model.
As will become clear in the following sections, deep interweaving typically yields the highest sampling efficiency gains and is thus the generally recommended strategy. However, also shallow interweaving has its merits. First, being conditionally conjugate, it is somewhat easier to implement. Second, it can be applied also to static factor models which are by construction not suited for deep interweaving bit-fru:ach.
As discussed in Section (ref), our baseline parameterization ((ref)) and ((ref)) is just one of several alternative ways to handle the scaling problem inherent in factor SV models and this identification issue is exploited by our schemes.
The parameterization underlying shallow interweaving constrains the diagonal elements of the factor loadings matrix to be equal to 1, whereas the variances of the factors depend on $r$ unknown scaling parameters ${\bm{ D}}=\mbox{\text diag}\!\left(\Lambda_{11}, \ldots, \Lambda_{r r}\right)$. The latent volatility processes are modeled as in the baseline parameterization ((ref)), whereas the factor model takes a different form:
with a lower triangular loading matrix $\bm{\Lambda} ^\star$ where $\Lambda ^\star_{11}=1, \ldots, \Lambda ^\star_{r r}=1$. The idiosyncratic errors $\bm{\varepsilon}_t$ are distributed as in ((ref)). Factor model ((ref)) in the baseline parameterization can be transformed into factor model ((ref)) through a simple linear transformation:
Boosting through shallow interweaving consists of three parts. First, transformation ((ref)) is used to move the current posterior draws of the latent factors ${\bm{ f}}_t$ and the factor loading matrix $\bm{\Lambda}$ from the baseline parameterization to parameterization ((ref)). Second, the scale parameters $\Lambda_{11}, \ldots, \Lambda_{r r}$, contained in ${\bm{ D}}$, are resampled in parameterization ((ref)), conditionally on the transformed values ${\bm{ f}} ^\star$, from $p(\Lambda_{11}, \ldots, \Lambda_{r r}|{\bm{ f}} ^\star, \bm{\Lambda} ^\star, {\bm{h}})$. Finally, the new values $\Lambda_{11} ^{\mbox{\rm \tiny new}}, \ldots, \Lambda_{r r}^{\mbox{\rm \tiny new}} $ are used in transformation ((ref)) to move ${\bm{ f}} ^\star_t$ and $\bm{\Lambda} ^\star$ back to new draws ${\bm{ f}}_t ^{\mbox{\rm \tiny new}}$ and $\bm{\Lambda} ^{\mbox{\rm \tiny new}}$ in the baseline parameterization.
It is evident from transformation ((ref)) that shallow interweaving only affects the factors and the factor loading matrix, whereas the latent volatilities remain untouched. This is the feature that makes shallow interweaving also applicable to static factor models. However, to achieve boosting also for the $r$ factor volatilities, deep interweaving is based on an alternative SV model for the factor volatilities where the level is assumed to be unknown. The parameterization underlying deep interweaving relies on the factor model
where $\bm{\Lambda} ^\star$ has the same structure as in the factor model ((ref)) for shallow interweaving and the idiosyncratic errors $\bm{\varepsilon}_t$ are distributed as before, with the univariate SV models for the $m$ underlying volatilities following ((ref)). However, the $r$ latent factor volatilities $h ^\star_{m+j,t}$ follow alternative univariate SV models where the level is $\mu_{m+j} =\log \Lambda_{jj}^2$ rather than zero:
This parameterization can be motivated by moving the parameters $\Lambda_{11}, \ldots, \Lambda_{r r}$ from factor model ((ref)) into SV model ((ref)), since:
Hence, the baseline parameterization can be transformed into parameterization ((ref)) and ((ref)) by applying transformation ((ref)) to the factors and the factor loadings, as well as the following transformation to the factor volatilities:
Boosting through deep interweaving also consists of three parts: first, transformations ((ref)) and ((ref)) are used to move from the current draws of ${\bm{ f}}_t$, $\bm{\Lambda}$ and the factor log-variances $h_{m+j,t}$ from the baseline parameterization to parameterization ((ref)) and ((ref)). Second, the scale parameters $\Lambda_{11}, \ldots, \Lambda_{r r}$ are resampled in parameterization ((ref)) conditionally on the transformed values ${\bm{h}} ^\star_{m+j,\bullet}=(h ^\star_{m+j,0}, \ldots, h ^\star_{m+j,T})'$ from $p(\Lambda_{11}, \ldots, \Lambda_{r r}|{\bm{h}} ^\star_{m+1,\bullet}, \cdots, {\bm{h}} ^\star_{m+r,\bullet}, \bm{\Lambda} ^\star)$. Based on the new values $\Lambda_{11}^{\mbox{\rm \tiny new}}, \ldots, \Lambda_{r r}^{\mbox{\rm \tiny new}}$, transformations ((ref)) and ((ref)) are inverted to move ${\bm{ f}} ^\star_t$, $\bm{\Lambda} ^\star$, $h ^\star_{m+j,t}$ back to new draws ${\bm{ f}}_t ^{\mbox{\rm \tiny new}}$, $\bm{\Lambda} ^{\mbox{\rm \tiny new}}$, $h_{m+j,t} ^{\mbox{\rm \tiny new}}$ in the baseline parameterization.
Both interweaving strategies are summarized in Algorithm (ref). Details on resampling $\Lambda_{jj}^{\mbox{\rm \tiny new}}$ are provided in Section (ref). It is evident that deep interweaving affects the factors, the factor loading matrix as well as the latent factor volatilities and for this reason is more effective in boosting MCMC for factor SV models than shallow interweaving.
To derive the full conditional posterior distribution of $\Lambda_{jj}$, we combine the appropriate full conditional likelihood function with the Gaussian prior $\Lambda_{jj} \sim \mathcal{N}\!\left(0, B_\Lambda\right)$. In addition, the prior $\bm{\Lambda} ^\star_{\bullet, j}|\Lambda_{jj}^2 \sim \mathcal{N}_{k_j}\!\left(0, B_\Lambda /\Lambda_{jj}^{2} {{\bm{I}}}_{k_j}\right)$ of the transformed factor loadings in column $j$ contributes to the posterior distribution of $\Lambda_{jj}^2$ because its scale depends on $\Lambda_{jj}^2$.
For shallow interweaving, we sample $\Lambda_{jj}^2$ and define $\Lambda_{jj}^{\mbox{\rm \tiny new}}$ as the square root of $\Lambda_{jj}^2$. Combining the likelihood obtained from factor model ((ref)) with the implied prior $\Lambda^2_{jj} \sim \mathcal{G}\!\left(1/2, 1/(2B_\Lambda)\right)$ and $p(\bm{\Lambda} ^\star_{\bullet, j}|\Lambda_{jj}^2)$ yields
which is the product of $T$ univariate Gaussian densities with $\Lambda_{jj}^2$ appearing as part of the variance, $k_j$ univariate Gaussian densities with $\Lambda_{jj}^2$ appearing as part of the precision, and one Gamma density with $\Lambda_{jj}^2$ appearing as argument. Thus, the resulting posterior distribution of $\Lambda_{jj}^2$ is Generalized Inverse Gaussian, i.e.
where $\text{GIG}\!\left(p,a,b\right)$ has a density proportional to $ x^{p-1}\exp\left\{-\frac{1}{2}(ax+b/x)\right\}$. Given an efficient method to draw from the GIG such as the adaptive rejection sampling algorithm provided by hoe-ley:gen, sampling from ((ref)) is straightforward. For practical implementation, we use the \proglang{R} package {\normalfont\fontseries{b}\selectfont GIGrvg} r:gig which provides a \proglang{C/C++} interface to avoid the cost of interpreting code at every MCMC iteration, thereby rendering the re-updating negligible in terms of overall computation time.
For deep interweaving, we sample $\Lambda_{jj}$ indirectly through $\mu_{m+j} =\log \Lambda_{jj}^2$. Combining the implied prior $p(\mu_{m+j}) \propto \exp\left\{\mu_{m+j}/2 - e^{\mu_{m+j}}/{(2B_\Lambda)}\right\}$ with the likelihood obtained from SV model ((ref)) and the priors $h ^\star_{m+j,0}| \mu_{m+j}, \phi_{m+j}, \sigma^2_{m+j} \sim \mathcal{N}\!\left(\mu_{m+j},\sigma_{m+j}^2/(1-\phi_{m+j}^2)\right)$ and $\bm{\Lambda} ^\star_{\bullet, j}|\mu_{m+j} \sim \mathcal{N}_{k_j}\!\left(0, B_\Lambda e^{-\mu_{m+j}} {{\bm{I}}}_{k_j}\right)$ yields the posterior
which has a non-standard form. To generate draws from this density, we consider an independence Metropolis-Hastings update in the spirit of kas-fru:anc. Since the likelihood $p(h ^\star_{m+j,1}, \ldots, h ^\star_{m+j,T} |h ^\star_{m+j,0}, \mu_{m+j}, \phi_{m+j}, \sigma^2_{m+j})$ is the kernel of a Gaussian density in $\mu_{m+j}$, it can be used to construct an auxiliary posterior under a conjugate auxiliary prior $p_\text{aux}(\mu_{m+j}|\sigma^2_{m+j},\phi_{m+j}) \sim \mathcal{N}\!\left(0,B_0\sigma^2_{m+j}/(1-\phi_{m+j})^2\right)$ with $B_0$ large. Consequently, we draw a proposal $\mu_{m+j}^\text{prop}$ from the $\mathcal{N}\!\left(m_j^\mu,S_j^\mu\right)$ distribution with:
Denoting the old value of $\mu_{m+j}$ by $\mu_{m+j}^{\mbox{\rm \tiny old}}$, this proposal gets accepted with probability $\min(1,R)$, where \[R= \frac{ p(\bm{\Lambda} ^\star_{\bullet, j}|\mu^\text{prop}_{m+j}) p(h ^\star_{m+j,0}|\mu_{m+j}^\text{prop}, \phi_{m+j}, \sigma^2_{m+j}) p(\mu_{m+j}^\text{prop})} { p(\bm{\Lambda} ^\star_{\bullet, j}| \mu^{\mbox{\rm \tiny old}}_{m+j}) p(h ^\star_{m+j,0}|\mu_{m+j}^{\mbox{\rm \tiny old}}, \phi_{m+j}, \sigma^2_{m+j}) p(\mu_{m+j}^{\mbox{\rm \tiny old}}) } \times \frac{p_\text{aux}(\mu_{m+j}^{\mbox{\rm \tiny old}}|\sigma^2_{m+j},\phi_{m+j})} {p_\text{aux}(\mu_{m+j}^\text{prop}|\sigma^2_{m+j},\phi_{m+j})}.\] In case of acceptance, set $\Lambda_{jj}^\text{new} = e^{\mu_{m+j}^\text{prop}/2}$; otherwise, let $\Lambda_{jj}^\text{new} = \Lambda_{jj}^\text{old}$.
To conclude, we add two remarks. First, note that it is easy to combine both interweaving schemes within the MCMC sampler by daisy-chaining the corresponding steps. Second, note that any nonzero element of the $j$th factor column $\bm{\Lambda}_{\bullet, j}$ can be used to boost its mixing (not only the diagonal element $\Lambda_{jj}$). This is useful in particular when no loading matrix restrictions are enforced as it cannot be guaranteed that the diagonal elements are nonzero. Thus, in such situations, one could use a randomly selected (nonzero) element of each loadings column instead. Alternatively, one could also use the element whose absolute value is maximal.
In order to compare the different algorithms in terms of sampling efficiency, a simple simulation experiment is conducted. We use $m=10$ (simulated) series and $r=2$ (simulated) factors to generate $T=1000$ observations, thereby imposing the usual lower triangular constraint. The data generating parameter values -- listed in Table (ref) in Appendix (ref) -- are kept constant, whereas the data generating process as well as the estimation procedure is repeated $100$ times. Each time, the draws are initialized at the data generating values; then, $5100\,000$ draws are obtained of which $100\,000$ are discarded as burn-in. Prior hyperparameters are set as follows: $B_\Lambda = 1$, $b_\mu = 0$, $B_\mu = 100$, $a_0 = 20$, $b_0 = 1.5$, and $B_\sigma = 1$.
To gain insight about the mixing behavior of the different sampling strategies, trace plots (i.e.\ time series plots of the MCMC draws) for $\Lambda_{11}$ are displayed in the left hand panel of Figure (ref). Even though the plots depict only the first $10\,000$ iterations after burn-in, it becomes very clear that the mixing of the non-interwoven sampler is extremely slow. The algorithm doesn't seem to explore the posterior distribution within a reasonable amount of draws which renders this output practically useless in terms of posterior inference. Moreover, the burn-in period for this sampler would need to be chosen extremely long to avoid strong dependence on the starting values. This situation is slightly mitigated when using shallow interweaving; nevertheless, mixing is still poor and for reliable posterior inference many draws are required. Turning towards the deeply interwoven sampler, one can observe quick mixing and hardly any visible autocorrelation.
Investigating autocorrelations of the draws via the empirical autocorrelation function confirms this picture; the right hand panel of Figure (ref) shows that the empirical autocorrelation function for draws from $p(\Lambda_{11}|\bm{y})$ decays very quickly for the sampler using deep interweaving which is not the case for the other two samplers, where visible autocorrelation remains even at large lags.
A convenient and common way of measuring sampling (in)efficiency is by means of the inefficiency factor (IF), sometimes called the (integrated) autocorrelation time. It is defined as the ratio of the numerical variance of a statistic which is estimated from the Markov chain to the variance of that statistic when estimated from independent draws, thereby quantifying the relative loss of efficiency when inferring from correlated as opposed to independent samples. In other words, to achieve the same inferential accuracy about some posterior moment of some parameter as with $k$ independent samples, $\text{IF}\times k$ MCMC draws are required. For the paper at hand, we use the \proglang{R} package {\normalfont\fontseries{b}\selectfont coda} plu-etal:cod to estimate the inefficiency factors.
Moreover, when investigating performance of MCMC samplers through simulation studies, it is of great importance to take sample variation into account; even when identical parameter values are used for the generation of latent variables and data, sampling (in)efficiency may vary greatly. To illustrate this, we show box plots of the inefficiency factors stemming from repeated data generating processes in the left panel of Figure (ref). Note the enormous range for the standard sampler; depending on the data, IFs of $5000$ or more are not uncommon, while at the same time IFs of around $100$ can be observed. However, independently of the actual data, interweaving attenuates this effect drastically and increases efficiency uniformly. The right panel of Figure (ref) shows pairwise scatter plots of these IFs. Note that shallow interweaving yields efficiency improvements which are more or less independent of the actual data (around five-fold for all data sets), whereas deep interweaving IFs appear less clearly correlated.
To provide a more complete picture, we list the inefficiency factors for all elements of $\bm{\Lambda}$ for the various algorithms in Table (ref), averaged over all $100$ runs. Note that shallow interweaving permits efficiency gains of around two- to eight-fold as opposed to the standard sampler, whereas deep interweaving delivers gains up to around $400$-fold.
It comes as no surprise that sampling (in)efficiencies of draws for the volatility parameters $\mu_i$, $i \in \{1,\dots,m\}$ as well as $\phi_i$ and $\sigma_i$, $i \in \{1,\dots,m+r\}$ are not affected substantially by this interweaving strategy, thus they are not reported here. It is however worth noting that the inefficiency of factor ${\bm{f}}_{j,\bullet}$ and factor log-variance draws ${\bm{h}}_{m+j,\bullet}$, $j \in \{1,\dots,r\}$ may be influenced by bad mixing of $\bm{\Lambda}$. For illustration, IFs are reported for the final factors $f_{1T}$ and $f_{2T}$ and their log-variances $h_{m+1,T}$ and $h_{m+2,T}$ in Table (ref).
To conclude the simulation exercise, we investigate predictive performance of under- and overfitting models through cumulative log predictive Bayes factors in Figure (ref); see kas:spa for computational details. It stands out that the biggest predictive gain over a model that ignores contemporaneous correlations comes from introducing the first factor, i.e. allowing for co-volatility through one common factor. Then, as expected, the second factor bumps the predictive score to its maximum. After that, it remains (almost) constant for three and more factors, irrespectively of whether the lower triangular restriction is enforced or not. This points out that underfitting models are severly worse in terms of prediction while overfitting models hardly suffer from the extra parameters introduced.
Note that, similar to a scree plot in principal component analysis, Figure (ref) can also be used as a graphical tool for finding the appropriate number of factors. For this exercise, we clearly find the true number of two factors.
In this section, we analyze exchange rates with respect to EUR. Data was obtained from the European Central Bank's Statistical Data Warehouse and ranges from April 1, 2005 to August 6, 2015. It contains $m=26$ (all which were available for this time frame) daily exchange rates on $2650$ days listed in Table (ref). For further analysis, we thus use $T=2649$ demeaned log returns. The data is displayed in Figure (ref) in Appendix (ref). Common “stylized facts” of financial time series are clearly visible; note e.g.\ the obvious volatility clustering during 2008 and 2009 and again throughout late 2014 and early 2015. To put the robustness of our sampler to the test, we use the data as-is, i.e.\ without excluding series containing extreme outliers such as the CHF spike on January 14, 2015 or the near collapse of RUB around December 16, 2014.
For selecting the number of factors in this application, it is important to keep in mind the primary purpose of the analysis. Different sampling strategies are applied, depending on whether identification of $\bm{\Lambda}$ is of no concern (e.g.\ for covariance matrix prediction only), or whether identification is instrumental for understanding the unobserved, underlying factors.
For the first case, we experimented with fitting unrestricted models to the exchange rates data. This implies that the method is completely invariant to series ordering and there are no model-implied “leading factors” as is usually the case. With respect to selecting the number of factors, we found that higher-order models without any restriction on the factor loadings matrix yield higher marginal likelihoods and are thus recommended. Figure (ref) illustrates this via log predictive Bayes factors.
If identification is warranted, an important step is the appropriate ordering of the variables, before the usual lower triangular structure is imposed on the factor loadings matrix to guarantee mathematical identifiability, as outlined in Section (ref). This, however, makes inference on the factor loadings matrix dependent on the appropriate ordering of the variables. We exemplify this by the predictive Bayes factors for models where the component series are ordered alphabetically and the lower triangular structure is imposed on the first three series appearing in Table (ref), namely AUD, CAD and CHF. As shown in Figure (ref), for a given number of factors $r$, the log predictive Bayes factors of the constrained models are consistently smaller than for the unrestricted models, indicating that a purely mathematical identifiability constraint may be in conflict with the data.
Also for the constrained models, the log predictive Bayes factors are ever increasing for the exchange rate data. However, it stands out that the relative gain per additional factor is highest for few factors, flattening out quickly. Furthermore, draws of the factor loading matrix in these higher-order models are difficult to identify, in particular, if the lower triangular constraint is in conflict with the data and spurious factors, i.e.\ factors which are significantly loaded on by only few series, are present fru-lop:par. Thus, to keep presentation feasible and to avoid spurious factors, we restrict ourselves to a model with $r=4$ factors for the following in-depth discussion.
One way to find constraints that are not in conflict with the data is to post-process the MCMC draws of an unrestricted sampler with $r=4$ factors, a method that has been applied in con-etal:bay and ass-etal:bay. Note that rather than reordering the variables before imposing a lower triangular constraint, we can choose three (out of the 26) currencies, and impose the $r(r-1)/2=6$ zero restrictions on the corresponding factor loadings, see e.g.\ dun:not. While the choice of these currencies is not unique, inference is robust to specific choices, as long as the corresponding currencies serve as “leaders” for specific factors. This is exemplified by Figure (ref) which displays the posterior median of the MCMC draws of the factor loadings obtained from an unrestricted sampler with $r=4$. To solve column switching, the columns of $\bm{\Lambda}$ are rearranged by the size of their maximum median loading. According to Figure (ref) in Appendix (ref), USD is a definite candidate to lead factor one, PLN leads a second factor, and AUD leads a third factor. Alternatively, identification could be based on any other currency strongly loading on factor 1 (such as HKD or CNY), in combination with HUF (instead of PLN) and NZD (instead of AUD).
Prior hyperparameters are the same as for the simulation study in Section (ref). A sensitivity analysis shows that none of the hyperparameter choices turn out to be very influential in this particular application with the exception of the prior factor loadings variances $B_\Lambda$. These, however, are only important for the absolute scaling of the factors and do not notably influence the relative loadings sizes or predictive Bayes factors. We run each sampler for $550\,000$ iterations, discard the first $50\,000$ draws as burn-in, leaving $500\,000$ which we use for posterior inference. Even after this substantial amount of iterations, it is not clear that the sampler without interweaving has properly converged; we therefore omit its presentation. IFs from the interwoven samplers are presented in Table (ref) in Appendix (ref).
Finally, we identify the signs of the loadings in the post-processing phase by investigating the MCMC draws. For each factor, the series whose posterior absolute loadings distribution is furthest away from zero is assigned a positive sign, the other loadings are aligned thereafter, see also Section (ref) in Appendix (ref).
We begin by discussing the log-variances of the latent factors, visualized in Figure (ref), alongside the corresponding factor loadings whose marginal posterior distributions are depicted in Figure (ref) and whose posterior means are listed in Table (ref).
The first factor can clearly be interpreted as the USD-driven one, as the pegged triplet USD, CNY and HKD loads very highly on this factor, alongside many other currencies. Its volatility is generally very smooth, rising in the aftermath of the 2008 financial crisis and going down again after 2009; a second increase can be seen in the second half of 2014, possibly in connection with the Greek government-debt crisis. Factor 2's log-variance appears slightly less persistent and more volatile, it is driven by ZAR, the only African currency in the sample, alongside Eastern Europe's / Southwestern Asia's HUF, PLN and TRY. Interestingly, JPY loads negatively on this factor. The third factor shows a similar overall pattern as the first. The highest loading series for this factor are AUD and NZD, emphasizing the Trans-Tasman relations. Other commodity currencies such as ZAR and CAD also load highly on this factor. Factor 4 is clearly driven by the currencies of the Tiger Cub economies such as MYR, KRW, PHP and SGD.
In order not to overload the graphical displays used to visualize the results of the analysis, we display the results for a two-year period only for the rest of this section. More specifically, we look at the years 2008 and 2009, covering the most volatile span during the financial crisis. Irrespectively of that, the full data set has been used for estimation and other time spans could be displayed analogously.
We start out by visualizing the marginal posterior means of univariate volatilities for all $26$ currencies in Figure (ref) from the last day of 2007 until the last day of 2009. Series such as DKK or HRK are (very) closely pegged to EUR and unsurprisingly show very low volatility throughout the crisis. Other European currencies (CHF\footnote{It is interesting to note that the Swiss franc stays comparably stable from a EUR perspective throughout 2007-2009. Very differently during summer 2011, where CHF shows atypical and very high volatility until the Swiss Central Bank sets the minimum exchange rate at CHF 1.20 per EUR 1 on September 6.}, RON, SEK, CZK) follow suit. Tiger Cub economies such as PHP, HKD, THB, and MYR align very closely with USD and CNY. The most volatile currencies during this period are KRW, ZAR, IDR, JPY, and also TRY, followed by NZD, AUD, and CAD. Overall, it stands out that even though some series-specific ups and downs can be spotted, a common trend is clearly visible.
Next, implied correlation matrices are displayed in Figure (ref), exemplified for the last day of 2007, 2008, and 2009. Additionally to displaying the mean posterior pairwise correlations (at the given dates) via color and shading, these plots visualize posterior uncertainty; the outer and inner circles' sizes correspond to $\text{posterior mean} \pm 2\text{ standard deviations}$, respectively. The images were generated using the \proglang{R} package {\normalfont\fontseries{b}\selectfont corrplot} r:corrplot; its option \code{hclust} (hierarchical clustering) was used for ordering the series to emphasize the blocks of currencies.
To further illustrate variability over time, we determine the posterior means of the pairwise time-varying correlations of USD against the other currencies which are plotted in Figure (ref) in Appendix (ref). As was to be expected, correlations of CNY and HKD with USD are almost always very close to one; IDR, THB, SGD, MYR, and PHP show rather high correlation throughout. The correlation between USD and RUB on the other hand falls from around $0.9$ in early 2008 to around $0.4$ in late 2009, whereas THB moves in the opposite direction; its correlation with USD is around $0.5$ at the beginning of the time window and increases quickly to around $0.9$. Eastern European non-euro currencies, in particular PLN and HUF, appear to be slightly negatively correlation with USD throughout the entire period.
Estimating time-varying (dynamic) covariance and correlation matrices of financial and economic time series constitutes a current and active area of research. One of the main challenges thereby is the curse of dimensionality, i.e.\ the fact that the number of elements of these matrices grows quadratically with the number of observed series. We address this issue by imposing a low-dimensional latent factor structure where the factors are allowed to exhibit stochastic volatility and thereby govern co-movement of volatility over time. To conduct reliable statistical inference, we propose novel Bayesian MCMC algorithms which exploit the model-inherent identifiability constraints. By interweaving different (but mathematically equivalent) parameterizations, the proposed strategies substantially improve mixing of draws obtained from the posterior distribution, in particular for the factor loadings matrix. The method proposed is fully automatic in the sense that the end-user is not required to manually adjust any tuning parameters.
In an extensive case study discussing exchange rates with respect to EUR we show that the algorithm plays well with real-world data that exhibits a fair degree of outliers (e.g.\ CHF, RUB) which are captured through the idiosyncratic stochastic volatility components. The model structure allows for a covariance decomposition in four interpretable factors (USD/CNY driven, Eastern Europe, commodity currencies, Tiger Cub economies). These, alongside the idiosyncratic volatilities, drive the dynamics of the joint correlation structure. The pairwise correlations with USD range from “almost perfect” (CNY, HKD) over “hardly existent” (CHF, HRK) to “slightly negative” (PLN, HUF) with a varying and time-dependent degree of variability.
Concerning extensions of the model, we point out that due to the modular nature of MCMC, all the ideas of this paper can be straightforwardly generalized to models that independently model the mean, be it through a simple nonzero mean vector, a local level model, external regressors, or via (vector) autoregressive processes. For models where the level of the returns explicitly depends on the (co-)volatilities cha:sto or the returns are assumed to be correlated with the (co-)volatilities ish-omo:por, more involved estimation methods are required; this unfortunately places their discussion outside the scope of this paper. Nevertheless, due to the growing number of successful applications of interweaving methods in different contexts, there is good reason to hope for similar effects when they are used for these type of factor SV model extensions.