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.
116,334 characters · 33 sections · 82 citation commands
Modeling Univariate and Multivariate Stochastic Volatility in R with stochvol and factorstochvol
This preprint corresponds to the forthcoming article of the same name thisjss.
Time dependent variance is an indispensable ingredient of financial and economic time series modeling. Already markowitz1952portfolio concerns himself with methods that take into account heteroskedasticity in a better way than a rolling window estimation. By 1982, two fundamentally different approaches had been developed to cater to these needs. On the one hand, engle1982arch lays the groundwork for a family of time varying volatility models, most notably the generalized autoregressive conditional heteroskedasticity model bollerslev1986garch. These models feature conditionally deterministic changes in the variance. taylor1982sv, on the other hand, addresses heteroskedasticity in his seminal work with a non-linear latent state space model, later coined the stochastic volatility (SV) model. There, the volatility process evolves in a stochastic manner. Despite some empirical evidence in favor of SV models over their corresponding GARCH counterparts jacquier1994bayesian,ghysels1996stochastic,kim1998stochastic,nakajima2012bayesian, SV and its variants enjoy little publicity among practitioners. As bos2012relating underlines, one reason for this might be the lack of standard software. In response, kastner2016dealing provides a first version of the \proglang{R} rlanguage package \pkg{stochvol} but fails to feature conditional non-Gaussianity, asymmetry (the so-called leverage effect), and multivariate generalizations.
We address these shortcomings in the manuscript at hand. First, we extend \pkg{stochvol} kastner2020stochvol with several practically relevant univariate methods. Second, we introduce the new companion package \pkg{factorstochvol} factorstochvol which focuses on the multivariate case. The extended \pkg{stochvol} now provides the means for the Bayesian estimation of vanilla SV, heavy-tailed SV, SV with leverage, and heavy-tailed SV with leverage harvey1996estimation,omori2007stochastic,nakajima2012stochastic. Moreover, the package also handles these models naturally when embedded into a linear model or an autoregressive (AR) context. The \pkg{factorstochvol} package implements an efficient method for the Bayesian estimation of the factor SV model kastner2017efficient. Among other features, the package provides several automatic factor identification schemes, hierarchical shrinkage priors griffin2010inference, and an array of intuitive visualization methods for the high-dimensional posteriors.
The remainder of this paper is structured as follows. We formally introduce the univariate and the multivariate models in Sections (ref) and (ref), respectively, including a discussion about prior distributions and a brief overview of the estimation methods. In Section (ref), we unveil the new samplers of the \pkg{stochvol} package through three example models. We describe the \pkg{factorstochvol} package in Section (ref), and then we conclude.
We begin by introducing the vanilla SV model with linear regressors, henceforth simply called the SV model. This is a minor but important extension of the SV model without regressors. We also settle the notation and establish a baseline model that we generalize and reuse throughout the manuscript. Consequently, we proceed with three generalized models: the SV model with Student's $t$ errors (SVt), the SV model with leverage (SVl), and their combination, the SV model with Student's $t$ errors and leverage (SVtl). Finally, we close the section after discussing prior distributions and Markov chain Monte Carlo (MCMC) sampling.
The key feature of the SV model is its stochastic and time-varying specification of the variance evolution. In particular, the log-variance is assumed to follow an AR(1) process. This feature unites the following models.
Let $\bm{y}=(y_1,\dots,y_n)^\top$ denote a vector of observations. The SV model assumes the following structure for $\bm{y}$,
where $\Normal{b}{B}$ denotes the normal distribution with mean $b\in\mathbb{R}$ and variance $B\in\mathbb{R}^+$, and $\varepsilon_t$ and $\eta_t$ are independent. The log-variance process $\bm{h}=(h_1,\dots,h_n)^\top$ is initialized by $h_0\sim\Normal{\mu}{\sigma^2/(1-\varphi^2)}$. $\bm X=(\bm x_1^\top,\dots,\bm x_n^\top)^\top$ is an $n\timesK$ matrix containing in its $t$th row the vector of $K$ regressors at time $t$. The $K$ regression coefficients are collected in $\bm\beta=(\beta_1,\dots,\beta_K)^\top$. We refer to $\bm{\vartheta}=(\mu,\varphi,\sigma)$ as the SV parameters: $\mu$ is the level, $\varphi$ is the persistence, and $\sigma$ (also called volvol) is the standard deviation of the log-variance.
Several authors have suggested to use non-normal conditional residual distributions for stochastic volatility modeling. Examples include the Student's $t$ distribution harvey1994multivariate, the extended generalized inverse Gaussian silva2006extented, (semi-)parametric residuals jensen2010bayesian,delatola2011bayesian, or the generalized hyperbolic skew Student's $t$ distribution nakajima2012stochastic. We implement Student's $t$ errors for the observation equation in \pkg{stochvol}. Formally,
where $\varepsilon_t$ and $\eta_t$ are independent. $t_{\nu}(a,b)$ is the Student's $t$ distribution with $\nu$ degrees of freedom, mean $a$, and variance $b$. The single difference between Equation (ref) and Equation (ref) is that here the observations are conditionally $t$ distributed. Hence, Equation (ref) generalizes Equation (ref) through the new parameter $\nu$ as the Student's $t$ distribution converges in law to the standard normal distribution when $\nu$ goes to infinity.
Propositions for asymmetric innovations include non-parametric distributions jensen2014estimating, skewed distributions nakajima2012stochastic, and distributions featuring correlation with the variance process, also called the leverage effect harvey1996estimation,jacquier2004bayesian. We implement the leverage effect in the \pkg{stochvol} package. Formally,
where the correlation matrix of $(\varepsilon_t,\eta_t)$ is
The vector $\bm{\zeta}=(\mu,\varphi,\sigma,\rho)^\top$ collects the SV parameters. The new parameter compared to Equation (ref) is a correlation term $\rho$ which relates the residuals of the observations to the innovations of the variance process. Equation (ref) is therefore a special case of Equation (ref) with a pre-fixed $\rho=0$.
Some authors have proposed the combination of $t$ errors with the leverage effect jacquier2004bayesian,omori2007stochastic,nakajima2009leverage. We implement the common generalization of Equation (ref) and Equation (ref). Formally,
where the correlation matrix of $(\varepsilon_t,\eta_t)$ is $\bm\Sigma^\rho$ as in Equation (ref).
We a priori assume $\bm\beta\sim\Normal[K]{\bm{b_\beta}}{\bm{B_\beta}}$, where $\Normal[l]{\bm b}{\bm B}$ is the $l$-dimensional normal distribution with mean vector $\bm b$ and variance-covariance matrix $\bm B$. For small values in the diagonal of $\bm{B_\beta}$, this prior enforces shrinkage towards $\bm{b_\beta}$; for large values in the diagonal, the prior turns rather uninformative. By setting $\bm{b_\beta}$ to the zero vector and $\bm{B_\beta}$ to a scaled identity matrix, the prior distribution becomes the Bayesian analogue to ridge regression park2008bayesian.
The level $\mu\in\mathbb{R}$ is unrestricted, hence we can apply the common $\mu\sim\Normal{b_\mu}{B_\mu}$ prior. Depending on the application, a fairly uninformative distribution is the usual choice, e.g., setting $b_\mu=0$ and $B_\mu\ge100$ for daily asset log returns. In our experience, the exact values of the prior mean and prior variance of $\mu$ do not strongly affect the estimation results unless $B_\mu$ is small.
To achieve stationarity in the variance process, a restricted persistence $\varphi\in(-1,1)$ is needed. To this end, we assume $(\varphi+1)/2\sim\mathcal{B}(a_\varphi,b_\varphi)$, where $\mathcal{B}(a_\varphi,b_\varphi)$ is the beta distribution with shape parameters $a_\varphi$ and $b_\varphi$. The selection of the shape parameters may be relatively influential with many data sets. In financial applications with daily asset log returns, the variance tends to be highly persistent, i.e., $\varphi\approx1$. Such domain knowledge can be used as prior information by allocating more probability to positive high values of $\varphi$, e.g., by setting $a_\varphi\gtrsim5$ and $b_\varphi\approx1.5$. As an alternative, when stationarity is not assumed, the untruncated prior $\varphi\sim\Normal{b_\varphi}{B_\varphi}$ can also be applied.
The volvol is positive but we would like allow $\sigma$ to approach 0 as closely as needed -- that allows us to be less informative and to improve the estimates. Following fruhwirth2010stochastic and kastner2014ancillarity, we advocate $\sigma\sim\lvert\Normal{0}{B_\sigma}\rvert$ instead, where $\lvert\Normal{0}{B_\sigma}\rvert$ denotes the half normal distribution. It corresponds to $\sigma^2\sim\mathcal{G}(1/2,1/(2B_\sigma))$, where $\mathcal{G}(a,b)$ is the gamma distribution with shape parameter $a$ and rate parameter $b$. As an alternative, the commonly applied and convenient conjugate gamma prior on $\sigma^{-2}$ can be assumed. However, it bounds $\sigma$ away from 0 and it is therefore in our view an unsatisfactory choice.
As a last step in fully specifying the vanilla SV model in Equation (ref), the variance process is initialized a priori with its stationary distribution, i.e., $h_0\sim\Normal{\mu}{\sigma^2/(1-\varphi^2)}$. This consistently extends our prior assumptions about $\bm{h}$ following a stationary AR(1) process. As an alternative, when stationarity is not assumed, $h_0\sim\Normal{\mu}{B_h}$ can be applied with a constant variance $B_h$.
The SV models with Student's $t$ errors additionally require the prior specification of the degrees of freedom parameter $\nu$. To ascertain interpretability of the scaling $\exp(h_t/2)$, we ensure finite second moments of $\bm{y}$ by enforcing $\nu>2$. As a reviewer recommended, we follow geweke1993bayesian and equip $\nu$ with an exponential prior $\nu-2\sim\mathcal{E}(\lambda_\nu)$, where $\lambda_\nu$ is the rate of the exponential distribution.
Finally, in the case of the SV models with leverage, we employ the translated and scaled beta distribution for $\rho\in(-1,1)$ as in omori2007stochastic, i.e., $(\rho+1)/2\sim\mathcal{B}(a_\rho,b_\rho)$. We find that the posterior estimates of $\rho$ can be sensitive to its prior distribution, thus, some care is needed when setting the hyperparameters in practice. In our experience, slightly informative choices such as $a_\rho=b_\rho\approx4$ work well in financial applications.
All methods implemented in \pkg{stochvol} and \pkg{factorstochvol} rely on the Bayesian paradigm.\footnote{At this point we would also like to point out works aiming at estimating stochastic volatility and related models within the frequentist framework, see, e.g., abanto2017maximum, creal2017class, and in particular the recent \pkg{stochvolTMB} package stochvolTMB.} Bayesian analysis aims to estimate model parameters through Bayesian updating. By using probability distributions to represent information, Bayes' theorem can be employed to update the prior information to the posterior information by incorporating the observations. This approach has the advantage of providing full uncertainty quantification in a probabilistic framework without relying on asymptotic results; moreover, so-called shrinkage priors can be used to regularize the posterior and guard against overfitting. For an introductory textbook on Bayesian statistics, see, for instance, mcelreath2015statistical.
When the posterior distribution is not available analytically, one customarily resorts to approximations such as perfect simulation huber2015perfect, approximate Bayesian computation sisson2018handbook, adaptive Monte Carlo methods roberts2005coupling, or MCMC methods. When computationally feasible, MCMC is a valuable tool that provides draws from the posterior distribution in question. That way, MCMC approximates the posterior distribution similarly to a histogram approximating a density. For a more in-depth introduction on MCMC methods, see, for instance, brooks2011handbook.
The estimation algorithm of SV, SVt, SVl, and SVtl all resemble the original methodology developed in kastner2014ancillarity for the vanilla SV model. Namely, to draw from the posterior distribution of $\bm{h}$ efficiently, the MCMC sampler resorts to approximate mixture representations of Equations (ref), (ref), (ref) and (ref) similar to the ones in kim1998stochastic and omori2007stochastic. Doing so yields a conditionally Gaussian state space model for which efficient sampling methods are available fruhwirth1994data,carter1994gibbs. Following rue2001fast and mccausland2011simulation, we draw the full vector $\bm{h}$ “all without a loop” (AWOL).
When Student's $t$ errors with unknown degrees of freedom are used, we handle the added complication through the well-known representation of the $t$ distribution as a scale mixture of Gaussians. This requires additional Gibbs and independence Metropolis-Hastings steps documented in kastner2015heavy. Furthermore, we deal with the increased complexity in the posterior space of the leverage case by repeated ancillarity-sufficiency interweaving strategies yu2011to steps in the sampling scheme, see hosszejni2019approaches for details.
To verify the correctness of the implementation, unit tests are included in the package which can be run by \code{devtools::test()} devtools. In particular, a variant of Geweke's test geweke2004getting is part of the test suite. In this test, we exploit that the sampling distribution of the model parameters during the Geweke test is identical to their preset prior distribution. Therefore, the cumulative distribution function maps the sample to a uniform distribution, which in turn is mapped to a normal distribution using the normal distribution's quantile function. If the user chooses to execute the automated unit tests in \pkg{stochvol}, the system evaluates the thinned and transformed sample using the \code{shapiro.test()} function, where the thinning of the sample is done to approximate independent sampling.
For maximal computational effectiveness, all sampling algorithms are implemented in the compiled language \proglang{C++} cpplanguage with the help of the \proglang{R} package \pkg{Rcpp} eddelbuettel2011rcpp. Matrix computations make use of the efficient \proglang{C++} template library \pkg{Armadillo} sanderson2016armadillo through the \proglang{R} package \pkg{RcppArmadillo} eddelbuettel2014rcpparmadillo.\footnote{For explicit run time discussions please see kastner2014ancillarity and hosszejni2019approaches. For the possibility to use multi-core computing within a single MCMC chain and potential speed gains when doing so, please see kastner2019sparse.} After sampling, the resulting \proglang{R} objects make use of plotting and summary functions of the \proglang{R} package \pkg{coda} coda.
A key difficulty accompanying dynamic covariance estimation is the relatively high number of unknowns compared to the number of observations. More precisely, letting $m$ denote the cross-sectional dimension, the corresponding covariance matrix $\Sigmay_t$ contains $m(m+1)/2$ degrees of freedom, a quadratic term in $m$. Table (ref) illustrates the “curse of dimensionality” for various values of $m$. One way to break this curse is to use latent factors and thereby achieve a sparse representation of $\Sigmay_t$.
Latent factor models embody the idea that even high dimensional systems are driven by only a few sources of randomness. These few sources of randomness control a few factors, which in turn account for the interactions between the observations. Moreover, latent factor models provide an efficient tool for dynamic covariance matrix estimation. They allow for a reduction in the number of unknowns. A conventional latent factor model with $r$ factors implies the decomposition
where $\operatorname{rank}(\Sigmay[\check]_t)=r<m$, and $\Sigmaidi_t$ is the diagonal matrix containing the variances of the idiosyncratic errors. The rank assumption on the symmetric $\Sigmay[\check]_t$ gives rise to the factorization $\Sigmay[\check]_t=\Psi\Psi^\top$, where $\Psi\in\mathbb{R}^{m\timesr}$ contains $mr-r(r-1)/2$ free elements higham1990analysis. Hence, $m(r+1)-r(r-1)/2$ free elements remain in $\Sigmay_t$, now only linear in $m$. Table (ref) illustrates the “broken curse of dimensionality” for various values of $m$ and $r=4$.
In the following, we describe the factor SV model employed in the \pkg{factorstochvol} package. We model the observations $\bm{y}_t=(y_{t1},\dots,y_{tm})^\top$ as follows.
where $\bm{f}_t=(f_{t1},\dots,f_{tr})^\top$ is the vector of factors, $\bm\beta=(\beta_1,\dots,\beta_m)^\top$ is an observation-specific mean, and $\bm{\Lambda}\in\mathbb{R}^{m\timesr}$ is a tall matrix holding the factor loadings. The covariance matrices $\Sigmaidi_t$ and $\Sigmafac_t$ are both diagonal representing independent vanilla SV processes.
For a more theoretical treatment of factor SV from a Bayesian point of view, the reader is referred to, e.g., pitt1999time, aguilar2000bayesian, chib2006analysis, and han2006asset.
Based on Equation (ref), we can reformulate Equation (ref) as
from which several identification issues are apparent: the order, the sign, and the scale of the factors is unidentified. More specifically, for any generalized permutation matrix \footnote{A generalized permutation matrix has the zero--non-zero pattern of a permutation matrix, but it is allowed to have any non-zero values instead of just ones. Hence, a generalized permutation matrix not only permutes but also scales and switches the sign of its multiplier.} $\bm{P}$ of size $r\timesr$, we find another valid decomposition $\Sigmay_t=\bm{\Lambda}^\prime\Sigmafac^\prime_t(\bm{\Lambda}^\prime)^\top+\Sigmaidi_t$, where $\bm{\Lambda}^\prime=\bm{\Lambda}\bm{P}^{-1}$ and $\Sigmafac^\prime_t=\bm{P}\Sigmafac_t\bm{P}^\top$. We resolve the ambiguity in the scale of the factors by fixing the level of their log-variance to zero, i.e., $\mufac_j=0$ for $j=1,\dots,r$. Sign and order identification can be enforced through restrictions on the factor loadings matrix $\bm{\Lambda}$. Several options are available in \pkg{factorstochvol} for restricting $\bm{\Lambda}$, for details see Section (ref).
Priors need to be specified for the mean, the latent log-variance processes, and for the factor loadings matrix $\bm{\Lambda}$. We choose $\beta_j\sim\Normal{b_\beta}{B_\beta}$, independently for $j=1,\dots,m$. For small values of $B_\beta$, this shrinks $\beta_j$ toward $b_\beta$; for large values of $B_\beta$, the prior is fairly uninformative.
The log-variance processes have the same prior specification as in the univariate case in Section (ref). For $\bm{\Lambda}$, three types of priors are currently implemented in \pkg{factorstochvol}. All three can be written in the form $\Lambda_{ij}\sim\Normal{0}{\tau^2_{ij}}$ independently for each applicable $i\in\{1,\dotsm\}$ and $j\in\{1,\dots,r\}$. First, one can fix all the $\tau^2_{ij}$s -- not necessarily to the same value -- a priori. This results in a normal prior for each element of the loadings matrix.
The second type is a hierarchical prior which has been developed to induce more flexible and potentially stronger shrinkage,
This distribution is termed normal gamma prior by griffin2010inference and implies a conditional variance $\mathbb{V}(\Lambda_{ij}\mid\lambda_i^2)$ of $2/\lambda_i^2$ and an unconditional excess kurtosis of $3/a$. The value of $a$ is treated as a structural parameter to be fixed by the user, where choosing $a$ small ($\lesssim 1$) enforces strong shrinkage towards zero, while choosing $a$ large ($\gtrsim 1$) imposes little shrinkage. The case $a = 1$ is a special case termed the Bayesian Lasso prior park2008bayesian. The parameter $\lambda_i^2$ is estimated from the data with $\lambda_i^2 \sim \mathcal{G}(c,d)$.
The third type is a slight modification of the second. Because variances in each row of the factor loadings matrix $\bm{\Lambda}$ can be seen as “random effects” from the same underlying distribution, the prior in Equation (ref) induces row-wise shrinkage with element-wise adaption. Analogously, one could also consider column-wise shrinkage with element-wise adaption, i.e.,
with the corresponding prior $\lambda_j^2 \sim \mathcal{G}(c,d)$.
Bayesian estimation in the factor SV model builds on the univariate vanilla SV implementations in \pkg{stochvol} and features several levels of efficiency boosting. To alleviate the problem of potentially slow convergence in high dimensions, it is carried out via a sampler that utilizes several variants of ASIS. The sampling details implemented in \pkg{factorstochvol} are described in kastner2017efficient as well as kastner2019sparse.
Similarly to \pkg{stochvol} and in an attempt to make computation time bearable even in higher dimensions, \pkg{factorstochvol}'s main sampler is written in \proglang{C++}. It uses the \proglang{R} package \pkg{Rcpp} to ease communication between \proglang{R} and \proglang{C++}. The univariate SV parts are borrowed from \pkg{stochvol} and interfaced through its \proglang{C}/\proglang{C++}-level updating function \code{update_fast_sv()}. In doing so, moving between interpreted \proglang{R} code and compiled \proglang{C++} code at each MCMC iteration is avoided.
The \pkg{stochvol} package provides means for fitting univariate SV, SVt, SVl, and SVtl models via its sampling routines \code{svsample()}, \code{svtsample()}, \code{svlsample()}, and \code{svtlsample()}, respectively. In the following, we describe a recommended workflow with \pkg{stochvol}. First, we discuss estimation, visualization, and prediction using default settings. Then, we show how to adapt the values of the prior hyperparameters and how to configure the sampling mechanism.
We estimate three models that exemplify the features and the user interface of \pkg{stochvol}. Using the \code{exrates} data found in the package, we model the EURCHF exchange rate (the price of 1 euro in Swiss franc) in the period between March 1, 2008 and March 1, 2012 (1028 data points) in three different ways.
The first example is an AR(1) model with SV residuals, i.e., Equation (ref) turns into
Using this model, we test whether the exchange rate follows a random walk with SV. In this case, we expect the posteriors of $\beta_0$ and $\beta_1$ to concentrate around 0 and 1, respectively.
In order to estimate this AR(1)-SV model, we need to prepare the input $\bm{y}$ as a \code{numeric} sequence of length $n$ and pass it as the first input argument to \code{svsample()} as follows:
We set \code{designmatrix = "ar1"} to use the AR(1) specification. More generally, \code{designmatrix} can take \code{character} values of the form \code{"ar0"} for a constant mean model, or \code{"ar1"}, \code{"ar2"}, etc., for AR(1), AR(2), and so on.
The second example is a constant mean model with SVt residuals, i.e., Equation (ref) becomes
If the returns are heavy-tailed, most of the posterior mass of $\nu$ concentrates on low values, e.g., smaller than 20. Otherwise, there is little evidence for high kurtosis.
We compute the log returns by applying \code{logret()} on the previously calculated \code{CHF_price}. Then, to estimate the constant mean model with heavy tailed SV residuals, we pass the vector of log returns to \code{svtsample()} with \code{designmatrix} set to \code{"ar0"}.
The third example is a multiple regression model with an intercept, two regressors, and SVl residuals; that is, Equation (ref) turns into
For illustration, we regress EURCHF log returns onto the contemporaneous log returns on EURUSD and EURJPY, the value of 1 euro per US dollar and Japanese yen, respectively.
To estimate a multiple regression model using \pkg{stochvol}, we need to prepare a \code{numeric} matrix $\bm X$ of dimension $n\timesK$, where rows correspond to time points and columns to covariates. We create an intercept as the first column of $\bm X$, and we set the second and the third columns to the EURUSD log returns and the EURJPY log returns, respectively; finally, we use the columns of $\bm X$ as covariates in the multiple regression.
Often, the joint posterior distribution of model parameters and latent quantities mark the goal of a Bayesian analysis. To inspect it, one can look at summary statistics and various types of visualizations of marginal posterior distributions. Also, it is recommended to examine the Markov chain for possible convergence issues -- this happens usually by investigating trace plots of posterior quantities. For this reason, inspired by the \pkg{coda} package, \pkg{stochvol} provides its own instances of the \proglang{R} generic functions \code{plot()} and \code{summary()}. In order to introduce the tools that \pkg{stochvol} provides for analyzing MCMC output, we briefly examine the results of the third example (multiple regression with SVl errors) in the remaining part of the section.
First, we plot the output of the estimation.
The result is shown in Figure (ref). We see in the first row a summary of the posterior density of the volatility. Apart from its median, we also receive a quantification of the uncertainty through the 5% and the 95% quantiles at each time point. In the second row, we can follow the evolution of the Markov chain of the SV parameters. In this example, they are $\mu$, $\varphi$, $\sigma$, and $\rho$. Lastly, we see prior and posterior density plots of the parameters in the third row in gray and black, respectively. They show high persistence and significant leverage.
Next, we observe the AR coefficients.
The result is shown in Figure (ref). On the left hand side, we do not spot any signs of convergence or mixing problems in the trace plots. On the right hand side, we see that none of the posterior densities of $\beta_0$, $\beta_1$, and $\beta_2$ concentrate around 0, hence the covariates seem to have an impact on the dependent variable.
As the final step, we print a numeric summary of the estimation results.
For brevity, we set \code{showlatent = FALSE} in order not to print all the 1027 latent states. The output shows the length of the burn-in and the number of draws, the prior specification of the parameters, and a concise summary of the marginal posterior distributions of the parameters $\mu$, $\varphi$, $\sigma$, and $\rho$, and additionally of the level of the volatility $\exp(\mu/2)$ and of $\sigma^2$, and of the vector of regression coefficients $\bm\beta$. This posterior summary is a table consisting of columns for the posterior mean and standard deviation, the 5%, 50%, and 95% quantiles. The user can influence the shown quantiles by passing a sequence of values between 0 and 1 to \code{svsample()}, \code{svtsample()}, \code{svlsample()}, or \code{svtlsample()} via the argument \code{quantiles}.
The last column in the table depicts the so-called effective sample size (ESS), a measure of the quality of a converged MCMC chain. Formally, ESS of a Markov chain $\mathrm{C}$ is defined through $M/(1+2\sum_{s=1}^\infty\rho^\text{eff}(s))$, where $M$ is the length of $\mathrm{C}$ and $\rho^\text{eff}(s)$ denotes the autocorrelation function for lag $s$ among the elements of $\mathrm{C}$. In principle, ESS is the sample size of a serially uncorrelated chain bearing the same Monte Carlo error as our (marginal) chain. Intuitively speaking, this means that ESS is the number of independent and identically distributed draws that were acquired and gives a sense of how well our chain has explored the posterior space. Higher values of ESS indicate better mixing.
We employ our estimated model to predict log returns for the remaining days in the data set. To do so, we first prepare the covariates for the next 24 days and pass them via the argument \code{newdata} of the generic \code{predict()} function along with the estimation output. Note that we need 25 days of price data to obtain 24 returns.
As we have access to the entire distribution of future log returns, we can quantify the uncertainty around our predictions through quantiles. In the following code snippet, we visualize the $k$-step-ahead predictive distributions for $k=1,\dots,24$, along with the truly observed values. The result is in Figure (ref).
Inspired by \code{ugarchroll()} in the \proglang{R} package \pkg{rugarch} rugarch, we introduce the suite of wrapper functions \code{svsample_roll()}, \code{svtsample_roll()}, \code{svlsample_roll()}, and \code{svtlsample_roll()}, built around their corresponding routines \code{svsample()}, \code{svtsample()}, \code{svlsample()}, and, respectively, \code{svtlsample()}, to simplify rolling window estimation of SV models. In this estimation method, either a fixed width time window is moving through the time series or a sequence of expanding time windows with the same starting time point covers larger and larger chunks of the observations, and the same model is estimated in all time windows independently. Next, each estimated model is employed for out-of-sample prediction, typically one day to one week ahead of the time window. Lastly, the set of predicted values might be used to evaluate the model fit.
In Bayesian statistics, a natural approach for assessing the predictive power of a model is through its posterior predictive distribution. Its density, also called the predictive density, is defined as
where $\bm\kappa$ collects all unobserved variables, i.e., $\bm\kappa=(\mu,\varphi,\sigma,\rho,\nu,\bm{h},\bm\beta)^\top$ in the most general case of SVtl, and the domain of integration $\bm{\mathcal{K}}$ is the set of all possible values for $\bm\kappa$. We follow geweke2010comparing in our notation by using a superscript $\bm o$ for the vector of observed values $\bm{y^o}_{[1:t]}=(y_1,y_2,\dots,y_t)^\top$. Equation (ref) can be seen as the integration of the predictive likelihood over the posterior distribution of all parameters and therefore it accounts for posterior parameter uncertainty for the predicted values.
The integral in Equation (ref) has no closed form and its dimensionality increases with $t$; it is intractable. Hence, we rely on Monte Carlo integration and we simulate from the posterior predictive distribution. For the evaluation of the predictive density at an observation $x=y^o_{t+1}$, called the predictive likelihood, we apply the computation
where $\bm{\kappa^{(m)}}$ denotes the $m$th posterior sample from the estimation procedure of the SV model. For other applications, the quantiles of the posterior predictive distribution, henceforth the predictive quantiles, might be of interest. We estimate the $q\%$ quantile through random variates simulated from $y_{t+1}\sim p(y_{t+1}\mid\bm{y^o}_{[1:t]})$, which we acquire by repeating two steps for $m=1,\dots,M$:
Lastly, we take the $q\%$ quantile of the sample vector $(y^{(1)}_{t+1}, y^{(2)}_{t+1}, \dots, y^{(M)}_{t+1})^\top$ as the approximate $q\%$ quantile of the predictive density. We implement the estimation of both the predictive likelihood and predictive quantiles in \pkg{stochvol}.
All four rolling window routines \code{svsample_roll()}, \code{svtsample_roll()}, \code{svlsample_roll()}, and \code{svtlsample_roll()} bear the same programming interface. They expect as their first argument the input data $\bm{y}_{[1:L]}$, which is of length $L$. For estimating the SV model in each time window $j=1,\dots,J$ in the moving or expanding window scheme, the sub-vector $\bm{y}_{[j:(t+j-1)]}$, or, respectively, $\bm{y}_{[1:(t+j-1)]}$, is taken as data and is used to predict $n_\text{ahead}\ge1$ time steps ahead. The width $t$ of the first time window can be determined from $L$, $J$, and $n_\text{ahead}$. The following example demonstrates how the rolling window sampling routines can be called in \pkg{stochvol}.
Argument \code{n_ahead} is used to set $n_\text{ahead}$, \code{forecast_length} is used to set $J$, and \code{refit_window} expects either \code{"moving"} or \code{"expanding"} to set the rolling window scheme to moving or expanding, respectively. Argument \code{calculate_quantile} expects a vector of numbers between 0 and 1; the numbers are interpreted as the quantiles to be predicted. Furthermore, if \code{calculate_predictive_likelihood} is set to \code{TRUE}, the function estimates the predictive likelihood. Lastly, the output \code{res} is a list of length $J$, i.e., one element for each time window. It contains the respective posterior quantile and predictive likelihood results together with all posterior parameter draws for $\bm\kappa$.
As discussed in Section (ref), the prior distributions need to be specified before the estimation process can start. Concerning the common model parameters $\mu$, $\varphi$, and $\sigma$, all of \code{svsample()}, \code{svtsample()}, \code{svlsample()}, and \code{svtlsample()} expect through their input arguments \code{priormu}, \code{priorphi}, and \code{priorsigma} values for $(b_\mu,\sqrt{B_\mu})$, $(a_\varphi,b_\varphi)$, and $B_\sigma$, respectively. Furthermore, all sampling functions accept the argument \code{priorbeta} to set an independent prior for the regression coefficients by providing $(b_\beta,s_\beta)$, where $b_\beta$ and $s_\beta$ are the common mean and, respectively, the common standard deviation. For a general multivariate normal distribution, the \code{specify_priors()} interface exists, which we detail later in this Section. The prior for $\nu$ can be influenced in \code{svtsample()} and \code{svtlsample()} by passing $\lambda_\nu$ as the argument \code{priornu}. Finally, \code{svlsample()} and \code{svtlsample()} take the \code{numeric} sequence $(a_\rho,b_\rho)$ through the input argument \code{priorrho}.
The code snippet below shows all the default values of the prior hyperparameters.
As an alternative to the concise interface above, a broader set of prior distributions can be specified via an object created by the \code{specify_priors()} function. The function has an input argument for each model parameter: \code{mu}, \code{phi}, \code{sigma2}, \code{nu}, \code{rho}, \code{beta}, and additionally one for the variance of $h_0$ called \code{latent0_variance}. There is a list of accompanying functions that create distributions in \pkg{stochvol}: \code{sv_beta()} has arguments \code{shape1} and \code{shape2} and it is accepted for \code{phi} and \code{rho}; \code{sv_constant()} has argument \code{value} and it is accepted for \code{mu}, \code{phi}, \code{sigma2}, \code{nu}, \code{rho}, and \code{latent0_variance}; \code{sv_normal()} has arguments \code{mean} and \code{sd} and it is accepted for \code{mu} and \code{phi}; \code{sv_multinormal()} has arguments \code{mean} and either \code{sd} and \code{dim} or \code{precision}, and it is accepted for \code{beta}; \code{sv_exponential()} has argument \code{rate} and it accepted for \code{nu}; \code{sv_gamma()} has arguments \code{shape} and \code{rate} and it is accepted for \code{sigma2}; \code{sv_inverse_gamma()} has arguments \code{shape} and \code{scale} and it is accepted for \code{sigma2}; and \code{sv_infinity()} has no arguments and it is accepted for \code{nu} hence turning the Student's $t$ distribution into a normal distribution. Additionally, \code{latent0_variance} accepts the character value \code{"stationary"}. All four sampling methods accept the prior specification object through the input argument \code{priorspec}.
All input arguments for \code{specify_priors} are optional, their default values and how they are used is seen below.
When conducting Bayesian inference using an MCMC sampling scheme, the number of draws from the posterior distribution, the length of the so-called burn-in phase, the initial values of the Markov chain, and the various strategies of storing the results are all of general interest. The input arguments \code{draws} and \code{burnin} settle the first two points. A sample size of \code{burnin + draws} is acquired from the posterior distribution out of which the first \code{burnin} number of draws are thrown away. The default is to draw 10000 elements after a burnin of 1000 for SV models without leverage, and draw 20000 elements after a burnin of 2000 for SV models with leverage, which in our experience is enough for most applications.
As for the initial values, \code{startpara} and \code{startlatent} provide a way to set them. The argument \code{startpara} is expected to be a named \code{list} mapping parameter names to starting values, and \code{startlatent} must be a sequence of length $m$ that contains starting values for $\bm{h}$. Default values are set to be the prior mean for $\varphi$, $\sigma$, $\nu$, and $\rho$, these have only minor influence on the Markov chain. The default value for $\bm\beta$ is the ordinary least squares estimator $\bm{(X^\top X)^{-1}\bm{X^\top}}\bm{y}$, where $\bm X$ denotes the regression design matrix and $\bm{y}$ denotes the vector of observations. After setting $\bm\beta$, the level of log-variance $\mu$ is initialized according to the Bayesian linear regression
where $\xi_t\sim\Normal{-1.27}{4.934}$. Equation (ref) results from the first line of Equation (ref) by fixing $h_t$ at its stationary expected value $\mu$ and then taking $x\mapsto\log(x^2)$ of both sides. The homoskedastic error term $\xi_t$ is acquired as the Laplace approximation to $\log(\varepsilon_t^2)$ harvey1996estimation. At the end, by default all values of the vector \code{startlatent} are set to the initial value of $\mu$.
It is customary to start independent Markov chains in parallel and \pkg{stochvol} provides facilities for that in all of its sampling procedures. The argument \code{n_chains} is expected to be a positive integer, it sets the number of independent chains. Additionally, arguments \code{parallel}, \code{n_cpus}, and \code{cl} can be used to control parallelism used by \pkg{stochvol}. To overwrite the default sequential execution strategy, \code{parallel} is to be set either to \code{"snow"}, to employ the so-called “SNOW” clusters, or to \code{"multicore"} to use the “multicore” type computation rlanguage. Next, argument \code{n_cpus} should be set to the physical number of parallel processing units to be used. Finally, in case “SNOW” is applied, the sampling routines optionally accept an already running “SNOW” cluster through argument \code{cl}.
As mentioned earlier, the sampling algorithms for the latent states $\bm{h}$ in \pkg{stochvol} rely on a Gaussian mixture approximation as in omori2007stochastic and kastner2014ancillarity. The approximation tends to be very good, therefore the default setting is not to correct for model misspecification. However, this correction can be enabled in all of the sampling routines through the \code{expert} argument as shown for \code{svsample()} in the following.
Lastly, \pkg{stochvol} provides three ways to economize storage during and after the execution of the sampler. Setting the \code{integer} argument \code{thinpara} to $\iota$ tells the sampler to store only every $\iota$th draw of the vector of parameters, and supplying a value for \code{thinlatent} does the same for $\bm{h}$. Finally, one has the opportunity not to store the full vector $\bm{h}$ but only its last value by setting \code{keeptime = "last"}. The default behavior is to store every draw after the burn-in phase.
The most common workflow of using \pkg{factorstochvol} for fitting multivariate factor SV models consists of the following steps: (1) Prepare the data, (2) decide on an identification structure, (3) specify the prior hyperparameters, (4) run the sampler, (5) investigate the output and visualize the results, and (6) predict (if required). These steps are described in detail in the following sections.
The workhorse in \pkg{factorstochvol} is the sampling function \code{fsvsample()}. It expects the data to come in form of a matrix $\bm Y=(\bm{y}_1,\dots,\bm{y}_n)^\top$ with $n$ rows and $m$ columns. For illustration, we use the exchange rate data set in \pkg{stochvol} which contains $3140$ daily observations of exchanges rates for $23$ currency pairs against EUR, ranging from March 3, 2000 to April 4, 2012. To keep the analysis simple and computation times moderate, we however only model the last $1001$ days of the first six series in alphabetical order (Australian dollar, Canadian dollar, Swiss franc, Czech koruna, Danish krone, Great British pound) for further analysis. Instead of using the nominal exchange rates we compute log returns. This leaves us with a data set of size $n = 1000$ and $m = 6$. The data is prepared using the code snippet below and visualized in Figure (ref) using the \pkg{zoo} package zoo.
The likelihood in factor models is invariant to certain factor transformations such as reordering of factors and their loadings or sign switches thereof. In addition to this, it is often multimodal. Consequently, identifying the factor loadings is far from trivial. The most common way to address this issue in factor SV models is to impose a lower-diagonal factor loadings matrix where all elements above the diagonal are set to zero aguilar2000bayesian, chib2006analysis, han2006asset, zhou2014bayesian. To use this constraint in \pkg{factorstochvol}, the argument \code{restrict = "upper"} can be passed to the main sampling function \code{fsvsample()}. Evidently, this practice imposes an order dependence, as, e.g., the first variable is not allowed to load on anything else but the first factor.
A rather ad hoc method for automatically ordering the data is implemented in the helper function \code{preorder()}. After a maximum likelihood factor model fit to the data (using \code{factanal()} from the \pkg{stats} package with the default \code{varimax} rotation), the series are ordered as follows: The variable with the highest loading on factor 1 is placed first, the variable with the highest loading on factor 2 second (unless this variable is already placed first, in which case the variable with the second highest loading is taken), et cetera. For the data set at hand, this would imply the following ordering for a two-factor model.
According to this algorithm, the second series should be placed first and the third series should be placed second. Thereafter, the alphabetical ordering remains.
To achieve this effect without reordering the data, a \code{logical} matrix of size $m \times r$ can be passed to \code{fsvsample()} via \code{restrict}, where the entry \code{TRUE} means that this element is restricted to zero; \code{FALSE} means that it is to be estimated from the data. Similarly to \code{preorder()}, the function \code{findrestrict()} tries to automate this procedure. Again, the maximum likelihood estimates from a static factor analysis are used; however, \code{findrestrict()} uses a slightly different algorithm than the one above: The variable with the lowest absolute loadings on factors $2, 3, \dots, r$ (relative to factor 1) is determined to lead the first factor, the variable with the lowest absolute loadings on factors $3, 4, \dots, r$ (relative to factors 1 and 2) is placed second, et cetera. Below is the result for the data set at hand.
If \code{fsvsample()} is called with the argument \code{restrict = "auto"}, it automatically invokes \code{findrestrict()} with the appropriate number of factors. Using \code{restrict = "none"} (the default) causes the sampler not to place any constraints on the loadings matrix; thus, the resulting posterior draws may be unstable or suffer from multiple local modes. If, however, inference on the factor loadings themselves is not the primary concern of the analysis, leaving the factor loadings unidentified may be the preferred option. This is in particular the case when inference for the covariance matrix is sought, as this only depends on $\bm{\Lambda}$ through the rotation-invariant transformation of Equation (ref). For a more elaborate discussion of these issues, we refer the reader to sentana2001identification who discuss automatic identification through hetero\-skedasticity. A comparison of log predictive scores under different identification schemes for factor SV models is given in kastner2017efficient; see also fruhwirth2018sparse for related issues in static factor models. To continue with the current example, we chose not to place any a priori restrictions on the factor loadings matrix while using a row-wise normal-gamma shrinkage prior on the factor loadings matrix kastner2019sparse.
Apart from the obvious prior choice about the number of factors and the identification scheme discussed above, a number of hyperparameter choices are available in \pkg{factorstochvol}. Regarding the log-variance processes, the interface is analogous to that of \code{svsample()}. In the following, $i=1,\dots,m$ and $j=1,\dots,r$ index the idiosyncratic and the factor log-variance processes, respectively. The pair of common prior hyperparameters $(b_\beta, B_\beta)$ can be passed as a sequence of length two to \code{priorbeta}. The common prior of $\muidi_i$ can be set by passing a sequence of length two -- the mean and the standard deviation of the normal distribution -- to \code{priormu}; the common priors of $\phiidi_i$ and $\phifac_j$ can also be set by passing sequences of length two -- the parameters of the corresponding beta distribution -- to \code{priorphiidi} and to \code{priorphifac}, respectively; similarly, the common priors of $\sigmaidi_i$ and $\sigmafac_j$ can be specified via the arguments \code{priorsigmaidi} and \code{priorsigmafac}, respectively, that accept as positive numbers the scale $B_\sigma$ of the corresponding gamma distribution.
As discussed in Section (ref), \pkg{factorstochvol} offers three specifications as priors for $\bm{\Lambda}$, controlled through the argument \code{priorfacloadtype}. To use the first option (\code{priorfacloadtype = "normal"}), one needs to fix the values of $\tau_{ij}$ a priori. The user can pass these fixed values to \code{fsvsample()} via the argument \code{priorfacload}, either as an $m\timesr$ matrix with positive entries or as a single positive number which will be recycled accordingly. For the second option, the normal gamma prior with row-wise or column-wise shrinkage (\code{priorfacloadtype = "rowwiseng"} and \code{priorfacloadtype = "colwiseng"}, respectively), the value of argument \code{priorfacload} is then interpreted as the shrinkage parameter $a$. Both specifications of the normal gamma prior need the values $c$ and $d$. They can be set as a two-element vector passed to the argument \code{priorng}.
Running the sampler corresponds to invoking \code{fsvsample()}. Apart from the prior settings discussed above, its most important arguments are listed below with the default value in brackets. For a complete list of all arguments and more details, see \code{?fsvsample}.
For our illustrative example, most settings are left at their default values. The number of factors is increased from one to two, instead of 1000 we sample 10000 draws, we estimate a constant mean, a thinning of 10 is used, and \code{quiet} is set to \code{TRUE}.
The resulting object
holds a rich amount of information. In particular, it contains
For more details, please investigate \code{str(res)} and/or \code{?fsvsample}.
Using \code{covmat()}, one can extract the MCMC draws of the implied covariance matrices for all points in time which have been stored during sampling. By default, this is the last point in time (\code{keeptime = "last"}), and thus
shows that we have stored $1000$ posterior draws of a $6\times6$ covariance matrix at one point in time, $t = n = 1000$. To check convergence, one can take a look at the trace plot and the autocorrelation function of the log determinant, i.e.,
The results are visualized in Figure (ref); decent mixing for this quantity is apparent.
To assess the mixing speed of each individual covariance matrix element, one can check, e.g., the estimated effective sample size (out of $1000$ draws kept) which is implemented in \pkg{coda}. Again, no major convergence problems are apparent.
Assuming that \code{runningstore} was set sufficiently high when sampling, several convenience functions can be used for quick visualizations without having to post-process the MCMC draws. For example, to visualize the time-varying correlation matrices, consider
which produces the three estimated posterior correlation matrices depicted in Figure (ref). Setting \code{plotCI = "circle"} visualizes posterior uncertainty -- inner and outer radii correspond to the posterior mean plus/minus two standard deviations, respectively.
To get an idea about how the marginal volatilities evolve over time, \code{voltimeplot()} can be used. To exemplify,
yields the estimated volatilities in Figure (ref). The financial crisis of 2008 and the capping of CHF's appreciation in September 2011 are clearly visible, while DKK's volatility (relative to EUR) is practically zero. Note that \code{voltimeplot()} respects palette changes. In the above example, \pkg{RcolorBrewer} rcolorbrewer is used. Moreover,
yields the estimated pairwise correlations in Figure (ref). While, relative to EUR, the estimated correlation between AUD and CAD appears to be relatively stable over time, correlations with CHF can become negative at times. To visualize the communalities, i.e., the proportions of variances explained through the latent factors, invoke
which yields the estimated communalities in Figure (ref).
To gain an even deeper understanding of the estimated model, we now turn towards examining the latent factors and their variances themselves. To visualize the loadings, the functions \code{facloadpairplot()}, \code{facloadcredplot()}, \code{facloadpointplot()}, \code{facloadtraceplot()}, and \code{facloaddensplot()} are available; the former two are exemplified in Figure (ref). Moreover, we can see the factor log variances produced through \code{logvartimeplot(res, show = "fac")}. Similarly, \code{logvartimeplot(res, show = "idi")} produces plots of the idiosyncratic log variances which are displayed in Figure (ref).
In addition to the above, there is the plotting function \code{paratraceplot()} which produces trace plots of all parameters associated with the log variances processes: mean, persistence, and volatility of log variances.
In order to provide some guidance when it comes to selecting the number of factors, \pkg{factorstochvol} ships with the function \code{evdiag()}. It computes and visualizes the eigenvalues of $\bm{\Lambda}^\top\bm{\Lambda}$ which can be used as a rough guide in analogy to a scree plot for static factor models. To use it, one can fit a model with a relatively large number of factors and assess the importance of each of these (in descending order). For example, the code below produces Figure (ref), hinting at a model with no more than four factors.
One of the main use cases of \pkg{factorstochvol} might be to predict covariance and correlation matrices of time series. To this end, \code{predcov()} and \code{predcor()} yield draws from the posterior predictive distribution of these, defined in analogy to Equation (ref). For instance, the code below can be used to obtain one-step-ahead posterior predictive means and standard deviations for the correlation matrix on April 5, 2012 (using data up to April 4 only).
To obtain draws from the posterior predictive distribution of new data points, one can simply draw from the corresponding mixture of multivariate normals. In Figure (ref), these draws are visualized via \code{heatpairs()} from \pkg{LSD} LSD.
To conclude, we note that convenience functions such as \code{predloglik()} or \code{predloglikWB()} may aid in approximating predictive likelihoods geweke2010comparing. For instance, assuming that the actually observed value of $\bm y_{n+1} = \bm y_{n+2} = (0,0,0,0,0,0)^\top$, we can approximate the one and two step ahead log predictive scores through code along the following lines.
For a more detailed analysis of the out-of-sample performance of \pkg{factorstochvol} for exchange rate data, we refer to kastner2017efficient; for a predictive exercise on stock market data, please see kastner2019sparse.
We extended the work of kastner2016dealing to other SV models, including the univariate heavy-tailed SV, the SV model with leverage, and the multivariate factor SV model. We showcased the features that are the most important to end users in \proglang{R}: estimation through the sampler functions, visualization, summary, and prediction methods. Due to its more involved nature, however, we did not include the description of the \proglang{C++} interface. Two functions called \code{update_fast_sv()} and \code{update_general_sv()} are exported and programmers have the possibility to access the samplers in \pkg{stochvol} directly from \proglang{C++} after linking to the compiled package. For usage examples, see the implementations of \pkg{factorstochvol} or \pkg{shrinkTVP} shrinktvp. Source code is available in the GitHub repositories \url{https://github.com/gregorkastner/stochvol} and \url{https://github.com/gregorkastner/factorstochvol}.
Gregor Kastner acknowledges funding from the Austrian Science Fund (FWF) for the project “High-dimensional statistical learning: New methods to advance economic and sustainability policy” (ZK 35), jointly carried out by the University of Klagenfurt, WU Vienna University of Economics and Business, Paris Lodron University Salzburg, TU Wien, and the Austrian Institute of Economic Research (WIFO).