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.
60,981 characters · 21 sections · 36 citation commands
Dealing with Stochastic Volatility in Time Series Using the R Package stochvol
\shortcites{r:mvt}
This vignette corresponds to an article of the same name which is published in the Journal of Statistical Software. At the time of this writing, the vignette and the published article essentially coincide (in some trace plots within this vignette, only every 10th value is shown to keep \pkg{stochvol} below the 5MB mark). The version at hand might receive minor updates as time goes by. To cite, please use kas:dea. Further information about citing \pkg{stochvol} can be obtained in \proglang{R} by installing the package, e.g., through \code{install.packages("stochvol")}, and calling \code{citation("stochvol")}.
Returns -- in particular financial returns -- are commonly analyzed by estimating and predicting potentially time-varying volatilities. This focus has a long history, dating back at least to mar:por who investigates portfolio construction with optimal expected return-variance trade-off. In his article, he proposes rolling-window-type estimates for the instantaneous volatilities, but already then recognizes the potential for “better methods, which take into account more information”.
One approach is to model the evolution of volatility deterministically, i.e., through the (G)ARCH class of models. After the groundbreaking papers of eng:aut and bol:gen, these models have been generalized in numerous ways and applied to a vast amount of real-world problems. As an alternative, tay:fin proposes in his seminal work to model the volatility probabilistically, i.e., through a state-space model where the logarithm of the squared volatilities -- the latent states -- follow an autoregressive process of order one. Over time, this specification became known as the stochastic volatility (SV) model. Even though several papers jac-etal:bayJBES, ghy-etal:sto, kim-etal:sto provide early evidence in favor of using SV, these models have found comparably little use in applied work. This obvious discrepancy is discussed in bos:rel who points out two reasons: the variety (and potential incompatibility) of estimation methods for SV models -- whereas the many variants of the GARCH model have basically a single estimation method -- and the lack of standard software packages implementing these methods.
In kas-fru:anc, the former issue is thoroughly investigated and an efficient MCMC estimation scheme is proposed. The paper at hand and the corresponding package \pkg{stochvol} r:sto for \proglang{R} r:r are crafted to cope with the latter problem: the apparent lack of ready-to-use software packages for efficiently estimating SV models.
We begin by briefly introducing the model and specifying the notation used in the remainder of the paper. Furthermore, an overview of Bayesian parameter estimation via Markov chain Monte Carlo (MCMC) methods is given.
Let $\bm{y} = (y_1, y_2, \dots, y_n)^\top$ be a vector of returns with mean zero. The intrinsic feature of the SV model is that each observation $y_t$ is assumed to have its “own” contemporaneous variance $e^{h_t}$, thus relaxing the usual assumption of homoskedasticity. In order to make the estimation of such a model feasible, this variance is not allowed to vary unrestrictedly with time. Rather, its logarithm is assumed to follow an autoregressive process of order one. Note that this feature is fundamentally different to GARCH-type models where the time-varying volatility is assumed to follow a deterministic instead of a stochastic evolution.
The SV model can thus be conveniently expressed in hierarchical form. In its centered parameterization, it is given through
where $\mathcal{N}\!\left(\mu, {\sigma_\eta^2}\right)$ denotes the normal distribution with mean $\mu$ and variance ${\sigma_\eta^2}$. We refer to $\bm{\theta} = (\mu, \phi, {\sigma_\eta})^\top$ as the vector of parameters: the level of log-variance $\mu$, the persistence of log-variance $\phi$, and the volatility of log-variance ${\sigma_\eta}$. The process $\bm{h}=(h_{0}, h_{1}, \ldots,h_{n})$ appearing in Equation (ref) and Equation (ref) is unobserved and usually interpreted as the latent time-varying volatility process (more precisely, the log-variance process). Note that the initial state $h_0$ appearing in Equation (ref) is distributed according to the stationary distribution of the autoregressive process of order one.
To complete the model setup, a prior distribution for the parameter vector $\bm{\theta}$ needs to be specified. Following kim-etal:sto, we choose independent components for each parameter, i.e., $p(\bm{\theta}) = p(\mu)p(\phi)p({\sigma_\eta})$.
The level $\mu \in \mathbb{R}$ is equipped with the usual normal prior $\mu \sim \mathcal{N}\!\left(b_{\mu}, B_{\mu}\right)$. In practical applications, this prior is usually chosen to be rather uninformative, e.g., through setting $b_\mu = 0$ and $B_\mu \geq 100$ for daily log returns. Our experience with empirical data is that the exact choice is usually not very influential; see also Section (ref).
For the persistence parameter $\phi \in (-1,1)$, we choose $(\phi+1)/2 \sim \mathcal{B}\left(a_0, b_0\right)$, implying
where $a_0$ and $b_0$ are positive hyperparameters and $B (x,y) = \int_0^1t^{x-1}(1-t)^{y-1}\,dt$ denotes the beta function. Clearly, the support of this distribution is the interval $(-1, 1)$; thus, stationarity of the autoregressive volatility process is guaranteed. Its expected value and variance are given through the expressions
This obviously implies that the prior expectation of $\phi$ depends only on the ratio $a_0:b_0$. It is greater than zero if and only if $a_0 > b_0$ and smaller than zero if and only if $a_0 < b_0$. For a fixed ratio $a_0:b_0$, the prior variance decreases with larger values of $a_0$ and $b_0$. The uniform distribution on $(-1,1)$ arises as a special case when $a_0=b_0=1$. For financial datasets with not too many observations (i.e., $n \lesssim 1000$), the choice of the hyperparameters $a_0$ and $b_0$ can be quite influential on the shape of the posterior distribution of $\phi$. In fact, when the underlying data-generating process is (near-)homoskedastic, the volatility of log-variance ${\sigma_\eta}$ is (very close to) zero and thus the likelihood contains little to no information about $\phi$. Consequently, the posterior distribution of $\phi$ is (almost) equal to its prior, no matter how many data points are observed. For some discussion about this issue, see, e.g., kim-etal:sto who choose $a_0=20$ and $b_0=1.5$, implying a prior mean of $0.86$ with a prior standard deviation of $0.11$ and thus very little mass for nonpositive values of $\phi$.
For the volatility of log-variance ${\sigma_\eta} \in \mathbb{R}^+$, we choose ${\sigma_\eta^2} \sim B_{{\sigma_\eta}}\times \chi^2_1= \mathcal{G}\left(1/2,1/2B_{{\sigma_\eta}}\right)$. This choice is motivated by fru-wag:sto who equivalently stipulate the prior for $\pm \sqrt { {\sigma_\eta^2}}$ to follow a centered normal distribution, i.e., $\pm \sqrt { {\sigma_\eta^2}} \sim \mathcal{N}\!\left(0,B_{{\sigma_\eta}}\right)$. As opposed to the more common Inverse-Gamma prior for ${\sigma_\eta^2}$, this prior is not conjugate in the usual sampling scheme. However, it does not bound ${\sigma_\eta^2}$ away from zero a priori. The choice of the hyperparameter $B_{{\sigma_\eta}}$ turns out to be of minor influence in empirical applications as long as it is not set too small.
An MCMC algorithm such as the one implemented in the package \pkg{stochvol} provides its user with draws from the posterior distribution of the desired random variables; in our case with the latent log-variances $\bm{h}$ and the parameter vector $\bm{\theta}$. Because these draws are usually dependent, Bayesian inference via MCMC may require careful design of the algorithm and attentive investigation of the draws obtained.
One key feature of the algorithm used in this package is the joint sampling of all instantaneous volatilities “all without a loop” (AWOL), a technique going back at least to rue:fas and discussed in more detail in mcc-etal:sim. Doing so reduces correlation of the draws significantly and requires auxiliary finite mixture approximation of the errors as in kim-etal:sto or omo-etal:sto.
In order to avoid the cost of code interpretation within each MCMC iteration, the core computations are implemented in \proglang{C}. Their output is interfaced to \proglang{R} via the \pkg{Rcpp} package edd-fra:rcp; there, the convenience functions and the user interface are implemented. This combination allows to make use of the well-established and widely accepted ease-of-use of \proglang{R} and its underlying functional programming paradigm. Moreover, existing frameworks for analyzing MCMC output such as \pkg{coda} plu-etal:cod as well as high-level visualization tools can easily be used. Last but not least, users with a basic knowledge of \proglang{R} can use the package with a very low entry cost. Nevertheless, despite all these convenience features, the package profits from a highly optimized machine code generated by a compiler at package build time, thus providing acceptable runtime even for larger datasets.
A novel and crucial feature of the algorithm implemented in \pkg{stochvol} is the usage of a variant of the “ancillarity-sufficiency interweaving strategy” (ASIS) which has been brought forward in the general context of state-space models by yu-men:cen. ASIS exploits the fact that for certain parameter constellations sampling efficiency improves substantially when considering a non-centered version of a state-space model. This is commonly referred to as a reparameterization issue and dealt with extensively in scholarly literature; for an early reference see, e.g., hil-smi:par. For the model at hand, a move of this kind can be achieved by transferring the level of log-variance $\mu$ and/or the volatility of log-variance ${\sigma_\eta}$ from the state process (Equations (ref) and (ref)) to the observation process (Equation (ref)) through a simple reparameterization of $\bm{h}$. However, in the case of the SV model, it turns out that no single superior parameterization exists. Rather, for some underlying processes, the standard parameterization yields superior results, while for other processes non-centered versions are better. To overcome this issue, the parameter vector $\bm{\theta}$ is sampled twice: once in the centered and once in a noncentered parameterization. This method of “combining best of different worlds” allows for efficient inference regardless of the underlying process with one algorithm. For more details about the algorithm and empirical results concerning sampling efficiency, see kas-fru:anc.
The usual stand-alone approach to fitting SV models with \pkg{stochvol} exhibits the following workflow: (1) Prepare the data, (2) specify the prior distributions and configuration parameters, (3) run the sampler, (4) assess the output and display the results. All these steps are described in more detail below, along with a worked example. For a stepwise incorporation of SV effects into other MCMC samplers, see Section (ref).
The core sampling function \code{svsample} expects its input data \code{y} to be a numeric vector of returns without any missing values (\code{NA}s) and throws an error if provided with anything else. In the case that \code{y} contains zeros, a warning is issued and a small offset constant of size \code{sd(y)/10000} is added to the squared returns before doing the auxiliary mixture sampling omo-etal:sto. However, we generally recommend to avoid zero returns altogether, e.g., by demeaning them beforehand.
Below is an illustration of how to prepare data by using the sample dataset \code{exrates}\footnote{The dataset -- obtained from the European Central Bank's Statistical Data Warehouse -- contains the daily bilateral prices of one euro in 23 currencies from January 3, 2000 until April 4, 2012. Conversions to New Turkish lira and Fourth Romanian leu have been incorporated. See \code{?exrates} for more information.} included in the package. Figure (ref) provides a visualization of a time series from this dataset.
Additionally to real-world data, \pkg{stochvol} also has a built-in data generator \code{svsim}. This function simply produces realizations of an SV process and returns an object of class \code{svsim} which has its own \code{print}, \code{summary}, and \code{plot} methods. Exemplary code using \code{svsim} is given below and the particular instance of this simulated series is displayed in Figure (ref).
After preparing the data vector \code{y}, the user needs to specify the prior hyperparameters for the parameter vector $\bm{\theta} = (\mu, \phi, {\sigma_\eta})^\top$ -- see also Section (ref) -- and some configuration parameters. The appropriate values are passed to the main sampling function \code{svsample} as arguments which are described below.
The argument \code{priormu} is a vector of length $2$, containing mean and standard deviation of the normal prior for the level of the log-variance $\mu$. A common strategy is to choose a vague prior here, e.g., \code{c(0, 100)}, because the likelihood usually carries enough information about this parameter. If one prefers to use (slightly) informative priors, e.g., to avoid outlier draws of $\mu$, one must pay attention to whether log returns or percentage log returns are analyzed. Assuming daily data, log returns commonly have an unconditional variance of $0.0001$ or less and thus the level on the log scale $\mu$ lies around $\log(0.0001) \approx -9$. Percentage log returns, on the other hand, have the $100^2$-fold unconditional variance (around $1$) which implies a level of $\log(1) = 0$. Choices in the literature include \code{c(0, 10)} jac-etal:bayJE, \code{c(0, 5)} yu:on, \code{c(0, sqrt(10))} kim-etal:sto, mey-yu:bug or \code{c(0, 1)} omo-etal:sto. Note that most of these choices are quite informative and clearly designed for percentage log returns.
For specifying the prior hyperparameters for the persistence of log-variance, $\phi$, the argument \code{priorphi} may be used. It is again a vector of length $2$, containing $a_0$ and $b_0$ specified in Equation (ref). As elaborated in Section (ref), these values can possibly be quite influential, thus we advise to choose them carefully and study the effects of different choices. The default is currently given through \code{c(5, 1.5)}, implying a prior mean of $0.54$ and a prior standard deviation of $0.31$.
The prior variance of log-variance hyperparameter $B_{\sigma_\eta}$ may be controlled through \code{priorsigma}. This argument defaults to \code{1} if not provided by the user. As discussed in Section (ref), the exact choice of this value is usually not very influential in typical applications. In general, it should not be set too small unless there is a very good reason, e.g., explicit prior knowledge, to do so.
For specifying the size of the burn-in, the parameter \code{burnin} is provided. It is the amount of MCMC iterations that are run but discarded to ensure convergence to the stationary distribution of the chain. The current default value for this parameter is $1000$ which has turned out to suffice in most situations. Nevertheless, the user is encouraged to check convergence carefully; see Section (ref) for more details. The amount of iterations which are run after burn-in can be specified through the parameter \code{draws}, currently defaulting to $10\,000$. Consequently, the sampler is run for a total of \code{burnin + draws} iterations.
Three thinning parameters are available which all are $1$ if not specified otherwise. The first one, \code{thinpara}, is the denominator in the fraction of parameter draws (i.e., draws of $\bm{\theta}$) that are stored. E.g., if \code{thinpara} equals $10$, every $10$th draw is kept. The default parameter thinning value of $1$ means that all draws are saved. The second thinning parameter, \code{thinlatent}, acts in the same way for the latent variables $\bm{h}$. The third thinning parameter, \code{thintime}, refers to thinning with respect to the time dimension of the latent volatility. In the case that \code{thintime} is greater than 1, not all elements of $\bm{h}$ are stored, e.g., for \code{thintime} equaling $10$, only the draws of $h_1, h_{11}, h_{21}, \dots$ (and $h_0$) are kept.
Another configuration argument is \code{quiet} which defaults to \code{FALSE}. If set to \code{TRUE}, all output during sampling (progress bar, status messages) is omitted. The arguments \code{startpara} and \code{startlatent} are optional starting values for the parameter vector $\bm{\theta}$ and the latent variables $\bm{h}$, respectively. All other configuration parameters are summarized in the argument \code{expert}, because it is not very likely that the end-user needs to mess with the defaults.\footnote{Examples of configurations that can be changed with the \code{expert}-argument include the specification of the (baseline) parameterization (either centered or noncentered) and the possibility to turn off interweaving. Moreover, some algorithmic details such as the number of blocks used for the parameter updates or the possibility of using a random walk Metropolis-Hastings proposal (instead of the default independence proposal) can be found here.} Please refer to the package documentation and kas-fru:anc for details.
Any further arguments (\code{\dots}) are forwarded to \code{updatesummary}, controlling the type of summary statistics that are calculated for the posterior draws.
At the heart of the package \pkg{stochvol} lies the function \code{svsample} which serves as an \proglang{R}-wrapper for the actual sampler coded in \proglang{C}. Exemplary usage of this function is given in the code snipped below, along with the default output.
As can be seen, this function calls the main MCMC sampler and converts its output to \pkg{coda}-compatible objects. The latter is done mainly for reasons of compatibility and in order to have straightforward access to the convergence diagnostics checks implemented there. Moreover, some summary statistics for the posterior draws are calculated. The return value of \code{svsample} is an object of type \code{svdraws} which is a named list with eight elements, holding (1) the parameter draws in \code{para}, (2) the latent log-volatilities in \code{latent}, (3) the initial latent log-volatility draw in \code{latent0}, (4) the data provided in \code{y}, (5) the sampling runtime in \code{runtime}, (6) the prior hyperparameters in \code{priors}, (7) the thinning values in \code{thinning}, and (8) summary statistics of these draws, alongside some common transformations thereof, in \code{summary}.
Following common practice, \code{print} and \code{summary} methods are available for \code{svdraws} objects. Each of these has two optional parameters, \code{showpara} and \code{showlatent}, specifying which output should be displayed. If \code{showpara} is \code{TRUE} (the default), values/summaries of the parameter draws are shown. If \code{showlatent} is \code{TRUE} (the default), values/summaries of the latent variable draws are shown. In the example below, the summary for the parameter draws only is displayed.
There are several plotting functions specifically designed for objects of class \code{svsample} which are described in the following paragraphs.
The generic \code{plot} method for \code{svdraws} objects combines all above plots into one plot. All arguments described above can be used. See \code{?plot.svsample} for an exhaustive summary of possible arguments and Figure (ref) for an example.
For extracting standardized residuals, the \code{residuals}/\code{resid} method can be used on a given \code{svdraws} object. With the optional argument \code{type}, the type of summary statistic may be specified. Currently, \code{type} is allowed to be either \code{"mean"} or \code{"median"}, where the former corresponds to the default value. This method returns a real vector of class \code{svresid} which contains the requested summary statistic of standardized residuals for each point in time. There is also a \code{plot} method available, providing the option of comparing the standardized residuals to the original data when given through the argument \code{origdata}. See the code below for an example and Figure (ref) for the corresponding output.
We demonstrate how the \pkg{stochvol} package can be used to incorporate stochastic volatility into any given MCMC sampler. This is particularly easy when the sampler itself is coded in \proglang{R} or \proglang{C}/\proglang{C++} and applies, e.g., to many of the specialized procedures listed in the CRAN task view about Bayesian inference, available at \url{http://cran.r-project.org/web/views/Bayesian.html}. For the sake of simplicity, we explain the procedure using a “hand-coded” Gibbs-sampler for the Bayesian normal linear model with $n$ observations and $k=p-1$ predictors, given through
Here, $\bm{y}$ denotes the $n\times1$ vector of responses, $\bm{X}$ is the $n \times p$ design matrix containing ones in the first column and the predictors in the others, and $\bm{\beta} = (\beta_0, \beta_1, \dots, \beta_{p-1})^\top$ stands for the $p\times 1$ vector of regression coefficients. In the following sections, we discuss two specifications of the $n\times n$ error covariance matrix $\bm{\Sigma}$.
The arguably simplest specification of the error covariance matrix in Equation (ref) is given by $\bm{\Sigma}\equiv\sigmaepsextrasuper2\bm{I}$, where $\bm{I}$ denotes the $n$-dimensional unit matrix. This specification is used in many applications and commonly referred to as the linear regression model with homoskedastic errors. To keep things simple, let model parameters $\bm{\beta}$ and $\sigmaepsextrasuper2$ be equipped with the usual conjugate prior $p(\bm{\beta}, \sigmaepsextrasuper2) = p(\bm{\beta}|\sigmaepsextrasuper2)p(\sigmaepsextrasuper2)$, where
Commonly, samples from the posterior distribution of this model are obtained through a Gibbs-algorithm, where draws are generated in turn from the full conditional distributions $\bm{\beta}|\bm{y},\sigmaepsextrasuper2 \sim \mathcal{N}\!\left(\bm{b}_T, \bm{B}_T\right)$ with \[ \bm{b}_T = \left(\bm{X}^\top \bm{X} + \Bb0^{-1} \right)^{-1}\left( \bm{X}^\top\bm{y} + \Bb0^{-1}\bb0 \right), \qquad \bm{B}_T = \sigmaepsextrasuper2\left(\bm{X}^\top \bm{X} + \Bb0^{-1} \right)^{-1}, \]
and $\sigmaepsextrasuper2|\bm{y},\bm{\beta} \sim \mathcal{G}^{-1}\!\left(c_n, C_n\right)$ with \[ c_n = c_0 + \frac{n}{2} + \frac{p}{2}, \qquad C_n = C_0 + \frac{1}{2} \left( (\bm{y} - \bm{X}\bm{\beta})^\top(\bm{y} - \bm{X}\bm{\beta}) + (\bm{\beta} - \bb0)^\top\Bb0^{-1}(\bm{\beta} - \bb0) \right). \]
In \proglang{R}, this can straightforwardly be coded as follows:
Instead of homoskedastic errors, we now specify the error covariance matrix in Equation (ref) to be $\bm{\Sigma} \equiv \text{diag}({e^{h_1}, \dots, e^{h_n}})$, thus introducing nonlinear dependence between the observations due to the AR(1)-nature of $\bm{h}$. Instead of cooking up an entire new sampler, we adapt the code from above utilizing the \pkg{stochvol} package. To do so, we simply replace the sampling step of $\sigma^{2}_\epsilon$ from an Inverse-Gamma distribution with a sampling step of $\bm{\theta}$ and $\bm{h}$ through a call to \code{svsample2}.\footnote{In earlier version of the \pkg{stochvol} package, this function was called \code{.svsample}.} This function is a minimal-overhead version of the regular \code{svsample}. It provides the full sampling functionality of the original version but has slightly different default values, a simplified return value structure, and it does not perform costly input checks. Thus, it executes faster and is more suited for repeated calls. The drawback is that it needs to be used with proper care.\footnote{Erroneous or incompatible input values most likely result in run-time errors of compiled code, often implying a segmentation fault and the consecutive abnormal termination of \proglang{R}. This can render debugging tedious.} Note that the current draws of the variables need to be passed to the function through \code{startpara} and \code{startlatent}.
Here is how it goes:
It should be noted that even though \code{svsample2} is considerably faster than \code{svsample}, the cost of interpreting this function in each MCMC iteration is still rather high (which lies in the nature of \proglang{R} as an interpreted language). Thus, as of package version \code{0.8-0}, a single-step \code{update} is also made available at the \proglang{C}/\proglang{C++} level for samplers coded there. For details, please consult the \code{NEWS} file in the package source; for an application case using this approach, see kas-etal:ana.
In the following, we compare the performance of the Bayesian normal linear model with homoskedastic errors from Section (ref) with the Bayesian normal linear model with SV errors from Section (ref) using the \code{exrates} dataset introduced in Section (ref). As a benchmark, we also include the Bayesian normal linear model with GARCH(1,1) errors given through Equation (ref) with
where time index $t=1,\dots,n$. In the second equation, $\tilde y_{t-1}$ denotes the past “residual”, i.e., the $(t-1)$th element of $\tilde \bm{y} = \bm{y} - \bm{X}\bm{\beta}$.
We again use the daily price of 1 EUR in USD from January 3, 2000 to April 4, 2012, denoted by $\bm{p} = (p_1, p_2, \dots, p_n)^\top$. This time however, instead of modeling log returns, we investigate the development of log levels by regression. Let $\bm{y}$ contain the logarithm of all observations except the very first one, and let $\bm{X}$ denote the design matrix containing ones in the first column and lagged log prices in the second, i.e.,
Note that this specification simply corresponds to an AR(1) model for the log prices with three different error specifications: homoskedastic, SV, GARCH(1,1). It is a generalization of directly dealing with demeaned log returns where $\beta_0$ is a priori assumed to be equal to the mean of the log returns and $\beta_1$ is held constant at $1$.
Irrespectively of the error specification, we expect the posterior distribution of $\bm{\beta}$ to spread around $(0,1)^\top$ which corresponds to a random walk. A scatterplot of $\log p_t$ against $\log p_{t+1}$, displayed in Figure (ref), confirms this.
To obtain draws from the posterior distributions for the models with homoskedastic/SV errors, samplers developed in Section (ref) are used. For the GARCH-errors, a simple random walk Metropolis-Hastings sampler is employed. We run each of the three samplers for 100\,000 iterations and discard the first 10\,000 draws as burn-in. Starting values and tuning parameters are set using maximum likelihood estimates obtained through the \proglang{R} package \pkg{fGarch} r:fGarch.
Aiming for comparable results with minimal prior impact on predictive performance, the hyperparameters are chosen to yield vague priors: $\bm{b}_{0} = (0, 0)^\top$, $\bm{B}_{0} = \text{diag}({10^{10}, 10^{10}})$, $c_0 = C_0 = 0.001$, $b_\mu = 0$, $B_\mu=10^4$, $a_0=1$, $b_0 = 1$, $B_{\sigma_\eta}=1$. For the GARCH(1,1) parameter vector $\bm{\alpha} = (\alpha_0, \alpha_1, \alpha_2)^\top$, we pick independent flat priors on $\mathbb{R}^+$ for the components; the initial variance $\sigma_0^2$ is fixed to the empirical residual variance and $\tilde y_0$ is assumed to be zero. Due to the length of the dataset (and its obvious heteroskedasticity), none of these specific choices seem to particularly influence the following analysis. Note, however, that for shorter series or series with less pronounced heteroskedasticity, sensible prior choices are of great importance as the likelihood carries less information in these cases.
The samplers yield slightly different posteriors for $\bm{\beta}$, visualized in Figure (ref). In the top panels, the estimated marginal posterior densities $p(\beta_0|\bm{y})$ and $p(\beta_1|\bm{y})$ are displayed; the bottom panel depicts a scatterplot of draws from the joint posterior $\bm{\beta}|\bm{y}$.
To assess the model fit, mean standardized residuals are depicted in Figure (ref). The model with homoskedastic errors shows deficiencies in terms of pronounced dependence amongst the residuals. This can clearly be seen in the top left panel, where mean standardized residuals are plotted against time. The middle left panel shows the same plot for the model with GARCH errors where this effect is greatly diminished. The bottom left panel pictures that plot for the model with SV errors; here, this effect practically vanishes. Moreover, in the model with homoskedastic errors, the normality assumption about the unconditional error distribution is clearly violated. This can be seen by inspecting the quantile-quantile plot in the top right panel, where observed residuals exhibit much heavier tails than one would expect from a normal distribution. The model with GARCH errors provides a better fit, however, heavy tails are still visible. Standardized residuals from the model with SV errors appear to be approximately normal with only few potential outliers.
Within a Bayesian framework, a natural way of assessing the predictive performance of a given model is through its predictive density (sometimes also referred to as posterior predictive distribution). It is given through
where $\bm{\kappa}$ denotes the vector of all unobservables, i.e., parameters and possible latent variables. Note that for the model with homoskedastic errors, $\bm{\kappa}=(\bm{\beta},{\sigma_\epsilon})^\top$; for the model with SV errors, $\bm{\kappa}=(\bm{\beta}, \bm{\theta}, \bm{h})^\top$; in the GARCH(1,1) case, $\bm{\kappa}=(\bm{\beta},\bm{\alpha},\sigma_0^2,\tilde y_0^2)^\top$. By using the superscript $o$ in $\bm{y}^o_{[1:t]}$, we follow gew-ami:com and denote ex post realizations (observations) for the set of points in time $\{1,2,\dots ,t\}$ of the ex ante random values $\bm{y}_{[1:t]} = (y_1, y_2, \dots,y_t)^\top$. Integration is carried out over ${\bm K}$ which simply stands for the space of all possible values for $\bm{\kappa}$. Equation (ref) can be viewed as the integral of the likelihood function over the joint posterior distribution of the unobservables $\bm{\kappa}$. Thus, it can be interpreted as the predictive density for a future value $y_{t+1}$ after integrating out the uncertainty about $\bm{\kappa}$, conditional on the history $\bm{y}^o_{[1:t]}$.
In the SV errors case, Equation (ref) is a $(n+p+3)$-dimensional integral which cannot be solved analytically. Nevertheless, it may be evaluated at an arbitrary point $x$ through Monte Carlo integration,
where $\bm{\kappa}^{(m)}_{[1:t]}$ stands for the $m$\textsuperscript{th} draw from the respective posterior distribution up to time $t$. If Equation (ref) is evaluated at $x=y_{t+1}^o$, we refer to it as the (one-step-ahead) predictive likelihood at time $t+1$, denoted $P\!L_{t+1}$. Moreover, draws from the posterior predictive distribution can be obtained by simulating values $y_{t+1}^{(m)}$ from the distribution given through the density $p(y_{t+1}|\bm{y}^o_{[1:t]}, \bm{\kappa}^{(m)}_{[1:t]})$, the summands of Equation (ref).
For the model at hand, the predictive density and likelihood can thus be computed through the following
It is worth pointing out that log predictive likelihoods also carry an intrinsic connection to the log marginal likelihood, defined through
This real number corresponds to the logarithm of the normalizing constant which appears in the denominator of Bayes' law and is often used for evaluating model evidence. It can straightforwardly be decomposed into the sum of the logarithms of the one-step-ahead predictive likelihoods: \[ \log M\!L = \log p(\bm{y}^o) = \log \prod_{t=1}^n p(y_t^o|\bm{y}_{[1:t-1]}^o) = \sum_{t=1}^n \log P\!L_{t}. \] Thus, Algorithm (ref) provides a conceptually simple way of computing the marginal likelihood. However, these computations are quite costly in terms of CPU time, as they require an individual model fit for each of the $n$ points in time. On the other hand, due to the embarrassingly parallel nature of the task and because of today's comparably easy access to distributed computing environments, this burden becomes manageable. E.g., the computations for the stochastic volatility analysis in this paper have been conducted in less than one hour, using 25 IBM dx360M3 nodes within a cluster of workstations. Implementation in \proglang{R} was achieved through the packages \pkg{parallel} r:r and \pkg{snow} r:sno.
Cumulative sums of $\log P\!L_t$ also allow for model comparison through cumulative log predictive Bayes factors. Letting $P\!L_t(A)$ denote the predictive likelihood of model $A$ at time $t$, and $P\!L_t(B)$ the corresponding value of model $B$, the cumulative log predictive Bayes factor at time $u$ (and starting point $s$) in favor of model $A$ over model $B$ is straightforwardly given through
When the cumulative log predictive Bayes factor is positive at a given point in time, there is evidence in favor of model $A$, and vice versa. In this context, information up to time $s$ is regarded as prior information, while out-of-sample predictive evaluation starts at time $s+1$. Note that the usual (overall) log Bayes factor is a special case of Equation (ref) for $s=0$ and $u=n$.
In order to reduce prior influence, the first 1000 days are used as a training set only and the evaluation of the predictive distribution starts at $t=1001$, corresponding to December 4, 2003. The homoskedastic model is compared to the model with SV residuals in Figure (ref). In the top panel, the observed series along with the 98% one-day-ahead predictive intervals are displayed. The bottom panel shows the log one-day-ahead predictive likelihood. According to these values, SV errors can handle the inflated volatility during that time substantially better. In the course of 2009, the width of the intervals as well as the predictive likelihoods consolidate again. Figure (ref) also contains a close-up of the one-year time span from September 2008 to August 2009. Both credible intervals are similar at the beginning and at the end of this interval. However, there is a substantial difference in early 2009, where SV intervals become around twice as large compared to the corresponding homoskedastic analogs.
A visually barely distinguishable picture emerges when GARCH(1,1) errors are employed instead of SV errors, thus no separate figure is included in this article.
The top panel of Figure (ref) visualizes the performance of the three models; it depicts observed regression residuals against their (one-day-ahead) predicted distributions. For the purpose of this image, observed regression residuals are simply defined as the deviation of the data from the median of the predicted distribution. It stands out that predictive quantiles arising from the models with heteroskedastic errors exhibit much more flexibility to adapt to the “current state of the world”, while the simpler homoskedastic model barely exhibits this feature.
Generally speaking, there is little difference in predicted residuals until the beginning of 2007. However, during the pre-crisis era (less volatility) and during the actual crisis (more volatility), the models catering for heteroskedasticity perform substantially better. It is interesting to see that predictive intervals for the models with SV errors and GARCH errors are very similar.
In the bottom panel of Figure (ref), the cumulative sums of the log predictive Bayes factors are displayed. The last points plotted equal to $176.52$ in the SV case and $166.37$ in the GARCH case; this provides overwhelming overall evidence in favor of a model catering for heteroskedasticity and “decisive” jef:the3ed respectively “very strong” kas-raf:bay evidence in favor of SV over vanilla GARCH(1,1) with a final cumulative predictive Bayes factor around $25\,000:1$.
It is worth noting that for other exchange rates contained in \code{exrates}, a qualitatively similar picture emerges. For an overview of cumulative log predictive likelihoods, see Table (ref).
As pointed out above, the analysis of log returns can be viewed upon as a special case of the analysis of log levels where $\beta_1 \equiv 1$ is fixed a priori. When doing so, evidence in favor of heteroskedasticity is again striking. Once more, the model with SV scores highest; its demeaned predictive quantiles are almost indistinguishable from the “varying $\beta_1$ case”. Predictive quantiles of the model with GARCH residuals again resemble those of the model with SV residuals very closely. The sum of log predictive likelihoods for $t \in \{1001,1002,\dots,n\}$ is given in the last line of Table (ref). The values are slightly but consistently lower than for the models where $\beta_1$ is estimated from the data. The display of predictive intervals for the log returns is practically identical to Figure (ref) and thus omitted here.
Concluding this somewhat anecdotal prediction exercise, we would like to note that for a thorough and more universal statement concerning the real-world applicability and predictive accuracy of \pkg{stochvol}, further studies with different datasets and a richer variety of competing models, potentially including realized volatility measures, are called for. Such a voluminous undertaking is, however, beyond the scope of this paper.
The aim of this article was to introduce the reader to the functionality of \pkg{stochvol}. This \proglang{R} package provides a fully Bayesian simulation-based approach for statistical inference in stochastic volatility models. The typical application workflow of \pkg{stochvol} was illustrated through the analysis of exchange rate data contained in the package's \code{exrates} dataset. Furthermore, it was shown how the package can be used as a “plug-in” tool for other MCMC samplers. This was illustrated by estimating a Bayesian linear model with SV errors.
In the predictive example, both log levels of exchange rates from EUR to USD and log returns thereof were analyzed. For this dataset, out-of-sample analysis through cumulative predictive Bayes factors clearly showed that modeling regression residuals heteroskedastically substantially improves predictive performance, especially in turbulent times. A direct comparison of SV and vanilla GARCH(1,1) indicated that the former performs better in terms of predictive accuracy.
The author would like to thank Sylvia Fr\"{u}hwirth-Schnatter, Hedibert Freitas Lopes, Karin Dobernig, and two anonymous referees for helpful comments and suggestions.