Bayesian prediction of jumps in large panels of time series data

TL;DR

This paper introduces a Bayesian framework for modeling and predicting jumps in large panels of stock return data, improving forecast accuracy by jointly modeling volatility and jumps.

Contribution

It develops a novel computational approach for univariate stochastic volatility with Poisson jumps and extends it to large panels with dynamic jump intensities.

Findings

Joint modeling of jumps enhances predictive performance.

The proposed method outperforms existing tools in out-of-sample forecasts.

Dynamic factor models capture co-evolving jump intensities.

Abstract

We take a new look at the problem of disentangling the volatility and jumps processes of daily stock returns. We first provide a computational framework for the univariate stochastic volatility model with Poisson-driven jumps that offers a competitive inference alternative to the existing tools. This methodology is then extended to a large set of stocks for which we assume that their unobserved jump intensities co-evolve in time through a dynamic factor model. To evaluate the proposed modelling approach we conduct out-of-sample forecasts and we compare the posterior predictive distributions obtained from the different models. We provide evidence that joint modelling of jumps improves the predictive ability of the stochastic volatility models.

Tables

Table 1. Table 1: Mean, variance and mode of the G ​ a ​ m ​ ( 1 , 50 ) 𝐺 𝑎 𝑚 1 50 Gam(1,50) prior distribution assumed for the jump intensities λ i ​ t subscript 𝜆 𝑖 𝑡 \lambda_{it} in the independent model and for the prior induced by the dynamic factor model.

Prior for ${\lambda }_{it}$	Mean	Variance	Mode
$Gam(1,50)$	$0.020$	$0.00040$	0
Factor model	$0.003$	$0.00004$	$0.0001$

Table 2. Table 2: Mean, variance and mode of the density in ( 24 ) for different values of μ b subscript 𝜇 𝑏 \mu_{b} .

	Mean	Variance	Mode
${\mu }_b=-2.4$	$0.021$	$0.00054$	$0.002$
${\mu }_b=-5$	$0.003$	$0.00004$	$0.0001$
${\mu }_b=-10$	$2\times {10}^{-5}$	$4\times {10}^{-8}$	${10}^{-6}$

Table 3. Table 3: Mean, variance and mode of different gamma priors assumed for the jump intensities λ i ​ t subscript 𝜆 𝑖 𝑡 \lambda_{it} in the case of independent models.

Prior for ${\lambda }_{it}$	Mean	Variance	Mode
A: $Gam(1,250)$	$0.004$	$0.00016$	$0$
B: $Gam(1,50)$	$0.02$	$0.0004$	0
C: $Gam(1.05,52.63)$	$0.02$	$0.00038$	$0.001$
D: $Gam(1.33,66.67)$	$0.02$	$0.00029$	$0.005$
E: $Gam(2,100)$	$0.02$	$0.0002$	$0.01$

Table 4. Table 4: Seconds needed in order the proposed MCMC algorithm to complete one iteration for simulated datasets with different number p 𝑝 p of time series each consisted of 1 , 500 1 500 1,500 log-returns. The last line correspond to the time needed in the case of the real dataset in which the length of the time series with stock-returns varies from 1 , 000 1 000 1,000 to 1 , 917 1 917 1,917 . The timing conducted on a Laptop with a 1.6 GHz Dual-Core Intel Core i5 CPU running R 4.0.0 4.0.0 4.0.0 ( citeR ) .

$p$	Time in seconds
$50$	$0.25$
$100$	$0.52$
$200$	$0.98$
$600$	$2.1$
$571$ (real data)	$2.44$

Equations

$$r_{it}=\exp (h_{it}/2){ϵ}_{it}+\sum _{\kappa =1}^{n_{it}}{\xi }_{it}^{\kappa },{ϵ}_{it}\sim \mathcal{N}(0,1),t=1,\dots ,T$$

$$h_{it}={\mu }_i+{\varphi }_i(h_{i,t-1}-{\mu }_i)+{\sigma }_{i\eta }{\eta }_{it},{\eta }_{it}\sim \mathcal{N}(0,1),t=1,\dots ,T$$

$$n_{it}\sim Poisson({\Delta }_{it}{\lambda }_{it}),$$

$$p(n_{it})=\frac{\Gamma (\delta +n_{it})}{\Gamma (\delta )n_{it}!}{\beta }^{\delta }(1-\beta {)}^{n_{it}},$$

$${\lambda }_{it}={\lambda }^{\ast }(1+e^{-y_{it}}{)}^{-1},$$

$$Y_t$$

$$F_t$$

$$F_0$$

$$q_{m}(n)=\begin{cases}\mbox{\normalsize$\frac{\tilde{p}(n)}{c}$},&\mbox{if }n<m\\ \mbox{\normalsize$\frac{\tilde{p}(m)}{c}\left(\frac{\tilde{p}(m+1)}{\tilde{p}(m)}\right)^{n-m}$},&\mbox{if }n\geq m\end{cases}$$

$$\mathcal{N}_{n_{it}}\left(\xi_{it}\>\Big{|}\>(\mu_{i\xi}\,\frac{e^{h_{it}}}{e^{h_{it}}+n_{it}\sigma_{i\xi}^{2}}+\frac{r_{it}}{e^{h_{it}}+n_{it}\sigma_{i\xi}^{2}}n_{it}\sigma_{i\xi}^{2})\boldsymbol{1}_{n_{it}},\sigma_{i\xi}^{2}I_{n_{it}}-\frac{\sigma_{i\xi}^{4}}{e^{h_{it}}+n_{it}\sigma_{i\xi}^{2}}\boldsymbol{1}_{n_{it}}\boldsymbol{1}_{n_{it}}^{\prime}\right),$$

$$p(\mathrm{x},\theta )\propto \exp \{g(\mathrm{x})\}\mathcal{N}(\mathrm{x}|0,C_{\theta })\pi (\theta )$$

$$p(\mathrm{x},\theta ,\mathrm{u})\propto \exp \{g(\mathrm{x})\}{\mathcal{N}}_d(\mathrm{x}|0,C_{\theta }){\mathcal{N}}_d(\mathrm{u}|\mathrm{x},(\delta /2)I_d)\pi (\theta ).$$

$$q({\mathrm{x}}^{\prime },{\theta }^{\prime }|\mathrm{x},\theta ,\mathrm{u})=q({\mathrm{x}}^{\prime }|\mathrm{x},{\theta }^{\prime },\mathrm{u})q({\theta }^{\prime }|\theta ).$$

$$q({\mathrm{x}}^{\prime }|\mathrm{x},{\theta }^{\prime },\mathrm{u})\propto$$

$$\propto {\mathcal{N}}_d\left({\mathrm{x}}^{\prime }|\frac{2}{\delta }{Q}_{{\theta }^{\prime }}(\mathrm{u}+\frac{\delta }{2}\nabla g(\mathrm{x})),Q_{{\theta }^{\prime }}\right),$$

$$\mathrm{z}\equiv \mathrm{u}+(\delta /2)\nabla g(\mathrm{x})\sim {\mathcal{N}}_d(\mathrm{x}+(\delta /2)\nabla g(\mathrm{x}),(\delta /2)I_d).$$

$$p(\mathrm{x},\mathrm{z},\theta )\propto \exp \{g(\mathrm{x})\}{\mathcal{N}}_d(\mathrm{x}|0,C_{\theta }){\mathcal{N}}_d(\mathrm{z}|\mathrm{x}+(\delta /2)\nabla g(\mathrm{x}),(\delta /2)I_d)\pi (\theta ),$$

$$q({\mathrm{x}}^{\prime }|{\theta }^{\prime },\mathrm{z})={\mathcal{N}}_d\left({\mathrm{x}}^{\prime }|\frac{2}{\delta }{Q}_{{\theta }^{\prime }}\mathrm{z},Q_{{\theta }^{\prime }}\right)=\frac{1}{\mathcal{Z}(\mathrm{z},{\theta }^{\prime })}{\mathcal{N}}_d({\mathrm{x}}^{\prime }|0,C_{{\theta }^{\prime }}){\mathcal{N}}_d(\mathrm{z}|{\mathrm{x}}^{\prime },C_{{\theta }^{\prime }}),$$

$$q({\mathrm{x}}^{\prime },{\theta }^{\prime }|\mathrm{x},\theta ,\mathrm{z})=q({\mathrm{x}}^{\prime }|{\theta }^{\prime },\mathrm{z})q({\theta }^{\prime }|\theta )$$

$$\exp \left\{g({\mathrm{x}}^{\prime })-g(\mathrm{x})+\nu (\mathrm{z},{\mathrm{x}}^{\prime })-\nu (\mathrm{z},\mathrm{x})\right\}\times \frac{\mathcal{Z}(\mathrm{z},{\theta }^{\prime })\pi ({\theta }^{\prime })q(\theta |{\theta }^{\prime })}{\mathcal{Z}(\mathrm{z},\theta )\pi (\theta )q({\theta }^{\prime }|\theta )}$$

$$\log (p(R_{T+1:T+j}|R_{1:T},{\mathbb{M}}_{\mathrm{svj}}))-\log (p(R_{T+1:T+j}|R_{1:T},{\mathbb{M}}_{\mathrm{sv}})),j=1,\dots ,\ell ,$$

$$p(R_{T+1:T+j}|R_{1:T})=\prod _{i=1}^p\prod _{t=T+1}^{T+j}p(r_{it}|r_{i,1:t-1}).$$

$$(u_{it}-l_{it})+\frac{2}{\alpha }(l_{it}-r_{it})\text{mathbbm}1\{r_{it}<l_{it}\}+\frac{2}{\alpha }(r_{it}-u_{it})\text{mathbbm}1\{r_{it}>u_{it}\}.$$

$$p(h_{i,0:t}|r_{i,1:t})=\frac{p(h_{i,0:t-1}|r_{i,1:t-1})p(r_{it}|h_{it})p(h_{it}|h_{i,t-1})}{p(r_{it}|r_{i,1:t-1})},$$

$$g_j(H_{0:T+\ell },F_{0:T+\ell })$$

$$=p(H_{0:T+j},F_{0:T+j}|R_{1:T+j})\prod _{t=T+j+1}^{T+\ell }p(H_t|H_{t-1})p(F_t|F_{t-1}),$$

$$g_j^{\ast }(H_{0:T+\ell },F_{0:T+\ell })=p(H_{0:T+j},F_{0:T+j},R_{1:T+j})\prod _{t=T+j+1}^{T+\ell }p(H_t|H_{t-1})p(F_t|F_{t-1}),j<\ell ,$$

$${\Omega }_{T+j}^s={\Omega }_{T+j-1}^s\frac{g_j^{\ast }(H_{0:T+\ell }^s,F_{0:T+\ell }^s)}{g_{j-1}^{\ast }(H_{0:T+\ell }^s,F_{0:T+\ell }^s)}={\Omega }_{T+j-1}^{s}p(R_{T+j}|H_{T+j}^s,F_{T+j}^s),j\ge 2,$$

$$p(R_t|H_t^s,F_t^s)$$

$$=\prod _{i=1}^p\sum _{n_{it}=0}^{\infty }\mathcal{N}(r_{it}|n_{it}{\mu }_{i\xi },\exp (h_{it}^s)+n_{it}{\sigma }_{i\xi }^2)p(n_{it}|{\Delta }_{it}{\lambda }_{it}^s).$$

$$p(r_{it}|R_{1:t-1})$$

$$\times p(h_{i,0:t-1},F_{0:t-1}|R_{1:t-1})\mathrm{d}{h}_{i,0:t}\mathrm{d}{F}_{0:t}.$$

$${\omega }_{it}^s={\omega }_{i,t-1}^{s}p(r_{it}|h_{it}^s,F_t^s),$$

Full text

Keywords— dynamic factor model; forecasting stock returns; Markov chain Monte Carlo; stochastic volatility with jumps; sequential Monte Carlo

1 Introduction

It has been recognised in the financial literature that jumps in asset returns occur clustered in time and affect several stock markets within a few hours or days, see for example Aït-Sahalia et al. (2015). We work with a large panel of stocks from several European markets and we aim to identify and predict joint tail risk expressed as probabilities of jumps in their daily returns. Our modelling assumptions are based on the well-known paradigms of stochastic volatility (SV) models (Taylor, 1982) combined with Poisson-driven jumps (Andersen et al., 2002). The estimation of SV models is preferably conducted by using likelihood-based approaches (Harvey et al., 1994) while Bayesian methods, in particular, are notably desirable since they deal efficiently with likelihood intractability problems; see for example Jarquier et al. (1994), Chib et al. (2002) and Kastner and Frühwirth-Schnatter (2014) for detailed discussions. This motivates us to utilize the Bayesian approach both for modelling and inference purposes. Our purpose is to provide a general modelling approach for the time and cross-sectional dependence of jumps in multiple time series. The resulting Bayesian hierarchical models require careful prior specifications and sophisticated, modern Markov chain Monte Carlo (MCMC) inference implementation strategies. To that end, we build the necessary algorithmic framework by developing robust computational methods. To evaluate the proposed model we focus on forecasting future stock returns. This is a problem of primary interest in financial statistics since it is closely related with risk management, portfolio allocation and asset pricing; see for example Aguilar and West (2000), Rapach and Zhou (2013) and Clements and Liao (2017) for more detailed discussions. We test the proposed methods in an out-of-sample forecasting scenario with real data.

Our first point of departure is the univariate SV model with jumps. We provide a new efficient MCMC algorithm combined with careful prior specification that improves upon existing MCMC strategies. Next, we extend this approach by utilising the panel structure of stock returns so that unobserved jump intensities are assumed to propagate over time through a dynamic factor model. The resulting Bayesian hierarchical model can be viewed as a generalisation of the Bates (1996), univariate SV model with jumps since it can be obtained as a special case by just assuming that jump intensities are independent across time and stocks. Furthermore, our computational techniques are relevant to other scientific areas such as, for example, neuroscience, whereas simultaneous recordings of neural spikes are often modelled through latent factor models (Buesing et al., 2014).

To illustrate our motivation, consider the daily returns from the stocks of the STOXX Europe $600$ Index over a period $10/1/2007$–$11/6/2014$. The top panel of Figure 1 has been created by empirically identifying a jump when the return exceeds three standard deviations, where the estimators of mean and variance of each series were robustly estimated as the median and the robust scale estimator of (Rousseeuw and Croux, 1993) respectively; see the supplementary material for the details of the empirical method used for outlier detection. The middle panel depicts a summary statistic of a SV model with jumps estimated separately for each stock return series: each dot denotes a stock return in which that particular day the probability of a jump in the series has been estimated to be greater than $.5$; by choosing $.5$ as threshold we avoid to over- or underestimate the number of jumps for each stock at each day, see the supplementary material for more details. Compared with the top panel, it provides a more sparse jump process indicating a successful separation of jumps from persistent moves modelled by the volatility process, a topic discussed in detail by, for example, Aït-Sahalia (2004). The notable feature that inspired this work is that clustering and inter-dependence of jump intensities is evident in the middle panel and, interestingly, days in which jumps have been simultaneously identified in a large number of stocks coincide with days of events that affected the financial markets worldwide as pointed out at the x-axis of the plots. Thus, one may attempt a joint Markovian modelling of jump intensities across time which might reveal some predictability of jumps. Indeed, the results of our proposed joint modelling formulation provides strong evidence for better out-of-sample predictability with estimated jump probabilities higher than $.5$ depicted in the lower panel of Figure 1.

We propose a modelling approach for the jump processes of SV models in which the time and cross-sectional dependence of the jump intensities are driven by a latent, common across stocks, dynamic factor model. We carefully specify informative prior distributions to the parameters of the factor model so that the implied priors for the jump intensities in each stock are in concordance with the well-known (Eraker et al., 2003) a priori expectation for one jump every few months. To avoid identifiability problems of the factor models extra assumptions for the form of the loadings matrix are needed. Nevertheless, it is recognised (Ghosh and Dunson, 2009) that MCMC algorithms developed under identifiability constraints exhibit slow convergence to the target distribution of interest. To improve the mixing properties of the proposed MCMC algorithm we follow Bhattacharya and Dunson (2011) and we do not impose the usual (see, e.g., Aguilar and West (2000)) identifiability constraints on the factor loadings parameters. We also note that there is a considerable amount of literature dedicated to the modelling of time series count data; see for example Fokianos et al. (2009), Jung et al. (2011), Pillow and Scott (2012) and Buesing et al. (2014). Our modelling perspective is applied on unobserved count data and, thus, the aim of the proposed latent factor model is twofold. First, since a few jumps are expected to occur during the observation period we borrow information across all the stocks in order to learn the unobserved processes that drive the evolution of jumps across time and stocks. Second, we exploit the dimension reduction achieved by the latent factors to preserve the scalability of our method.

Our Bayesian inference implementation strategy is based on an MCMC algorithm that alternates sampling of the latent volatility process and its parameters with sampling of the jump process and its parameters. Both these steps require drawing of the high-dimensional paths of the latent volatilities and factors from their full conditional distributions. By noting that the target distributions are proportional to the product of intractable, non-linear likelihood functions with Gaussian priors, we utilize the sampler proposed by Titsias and Papaspiliopoulos (2018) to simultaneously draw the whole path of each latent process at each MCMC iteration. This is an important feature of our methods for the following reasons. First updating each state of a volatility process separately results in Markov chain samplers with slow mixing (Shephard and Kim, 1994). Second, although Chib et al. (2002) and Nakajima and Omori (2009) use a mixture of Gaussian distributions approximation to also simultaneously update the whole volatility path, their methods rely on an importance sampling step which is quite problematic, see Johannes et al. (2009) for a detailed discussion and Section 4.3 of this paper for a simulation that illustrates this issue.

To assess the predictive performance of the different models we compare the corresponding posterior predictive distributions by utilizing proper scoring rules (Gneiting and Raftery, 2007). We provide a full-fledged quantitative evaluation of the obtained forecasts by calculating logarithmic scores, such as predictive log Bayes factors (Geweke and Amisano, 2010), interval and continuous ranked probability scores (Gneiting and Katzfuss, 2014) as well as root mean squared errors. We estimate the posterior predictive distributions by employing sequential Monte Carlo methods such as the particle filters (Chopin and Papaspiliopoulos, 2020a) and the annealed importance sampling (Neal, 2001) algorithms.

The structure of the remaining of the paper is as follows. In Section 2 we provide a review of the literature related to univariate and multivariate SV models in econometrics as well as the modelling of jumps in financial applications. Section 3 presents our proposed modelling framework and Section 4 describes the methods that we develop to conduct Bayesian inference for the proposed model. In Section 5 we present the computational techniques that we use to assess the predictive performance of the proposed model. Section 6 presents results and insights from the application of our methods on the real dataset and Section 7 concludes with a small discussion.

2 Related work

There are two main classes of models that have been considered in the literature for the modelling of financial returns. The first class is based on observation-driven models where the time evolution of the variance of the returns is modelled by using past observations and variances. These models are known as autoregressive conditional heteroscedasticity (ARCH) models and have been developed by Engle (1982). A very popular parametrization, the generalized ARCH model, introduced later by Bollerslev (1986) while Bollerslev (1987) and Harvey and Chakravarty (2008) presented modifications that account for the heavy tails of the distributions of the observed returns; see also Harvey (2013) and references therein. Creal et al. (2011) extended this modelling approach in the multivariate case; for the Bayesian analysis of these models see, for example, Vrontos et al. (2000) and Vrontos et al. (2003). The second class is known as parameter-driven models where the time evolution of the variance of the returns is modelled by employing latent stochastic processes. The models that belong to this class are the SV models which constitute the base of our proposed modelling approach.

The SV model proposed by Taylor (1982) in order to account for the time evolution of the volatilities in financial time series while Hull and White (1987) and Chesney and Scott (1989) have employed the SV model as discretization of continuous time SV diffusion processes. As mentioned earlier some early works on the Bayesian estimation of SV models include, but are not limited to, the papers presented by Jarquier et al. (1994), Kim et al. (1998) and Chib et al. (2002). Omori et al. (2007) dealt with efficient Bayesian estimation of the SV model by taking into account the correlation between the observation error and the error of the stochastic variance process; see also Silva et al. (2006) and Delatola et al. (2011) for other extensions of the basic SV model. A more recent contribution on the Bayesian analysis of univariate SV models has been offered by Kastner and Frühwirth-Schnatter (2014) while Zhang et al. (2020) developed a new class of SV models with autoregressive moving average errors. Considering a jump component to model the empirically observed fat tails of stock returns proposed, in his pioneering paper, by Merton (1976) while Bates (1996) and Andersen et al. (2002) combined jump diffusions with SV models. Bayesian estimation of SV models with jumps has been considered, among others, by Chib et al. (2002), Nakajima and Omori (2009) and Eraker et al. (2003) and Johannes et al. (2009) where the latter two papers account for jumps in the volatility process as well; see also Bandi and Renò (2016) for more recent research on the occurrence of contemporaneous jumps in the returns and volatilities.

To conduct joint modelling of the observed financial time series, multivariate SV models have been considered more than twenty years ago; see for example Harvey et al. (1994) and Asai et al. (2006) for a review. The existing multivariate SV models assume either constant correlations over time or some form of dynamic correlation modelling through factor models with factors being independent univariate stochastic volatility models. The latter approach has been considered, for example, by Chib et al. (2006), Kastner et al. (2017) and Kastner (2019). In a similar framework Chan and Eisenstat (2018) and Kastner and Huber (2020) deal with Bayesian estimation of time varying parameter vector autoregressive models with stochastic volatility. An alternative multivariate SV modelling approach proposed recently by Dellaportas et al. (2015) where all the elements of the covariance matrix of the observed time series evolve over time. More closely related to the present paper, Bollerslev et al. (2008), Jacod et al. (2009) and Clements and Liao (2017) study the presence of common or not common jumps in multiple time series. Aït-Sahalia et al. (2015) model the dependencies of jumps by combining the SV model with a multivariate Hawkes process (Hawkes, 1971) that models the propagation of jumps over time and across assets by introducing a feedback from the jumps to their intensities and back. They employ the generalised version method of moments to estimate the parameters of the model. We close this literature review section by noting that tools for the identification of non-observed events that affect the evolution of time series have been also developed by Hamilton (1989) and Billio et al. (2012) that deal with the problem of identifying structural changes in time series of economic indices.

3 Latent jump modelling

3.1 Notation

We denote the univariate Gaussian distribution with mean $m$ and variance $s^{2}$ by $\mathcal{N}(m,s^{2})$, and its density evaluated at $x$ by $\mathcal{N}(x|m,s^{2})$; $\mathcal{N}_{d}(X|M,S)$ denotes the density of the $d$-variate normal distribution with mean $M$ and covariance matrix $S$ evaluated at $X$; $Gam(\alpha,\beta)$ and $IGam(\alpha,\beta)$ denote the gamma, with mean $\alpha/\beta$, and inverse gamma, defined as the reciprocal of gamma, distributions respectively; $U(\alpha,\beta)$ the uniform distribution on interval $(\alpha,\beta)$; $\Gamma(\cdot)$ denotes the gamma function. Index $i$ is used for stocks and index $t$ for times, e.g., $r_{it}$ is the $i$th stock return at day $t$; when we introduce factor models index $k$ denotes factor. Upper case letters denote vectors and matrices, and lower case letters denote scalars. For a collection of variables indexed by time or stock and denoted by a lower case character, the corresponding upper case character denotes their vector or matrix representation, e.g. $R$ denotes the matrix whose elements are $r_{it}$; $R_{i}=(r_{i1},\ldots,r_{iT})$, and $R_{t}=(r_{1t},\ldots,r_{pt})$ are then generic rows and columns of $R$; all vectors are understood as column vectors. A subscript, of the form $(t+1):(t+\ell)$, where $\ell$ is a positive integer, indicates the collection of vectors with subscripts $t+1,\ldots,t+\ell$; e.g. $R_{t+1:t+\ell}$ is the collection of vectors $R_{t+1},\ldots,R_{t+\ell}$. The transpose of a vector or matrix is denoted by a prime; e.g. $R^{\top}$ is the transpose of $R$; $\mathbbm{1}\{\mathord{\cdot}\}$ denotes an indicator function that takes value $1$ if the event in brackets is true and [math] otherwise; for a function $f(x)$ we set $\nabla\eqqcolon(\partial f/\partial x_{1},\ldots,\partial f/\partial x_{d})$.

3.2 Stochastic volatility with jumps

We model stock returns as

[TABLE]

[TABLE]

with $h_{i0}\sim\mathcal{N}(\mu_{i},\sigma_{i\eta}^{2}/(1-\phi_{i}^{2}))$ and $\epsilon_{it}$ and $\eta_{it}$ are independent. The number of jumps of the $i$th stock at time $t$ follow a Poisson distribution,

[TABLE]

where $\lambda_{it}$ denotes the corresponding jump intensity, $\Delta_{it}$ is a time-increment associated to each return and $\xi_{it}^{\kappa}\sim\mathcal{N}(\mu_{i\xi},\sigma^{2}_{i\xi})$ is the $\kappa$th jump size. We explicitly take into account that stock returns are computed over varying time increments $\Delta_{it}$, such as in-between weekends and holidays, hence $\lambda_{it}$ is the daily jump intensity.

We model the parameters $\mu_{i},\phi_{i}$ and $\sigma_{i\eta}^{2}$ of each log-volatility process and the parameters $\mu_{i\xi}$ and $\sigma^{2}_{i\xi}$ of the jump sizes as independent across stocks. For $\mu_{i}$, $\phi_{i}$ and $\sigma^{2}_{i\eta}$ we follow the related literature (Kim et al. (1998), Omori et al. (2007), Kastner and Frühwirth-Schnatter (2014)) and take $\mu_{i}\sim\mathcal{N}(0,10)$, $(\phi_{i}+1)/2\sim Beta(20,1.5)$ and $\sigma^{2}_{i\eta}\sim Gam(1/2,1/2)$.

For $\mu_{i\xi}$ and $\sigma_{i\xi}^{2}$ we choose proper priors. In fact, improper priors will lead to improper posterior. This is due to the fact that according to the model there is positive probability that there are no jumps for the whole period of time, i.e., the event $n_{it}=0$ for all $t$, has positive prior probability, and in such an event, $\mu_{i\xi}$ and $\sigma_{i\xi}^{2}$ are unidentifiable in the likelihood. Instead, we take $\sigma_{i\xi}^{2}\sim IGam(3,\textrm{range}_{i}^{2}/18)$, where $\textrm{range}_{i}=\max_{t}\{r_{it}\}-\min_{t}\{r_{it}\}$, and $\mu_{i\xi}\sim\mathcal{N}(0,5\textrm{range}_{i}^{2})$. In this specification, $\mathrm{E}[\sigma^{2}_{i\xi}]=\textrm{range}_{i}^{2}/36$, hence we match the variance of jump sizes with a quantity that relates to the sample variance of returns. A same reasoning is followed for setting the prior variance of $\mu_{i\xi}$ and we have found that the choice of the multiplier, here taken 5, is not crucial in the analysis.

3.2.1 Modelling jumps independently across stocks and time

By modelling $\Lambda_{i}$ as independent random vectors, we obtain $p$ independent univariate SV with jumps models. Additionally, each vector can be taken as one of independent elements, $\lambda_{it}\sim Gam(\delta,c)$, in which case one obtains the models described, among others, by Bates (1996), Chib et al. (2002) and Eraker et al. (2003). It follows that $n_{it}$ has, marginally with respect to $\lambda_{it}$, a negative binomial distribution with density

[TABLE]

where $\beta=c/(c+\Delta_{it})$. The mean of the distribution is $\delta(1-\beta)/\beta$ and the variance is $\delta(1-\beta)/\beta^{2}.$ Taking $\delta\leq 1$ ensures that the density is monotonically decreasing. This is a feature we are interested in: jumps are meant to capture infrequently large price movements, and a priori we expect that with highest probability there is no jump and it is more likely to have less than more jumps, if any. We can choose $c$ such that the probability of no jump at a given time increment $\Delta_{it}$, is at least $\gamma\in(0,1)$, a user-specified threshold. Taking for example $\delta=1$ and $c=50$ results in $\gamma=0.98$ with $\mathrm{E}(\lambda_{it})=0.02$. These choices are consistent with the prior expectation, in financial applications, for one jump every few months; see for example Eraker et al. (2003) for a more detailed discussion.

3.3 Joint modelling of jumps

To capture the dependence of the jumps over time and across stocks we model $\lambda_{it}$ by using a dynamic factor model which is specified as follows. We first transform $\lambda_{it}$ to $y_{it}$ via

[TABLE]

which implies that $\lambda_{it}<\lambda^{*}$. The parameter $\lambda^{*}$ is, thus, the upper bound for the jump intensities off all the stocks at any time point. This is an important feature of our modelling approach since it allows to incorporate the necessary (Chib et al., 2002; Eraker et al., 2003) prior information about the expected number of jumps which are typically described as large but infrequent moves. Furthermore, as discussed below, by carefully choosing this upper bound we are able to construct a prior for the jump intensities which is similar to the $Gam(1,50)$ distributions that we used in the case of univariate SV models with jumps. We also note that this type of information is not necessary only in Bayesian applications; see for example Honore (1998) for a similar discussion related to maximum likelihood estimation of models with jumps. Then, we model the time-dependence of jumps via $K$ latent factors $F_{t}$ modelled as independent autoregressive processes. The cross-sectional dependence of jumps is captured by a $p\times K$ matrix of factor loadings $W$ and a $p\times 1$ vector $B$. The joint model specification is

[TABLE]

where $A$ is a $K\times K$ diagonal matrix with elements $\alpha_{1},\ldots,\alpha_{K}$ and $\Sigma_{F}$ is a $K\times K$ diagonal matrix with elements $1/(1-\alpha_{1}^{2}),\ldots,1/(1-\alpha_{K}^{2})$.

It is known that parameters of latent factor models are only identifiable up to certain transformations, such as orthogonal rotations and sign changes (Aguilar and West, 2000). However, latent factor models are also used as low-rank predictive models in which case the out-of-sample prediction is of prior importance; see Bhattacharya and Dunson (2011). This is the perspective we adopt here, where we try to borrow strength from the information in large panels of stocks in order to predict with higher accuracy individual jumps by using a parsimonious factor model. Therefore, our prior specification imposes no identifiability constraints on the parameters of the factor model and are taken as $w_{ik}\overset{iid}{\sim}\mathcal{N}(0,\sigma^{2}_{w})$, $b_{i}\overset{iid}{\sim}\mathcal{N}(\mu_{b},\sigma^{2}_{b})$, and $\alpha_{k}\overset{iid}{\sim}U(-1,1)$.

The specification of the hyperparameters of the priors assigned on the parameters $A,B$ and $W$ require extra care because they affect, through (5), the imposed prior on $\lambda_{it}$. The fact that our modelling formulation of Section 3.2.1 was based on informative priors $Gam(1,50)$ provides a way to elicit informative prior specification of the hyperparameters $\sigma^{2}_{w},\sigma^{2}_{b}$ and $\mu_{b}$. We first set $\lambda^{*}=0.15$ which is the $0.0005$ upper quantile of the $Gam(1,50)$ density. As discussed earlier, this an important characteristic of our modelling approach since it penalizes the estimation of a large number of jumps which is not consistent with the prior information for one or two jumps every $100$ days. A small sensitivity analysis, not reported here, confirms that choosing a value considerably greater than $0.15$ does not allow the efficient separation of the jumps from the underlying volatility process of the stock daily returns. Moreover, by using values notably smaller than $0.15$ we estimate quite low jump intensities and this can result in less accurate prediction of future stock returns; in the supplementary material we compare predictions obtained under different prior assumptions for the expected number of jumps. We have concluded that a value for $\lambda^{*}$ between $0.1$ and $0.2$ does not have a substantial effect on the results presented in the rest of the paper. Then, we choose the hyperparameters $\sigma^{2}_{w}$, $\mu_{b}$ and $\sigma^{2}_{b}$ to be such that the mean, the variance and the mode for the resulting prior on each $\lambda_{it}$ in (5) are comparable with the corresponding quantities of the $Gam(1,50)$ distribution. Based on calculations presented in the supplementary material we set $\sigma^{2}_{w}=0.5$, $\sigma^{2}_{b}=1$ and $\mu_{b}=-5$. In Table 1 we display the quantities of interest in the two prior distributions. To evaluate the differences between the two priors we also examine the impact of different choices for their hyperparameters in the forecasts of future observations. See the supplementary material for the related results.

4 Sampling from the posterior of interest

Let $\Theta_{i}=(\mu_{i},\phi_{i},\sigma_{i\eta}^{2},\mu_{i\xi}$,$\sigma^{2}_{i\xi},b_{i})$. Our interest lies on the joint posterior of parameters and latent states $p(H,N,\Xi,W,\Theta,F,A|R)$. To draw samples from the posterior of interest we design an MCMC algorithm which proceeds as follows. At each MCMC iteration we perform stock-specific updates of $H_{i}$, $N_{i}$, $\Xi_{i}$, $W_{i}$, $\Theta_{i}$ and then we update the latent factors $F$ and their parameters $A$.

An important characteristic of the designed MCMC algorithm is that before updating the log-volatilities $H_{i}$ and the number of jumps $N_{i}$ we integrate out the jump sizes $\Xi_{i}$. This results in an efficient simultaneous update of the log-volatilities vector and direct sampling of the full conditional of the number of jumps $N_{i}$. Algorithm 1 summarizes the steps of the MCMC sampler. The algorithm presents how the full conditional densities of vectors of parameters are sampled by exploiting the conditional independence structure of the hierarchical model defined by equations (1)–(3) and (5)–(8). In Sections 4.1 and 4.2 we will describe in detail the more elaborate steps $4,6,7$ and $12$ of Algorithm 1 whereas the details of the other sampling steps are given in the supplementary material.

We also emphasize that by switching off certain steps of Algorithm 1 we obtain a novel MCMC sampler for existing models. More precisely, Bayesian inference for univariate SV with jumps models can be conducted if the step in the $10$th line is skipped and the step in the $12$th line replaced with the step of drawing the independent jump intensities of the model. Switching off the steps in lines $6$–$10$ and $12$ results in an MCMC algorithm for univariate SV without jumps models. By removing the steps in lines $4$–$9$ of Algorithm 1 we obtain a sampler for a Cox model (Cox, 1955) with intensities driven by latent factors.

4.1 Sampling the number of jumps and jump sizes

To sample each $n_{it}$ we first integrate out the jump sizes $\xi_{it}$ and then we construct a rejection sampling algorithm to draw from $p(n_{it}|h_{it},\mu_{i\xi},\sigma^{2}_{i\xi},b_{i},W_{i},F_{t},r_{it})$. This algorithm is based on Proposition 1. Notice also that a discrete distribution with pmf $p(n)$ is called log-concave if $p(n)^{2}\geq p(n-1)p(n+1)$ is satisfied for all $n$; see for example Devroye (1986) for more details. For the remaining of this subsection the unnormalized $p(n_{it}|h_{it},\mu_{i\xi},\sigma^{2}_{i\xi},b_{i},W_{i},F_{t},r_{it})$ is denoted by $\tilde{p}(n)$.

Proposition 1**.**

For each $i=1,\ldots,p$ and $t=1,\ldots,T$ the discrete distribution with probability mass function $p(n_{it}|h_{it},\mu_{i\xi},\sigma^{2}_{i\xi},b_{i},W_{i},F_{t},r_{it})$ is log-concave for $n_{it}\geq 1$.

Proof.

See the supplementary material. ∎

Due to Proposition 1, the target distribution is a unimodal distribution (Devroye, 1986). Let $m$ be an integer at the right of the mode of the distribution. We define the probability mass function (pmf)

[TABLE]

with $c=\sum_{n=0}^{m-1}\tilde{p}(n)+\tilde{p}(m)^{2}/(\tilde{p}(m)-\tilde{p}(m+1))$. The pmf $q_{m}(n)$ is proportional to $\tilde{p}(n)$ for $n<m$ and proportional to the geometric pmf with parameter, $(\tilde{p}(m)-\tilde{p}(m+1))/\tilde{p}(m)$, for $n\geq m$. By noting that if a random variable $X$ follows the exponential distribution with parameter, $\log(\tilde{p}(m))-\log(\tilde{p}(m+1))$, then $\lfloor X\rfloor$ follows the geometric distribution with parameter $(\tilde{p}(m)-\tilde{p}(m+1))/\tilde{p}(m)$, we observe that $q_{m}(n)$, for $n\geq m$, corresponds to an exponential density that goes through the points $(m,\tilde{p}(m))$ and $(m+1,\tilde{p}(m+1))$. Log-concavity of the distribution implies that $cq_{m}(n)\geq\tilde{p}(n)$, for $n\geq m$, since any straight line that goes through the points $(m,\log(\tilde{p}(m)))$ and $(m+1,\log(\tilde{p}(m+1)))$ bounds from above the logarithm of $\tilde{p}(n)$. See Figure 2 for an illustration. A rejection algorithm to sample from the distribution with pmf $p(n_{it}|h_{it},\mu_{i\xi},\sigma^{2}_{i\xi},b_{i},W_{i},F_{t},r_{it})$ proceeds as summarized by Algorithm 2. Note that Algorithm 2 implies the evaluation of $\tilde{p}(0),\dots,\tilde{p}(m+1)$ at each iteration of the proposed MCMC algorithm.

4.1.1 Jump sizes

Drawing the jump sizes $\xi_{it}$ given the number of jumps $n_{it}$ is based on the following proposition. Let $I_{n}$ and $\boldsymbol{1}_{n}$ be the $n\times n$ identity matrix and an $n$-dimensional vector with ones respectively.

Proposition 2**.**

According to the model defined by equations (1) and (2) we have that $p(\xi_{it}|n_{it},h_{it},\mu_{i\xi},\sigma^{2}_{i\xi},r_{it})$ is equal to the Gaussian density

[TABLE]

for each $i=1,\ldots,p$.

Proof.

See the supplementary material. ∎

4.2 Log-volatilities and factors: MCMC for latent Gaussian models

In the $4$th and $12$th lines of Algorithm 1 we sample the paths of the log-volatilities $H_{i}$ and factors $F$ respectively. Moreover, each path is updated jointly with the parameters of the corresponding autoregressive models defined by equations (2) and (8). These steps correspond to the joint sampling of latent paths and parameters of latent Gaussian models. A latent Gaussian model is defined in terms of a non-Gaussian likelihood and a latent Gaussian field controlled by some parameters. By denoting with $\mathrm{X}$ a $d$-dimensional random vector that corresponds either to a latent log-volatility or to a factor path, with $\mathrm{x}$ its realization and with $\theta$ any parameters we target the density

[TABLE]

where $g(\mathrm{x})$ is the log-likelihood function, $C_{\theta}$ denotes the covariance matrix of the latent Gaussian field, parametrized in terms of $\theta$ and $\pi(\theta)$ denotes the density of the prior for $\theta$. We have assumed, without loss of the generality, that the latent Gaussian field has zero mean. We also note that only throughout the present Section we violate our notation and we use lower-case letters to denote the realizations of a random vector; the dependence on any data has been suppressed.

To draw samples from the target in (9) two approaches have dominated in the literature. One of them is to alternate sampling from the conditionals $p(\mathrm{x}|\theta)$ and $p(\theta|\mathrm{x})$. Sampling from $p(\mathrm{x}|\theta)$ is quite challenging since usually corresponds to a high dimensional distribution. To draw the desired samples several well-established methods, based on the Metropolis-Hastings algorithm, have been used, see e.g. Roberts and Stramer (2002), Murray et al. (2010), Girolami and Calderhead (2011) and Cotter et al. (2013); in the latter approach the proposal distribution utilizes information both from the likelihood and the Gaussian prior of the latent states and this is an important feature of the proposed sampler (Titsias and Papaspiliopoulos, 2018). Sampling from $p(\theta|\mathrm{x})$ is also challenging since the parameters $\theta$ are highly correlated with the Gaussian latent states of the model. See, for example, Papaspiliopoulos et al. (2007), Murray and Adams (2010) and Yu and Meng (2011) for efficient methods that have been developed for the update of $\theta$, conditional on the latent Gaussian filed, in a more general class of Bayesian hierarchical models. In the second approach, joint sampling of the latent states $\mathrm{x}$ and the parameters $\theta$ is conducted. This perspective is considered, among others, by Knorr-Held and Rue (2002) where the latent states and the parameters are updated in, carefully selected, blocks and by Titsias and Papaspiliopoulos (2018) where a sampler which utilizes auxiliary variables to perform updates that are both prior and likelihood informed has been constructed.

In the present paper we follow the latter approach and we update jointly the latent states and the parameters of a given model by combining the methods proposed by Yu and Meng (2011) and Titsias and Papaspiliopoulos (2018). Following Yu and Meng (2011) we interweave the so-called (Papaspiliopoulos et al., 2007) centred and non-centred parametrization of the latent states and by utilizing the algorithm developed by Titsias and Papaspiliopoulos (2018) we sample jointly $\mathrm{x}$ and $\theta$ under each one of the two parametrizations. We show that by utilizing the proposed Metropolis-Hastings step we find efficiency gains relative to alternative samplers that can be used to update the latent log-volatilities, the factor paths and their parameters. In the rest of the Section we present the proposed methodology in the case of a general latent Gaussian model and in the supplementary material we provide the details for the application of the method in order to sample the latent factors $F$ as well as each log-volatility path $H_{i}$ and their parameters.

4.2.1 Metropolis-Hastings for the joint update of latent paths and parameters

To construct a Metropolis-Hastings step in order to sample jointly $\mathrm{x}$ and $\theta$ we work as in Titsias and Papaspiliopoulos (2018) and we introduce auxiliary variables $\mathrm{u}\sim\mathcal{N}_{d}(\mathrm{x},(\delta/2)I_{d})$. We consider, thus, the expanded target

[TABLE]

Samples from the expanded target can be drawn by first drawing $\mathrm{u}$ from $p(\mathrm{u}|\mathrm{x})=\mathcal{N}_{d}(\mathrm{u}|\mathrm{x},(\delta/2)I_{d})$ and then employing the Metropolis-Hastings algorithm to draw from the intractable $p(\mathrm{x},\theta|\mathrm{u})$. To perform the latter update we utilize a proposal distribution of the form

[TABLE]

To construct $q(\mathrm{x}^{\prime}|\mathrm{x},\theta,\mathrm{u})$ we approximate the intractable log-likelihood $g(\cdot)$ by employing a first order Taylor expansion of $g(\cdot)$ around the current value $\mathrm{x}$. We obtain, then, the proposal density

[TABLE]

where $Q_{\theta^{\prime}}=(C^{-1}_{\theta^{\prime}}+{2\over\delta}I_{d})^{-1}$. To reduce the costs of generating from the proposal in (12) as well as evaluating the resulting acceptance ratio we follow Titsias and Papaspiliopoulos (2018) and we define new auxiliary variables $\mathrm{z}$ by using the reparametrization

[TABLE]

The expanded target in (10) becomes now

[TABLE]

while the proposal in (12) can be written as

[TABLE]

where $\mathcal{Z}(\mathrm{z},\theta^{\prime})=\mathcal{N}_{d}(\mathrm{z}|0,C_{\theta^{\prime}}+(\delta/2)I_{d})$. Notice that, importantly, the proposed $\mathrm{x}^{\prime}$ becomes conditionally independent from the current $\mathrm{x}$ given the auxiliary variable $\mathrm{z}$. In particular, the proposal $q(\mathrm{x}^{\prime},\theta^{\prime}|\mathrm{x},\theta,\mathrm{u})$ in (11) becomes

[TABLE]

and the acceptance ratio of the move is

[TABLE]

where $\nu(\mathrm{z},\mathrm{x})=\left(\mathrm{z}-\mathrm{x}-(\delta/4)\nabla g(\mathrm{x})\right)^{\top}\nabla g(\mathrm{x})$. In our application we choose $q(\theta^{\prime}|\theta)$ to be a simple random walk proposal distribution with step-size $\kappa$. Algorithm 3 summarizes the proposed Metropolis-Hastings step for the joint update of $\mathrm{x}$ and $\theta$. Finally, we note that following Titsias and Papaspiliopoulos (2018) we learn $\kappa$ and $\delta$ during the burn-in period by decomposing step $2$ in Algorithm 3 in two sub-steps. We first update $\mathrm{x}$ alone and tune $\delta$ in order to achieve an acceptance ratio between $50\%$ and $60\%$. Then, we update $(\mathrm{x},\theta)$ jointly and we tune $\kappa$ to obtain acceptance ratio in the range of $20\%$ to $30\%$.

4.2.2 Interweaving strategy for the parameters

Although drawing the latent path $\mathrm{x}$ jointly with the parameters $\theta$ aims to overcome convergence problems that occurring due to their high dependence, it is well-known that the parametrization of latent Gaussian models as well as of more general Bayesian hierarchical models plays a key role in the efficiency of MCMC algorithms (Papaspiliopoulos et al., 2007). To address convergence issues that arise from the parametrizations of the latent log-volatility models and from the latent factor model for the jump intensities we employ the, so-called, ancillarity–sufficiency interweaving strategy (ASIS) developed by Yu and Meng (2011); see also Kastner and Frühwirth-Schnatter (2014) for the benefits of the ASIS in the case of univariate SV models and Simpson et al. (2017) where the ASIS has been applied in more general dynamic linear models. In the case of a latent Gaussian model with latent states $\mathrm{x}$ and parameters $\theta$ the application of ASIS is briefly described as follows. After the joint update of $\mathrm{x}$ and $\theta$ with the Metropolis-Hastings step described by Algorithm 3 we transform the latent path $\mathrm{x}$ to its non-centred parametrization $\tilde{\mathrm{x}}$ by using an one-to-one transformation. We re-draw then the parameters $\theta$ from $p(\theta|\tilde{\mathrm{x}})$ before transforming $\tilde{\mathrm{x}}$ back to $\mathrm{x}$ by using the inverse transformation.

To employ ASIS in our application we first note that equations (1)-(2) and (5)-(8) define the centred parametrizations (Papaspiliopoulos et al., 2007) of the corresponding models. The non-centred parametrization for the latent states of a latent Gaussian model can be found by following the reasoning developed by Yu and Meng (2011): under a non-centred parametrization $\mathrm{x}$ has to be an ancillary statistic for $\theta$, i.e., the distribution of $\mathrm{x}$ needs to be free of $\theta$. We note that there are cases where a many-to-one transformation is needed (Papaspiliopoulos et al., 2003) but this is not required in our application. It is easy to see that in the case of the latent log-volatilities $h_{it}$ defined in (2) the variables $\tilde{h}_{it}=(h_{it}-\mu_{i})/\sigma_{i\eta}$ are in a non-centred parametrization while in the case of the latent factors defined by (7) the transformation $\Gamma_{t}=F_{t}-AF_{t-1}$ results in a non-centred factor model for the jump intensities. ASIS in completed by re-drawing the log-volatility parameters $\mu_{i},\phi_{i}$ and $\sigma^{2}_{i\eta}$ and the factor parameters $A$ under the non-centred parametrizations and by using the inverse transformation to return back to $h_{it}$ and $F_{t}$.

In the supplementary material we provide all the details for our implementation of ASIS. To emphasize the gains obtained we conducted simulation studies which are also presented in the supplementary material. By describing briefly the results we note that combining ASIS with the auxiliary gradient-based sampler (Algorithm 3) improves substantially the efficiency of the proposed MCMC algorithm (Algorithm 1) in drawing samples from the parameters $\mu_{i}$ and $\sigma^{2}_{i\eta}$ of the log-volatility processes (see Figures S.$2$ and S.$3$ in the supplementary material) as well as from the parameters in the matrix $A$ of the latent factor model for the jump intensities; see Figure S.$4$ in the supplementary material.

4.3 Approximation by using a mixture of Gaussian distributions

For a given stock, sampling the whole path of the latent log-volatilities $H_{i}$ at each MCMC iteration could be performed by utilizing the methodology developed by Chib et al. (2002) and Nakajima and Omori (2009). Their methods are based on the idea of Kim et al. (1998) to approximate a SV without jumps model by using a mixture of Gaussian distributions.

The advantage of their approach is that the update of $H_{i}$ can be conducted by drawing samples directly from its full conditional distribution which is a $(T+1)$-variate Gaussian with tridiagonal inverse covariance matrix. The jumps sizes and their parameters should be updated after integrating out the mixture component indicators from the likelihood of the approximated model, otherwise convergence of the corresponding MCMC algorithm will be slow (Nakajima and Omori, 2009). Integrating out the mixture component indicators results in a partially collapsed Metropolis-Hastings within Gibbs sampler which has to be implemented with care since permuting the order of the updates can change the stationary distribution of the chain (Van Dyk and Jiao, 2015). If, for example, the update of the mixture indicators is conducted in the step which is intermediate between the update of the latent log-volatilities and the update of their parameters, then the Markov chain will not converge to the desired stationary distribution.

More importantly, by using the methods developed in Chib et al. (2002) and Nakajima and Omori (2009), samples from the posterior of the parameters of the exact model are obtained with an importance sampling step. In this step, weights that are equal to the ratio between the likelihood of the exact and the approximated model, are assigned to the MCMC samples obtained from the approximated model. Although in the case of SV models without jumps the described importance sampling step works without problems, it is known (Johannes et al., 2009) that in the case of models with jumps such weights may have large variance since a few or only one of them could be much larger than the others. Thus, the output of the importance sampling step will not be suitable for any Monte Carlo calculation. For an illustration of the problem we conducted a simulation experiment in which we calculated the required importance sampling weights for SV models with and without jumps. In particular, we simulated $p=1$ time series with $T=1,500$ log-returns from each model by utilizing equations (1)–(2) and omitting the jump component in (1) to simulate from the model without jumps. In the case of the SV model with jumps we assumed independent jump intensities over time following the modelling approach presented in Section 3.2.1. We chose values for the parameters of the log-volatility process which are common in the financial literature (Kim et al., 1998); $\mu_{1}=-0.85,\phi_{1}=0.98$ and $\sigma_{1\eta}=0.15$ and we simulated the sizes of possible jumps from the Gaussian distribution with $\mu_{1\xi}=0$ and $\sigma_{1\xi}=3.5$ (Eraker et al., 2003). Then, we drew $3,000$ (thinned) posterior samples of the parameters and latent states of the approximated models by using the corresponding MCMC algorithms.

Figure 3 displays the log-normalised importance sampling weights plus the logarithm of their number. The variance of the normalized importance sampling weights should become smaller as the importance sampling approximation improves. This implies that the histograms in Figure 3 should be centered at zero in the case of an accurate enough approximation; see the corresponding plots in Kim et al. (1998) and Omori et al. (2007). The left panel of Figure 3 indicates that there is little effect of the mixture approximation in the case of the model without jumps, but the right panel points out that this is not true in the case of the SV with jumps model. The bad performance of the approximation strategy in models with jumps is due to the fact that the mixture model has been proposed by Kim et al. (1998) and improved by Omori et al. (2007) for SV models without jumps. By adding a jump component in the model the mixture approximation has to take into account the uncertainty in the estimation of the jump process and this is not the case in the approximation used by Chib et al. (2002) and Nakajima and Omori (2009).

5 Out-of-sample forecasting

To compare the predictive performance of different SV models we develop a computational framework for out-of-sample forecasting. Based on $T$ in-sample and $\ell$ out-of-sample observations we assess the forecasts obtained from the various models by calculating proper scoring rules (Gneiting and Raftery, 2007). A scoring rule assigns a numerical score to a probabilistic forecast by taking into account the posterior predictive distribution and the realized observation and it is proper if the expected score is maximized when the estimated predictive distribution coincides with the true distribution. Proper scoring rules allow for the joint evaluation of the calibration (compatibility of forecasts and observations) and the sharpness (concentration of the predictive distributions); see, for example, Geweke and Amisano (2010) and Krüger et al. (2020) for more detailed discussions. We focus on three proper scoring rules: (i) logarithmic scores, and particularly, the predictive log Bayes factors which compare the models with respect to their marginal likelihoods, (ii) continuous ranked probability scores which are employed for the assessment of the predictive cumulative distribution functions (CDFs) and (iii) interval scores which are utilized for the evaluation of the prediction credible intervals. We complete a full-fledged assessment of the predictive performance of the different models by calculating root mean squared errors of the obtained forecasts.

To calculate these summary measures we draw samples from the posterior predictive distributions of interest and we estimate their densities by utilizing sequential Monte Carlo methods (Chopin and Papaspiliopoulos, 2020b). In the rest of the Section we present the computational framework that we built in order to conduct out-of-sample forecasts with univariate SV models with and without jumps as well as by using the proposed joint modelling approach for the jump intensities. To evaluate the accuracy of the obtained estimators we calculate the effective sample size of the samples drawn from the posterior predictive distributions. The effective sample size of a given sample is commonly employed (Neal, 2001; Chopin et al., 2013) to assess the variance of estimators obtained from weighted rather than independent samples; see in the supplementary material for more details. For the remaining of this section we assume that the parameters $\Theta$, $A$ and $W$ are fixed to their in-sample estimated posterior means and we omit them from our notation.

5.1 Prediction with independent SV models

Let $\mathbb{M}_{\mathrm{svj}}$ and $\mathbb{M}_{\mathrm{sv}}$ be the model that consists of $p$ independent, univariate SV models with and without jumps respectively. Let also $R_{1:T}$ and $R_{T+1:T+\ell}$ denote in-sample and out-of-sample observations. To calculate predictive Bayes factors and the continuous ranked probability and interval scores we consider the posterior predictive distributions with densities $p(r_{it}|r_{i,1:t-1})$, $t=T+1,\ldots,T+\ell$. The predictive log Bayes factors are defined as

[TABLE]

where for both $\mathbb{M}_{\mathrm{svj}}$ and $\mathbb{M}_{\mathrm{sv}}$ we have that

[TABLE]

Moreover, if $l_{it}$ and $u_{it}$ denote the $\alpha/2$ and $1-\alpha/2$ quantiles of a $(1-\alpha)\times 100\%$ prediction interval for the value of $r_{it}$ then its interval score is given by the formula

[TABLE]

We note, finally, that the calculation of the continuous ranked probability scores that we employ for the evaluation of predictive CDFs is also based on samples from $p(r_{it}|r_{i,1:t-1})$; see the supplementary material for details.

We draw the desired samples and we estimate the densities $p(r_{it}|r_{i,1:t-1})$ based on the formula

[TABLE]

and by utilizing sequential Monte Carlo methods also known as particle filters (Chopin and Papaspiliopoulos, 2020a). From the output of a particle filter algorithm that targets the posterior in (17) we are able to draw samples from the posterior predictive distribution of $r_{it}$ and to estimate the normalizing constants $p(r_{it}|r_{i,1:t-1})$. See in the supplementary material for more details and for the pseudocode of the particle filter algorithm.

5.2 Prediction by using joint modelled jump intensities

Let $\mathbb{M}_{\mathrm{mvj}}$ be our proposed SV with jumps model, defined by equations (1)–(3) and (5)–(8). To compare the predictive performance of $\mathbb{M}_{\mathrm{mvj}}$ with the performance of competing models we need to estimate the predictive densities $p(R_{T+j}|R_{1:T+j-1},\mathbb{M}_{\mathrm{mvj}})$, for each $j=1,\ldots,\ell$. This could be achieved with a particle filter algorithm, constructed similarly to the case of independent SV models, as described in section 5.1. However, in the case of $\mathbb{M}_{\mathrm{mvj}}$, particle filter algorithms deliver poor estimations of the marginal likelihoods since the importance weights, calculated at each of their steps, have large variance.

We develop an alternative sequential importance sampling algorithm to obtain the desired estimates which is based on the annealed importance sampling (AIS) proposed by Neal (2001). AIS utilizes the method of simulated annealing (Kirkpatrick et al., 1983) to construct an importance sampling algorithm for drawing weighted samples from an unnormalized target distribution. The proposal distribution of the algorithm is based on a sequence of auxiliary distributions and Markov transition kernels that leave them invariant. As in any importance sampling algorithm, the sample mean of the importance weights is an estimation of the ratio between the normalizing constants of the target and proposal distributions. Moreover, the output of AIS can be used for the estimation of expectations with respect to any of the auxiliary distributions used for the construction of the algorithm. Here we exploit these features of AIS to draw samples from the posterior predictive distributions of interest and to estimate their densities. In the rest of the section we provide a detailed description of how we apply AIS, $\mathbb{M}_{\mathrm{mvj}}$ is suppressed throughout except when it is necessary.

Let $g_{0},\ldots,g_{\ell}$ be the required sequence of auxiliary distributions. By noting that each $g_{j}$ must be proportional to a computable function and satisfies $g_{j}\neq 0$ wherever $g_{j-1}\neq 0$ we set $g_{\ell}(H_{0:T+\ell},F_{0:T+\ell})=p(H_{0:T+\ell},F_{0:T+\ell}|R_{1:T+\ell})$ and

[TABLE]

for $j=0,\ldots,\ell-1$. Samples from $g_{j}$, $j\geq 1$, are obtained by drawing samples from $g_{j-1}$ using a few iterations of Algorithm 1 and the transition densities $p(H_{T+j}|H_{T+j-1})$ and $p(F_{T+j}|F_{T+j-1})$. To sample from $g_{0}$ we utilize samples from $p(H_{0:T},F_{0:T}|R_{1:T})$. We compute the importance weights of the obtained samples by first noting that $g_{j}$ are proportional to the computable functions

[TABLE]

and $g^{*}_{\ell}(H_{0:T+\ell},F_{0:T+\ell})=p(H_{0:T+\ell},F_{0:T+\ell},R_{1:T+\ell})$. Therefore, the importance weight for every sampled point is

[TABLE]

and $\Omega^{s}_{T+1}=p(R_{T+1}|H_{T+1}^{s},F_{T+1}^{s})$, $s=1,\ldots,S$, where $S$ is the number of samples obtained by AIS and

[TABLE]

We also note that the infinite sum in (5.2) can be truncated without affecting the results of the algorithm (Johannes et al., 2009). We refer to the supplementary material for the pseudocode that we used to construct the described AIS algorithm.

The sample means of the weights calculated at each step of AIS provide estimators of the marginal likelihoods $p(R_{T+1:T+j}|R_{1:T})$. However, the resulting weights have large variance and the corresponding estimates will be inaccurate. Nevertheless, and in contrast with a particle filter that targets $p(H_{0:T+\ell},F_{0:T+\ell}|R_{1:T+\ell})$, we can use the samples $\{H^{s}_{0:T+j},F^{s}_{0:T+j}\}_{s=1}^{S}$ to estimate the quantities $p(r_{i,T+1:T+j}|R_{1:T})$ for each stock separately. We propose the following sampling procedure. First, note that

[TABLE]

This implies that by using only the $i$th element of the product in (5.2) we obtain an estimator of the $i$th marginal of $p(R_{T+1:T+j}|R_{1:T})$ as follows. We set $\hat{p}(r_{i,T+1:T+j}|R_{1:T})=\sum_{s=1}^{S}\omega^{s}_{i,T+j}/S$, where

[TABLE]

and $\omega^{s}_{i,T+1}=1$. Thus, $p(r_{i,T+1:T+j}|R_{1:T},\mathbb{M}_{\mathrm{mvj}})$ is estimated by $\hat{p}(r_{i,T+1:T+j}|R_{1:T})$ accurately since the latter is based on importance sampling weights with reduced variance and can be used for the estimation of log Bayes factors of the form

[TABLE]

where $p(r_{i,T+1:T+j}|R_{1:T},\mathbb{M}_{\mathrm{svj}})$ can be estimated with the methods presented in Section 5.1. We also note that by assuming that the joint predictive density of all returns is estimated by the product of their marginal predictive densities, as in the independent across stocks SV with jumps models, we estimate the approximate log Bayes factors

[TABLE]

Finally, $\{h^{s}_{i,0:T+j},F^{s}_{0:T+j}\}_{s=1}^{S}$ can be utilized in order to draw samples from the posterior predictive distributions needed to calculate the continuous ranked probability and the interval scores; see the supplementary material for details.

6 Real data application

Here we present results from the application of the developed methodology on real data. In the supplementary material we present results from simulation studies used to test our methodology. Further specifics of the data analysis, such as parameter estimates for the log-volatility processes, can be also found in the supplementary material as well as in Alexopoulos (2017).

6.1 Data and MCMC details

We applied the developed methods on time series consisted of daily log-returns from stocks of the STOXX Europe $600$ Index downloaded from Bloomberg between $10/1/2007$ to $11/6/2014$. We removed $29$ stocks by requiring each stock to have at least $1000$ traded days and no more than $10$ consecutive days with unchanged price. The final dataset consisted of $p=571$ stocks and $1947$ traded days. We used daily log-returns observed in the first $T=1917$ days as in-sample observations to estimate all models described through our article and the last $\ell=30$ observations to test their predictive ability. In our dynamic factor model we used $K=3$ factors. We ran all the MCMC algorithms for $250,000$ iterations using a burn-in period of $100,000$ sampled points and we collected $1,000$ (thinned) posterior samples. In the supplementary material we give the computing times for the real and simulated data analysis that we have implemented and we evaluate the accuracy of the sequential Monte Carlo methods that we used to conduct the assessment of the predictive performance of the different models.

6.2 Results

6.2.1 Bayesian inference for parameters and latent states

Figure 4 depicts the posterior means of the factor paths along with the posterior distributions of their persistent parameters $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$. The posterior densities of $\alpha_{1}$ and $\alpha_{2}$ indicate strong evidence of autocorrelation in the first two factor processes while that of $\alpha_{3}$ suggests that a model with only $K=2$ factors adequately captures the dependence of the jumps over times and across stocks. Therefore, all results in the following are based on $K=2$ factors.

By noting that we have applied our methods on a very high-dimensional dataset we present posterior distributions for summary statistics related to the latent volatility and jump processes. Figure 5 shows, for each stock, the $95\%$ credible interval of the temporal average of the volatility process. From the visual inspection of the Figure we note that the stocks can be separated in those with low, medium and high volatilities. More precisely, the posterior temporal average of the volatilities for the majority of the stocks is estimated lower than $5$, for a few number of stocks the temporal average ranges between $5$ and $15$ while for less than $10$ stocks the posterior temporal average has $95\%$ credible interval with lower limit greater than $15$. See also the supplementary material for additional results regarding Bayesian inference about the log-volatility processes and their parameters. Figure 6 presents the $95\%$ credible intervals for the posterior distribution of the temporal sum of jumps for each stock with respect to the economic sector that it belongs to. The Figure indicates that for each stock have been identified $100$ jumps on average between $10/1/2007$ and $11/6/2014$. We also note that the number of stocks for which almost no jumps have been estimated is very small and such stocks mainly belong to the financial, the industrial and the consumer- cyclical and non-cyclical, the energy and the communications sectors. In the supplementary material we present more results related to the Bayesian analysis of the jump processes of the SV models while Figure 1 offers a visual illustration of the estimated jumps underlying the time evolution of the stock returns.

As already noted in the introduction of the paper the top panel of Figure 1, which has been created by using an empirical method for jump identification, visualizes a dense estimation of the underlying jump process. The middle panel which corresponds to independent modelling of jumps depicts a very sparse jump process whereas the bottom panel indicates less number of jumps than the empirical method but more than the independent SV models. To evaluate the estimation of the jump processes in the three panels we first emphasize that the true jumps in the stock prices are unobserved events. Therefore, drawing conclusions for the accuracy of the jumps estimation by only visually inspecting the Figure is difficult. We note, however, that the dense jump structure obtained from the empirical method can be explained by the fact that the time evolution of the volatility of the stock returns is not taken into account; resulting in a possible overestimation of the number of jumps. The latter combined with the more narrow credible intervals for the log-volatilities delivered by SV models with jumps (see Figure S$.12$ in the supplementary material) indicates that an efficient jump identification results in accurate estimation of the log-volatility process without increasing unnecessarily its level due to the presence of large price movements and, thus, more precise forecasts of future returns are expected. See also the following Section where we compare the predictive performance of SV models with and without jumps. To assess the jump structures estimated by employing independent and joint models for the jump intensities of SV models we rely on their forecasting performance. In particular, as also described below, the proposed joint modelling approach for the jump intensities results in a more accurate prediction of future returns. We conclude, therefore, that the estimation of the underlying jumps process depicted by the third panel of Figure 1 is closer to the real jump activity of the stock prices allowing further improvement in the estimation of the log-volatilities and more accurate prediction of future stock returns.

6.2.2 Predictive performance

Here we evaluate the predictive performance of the proposed SV model with jointly modelled jump intensities as well as the corresponding performance of univariate SV models with and without jumps which are also specified by equations (1)-(3) by making suitable assumptions. In the supplementary material we compare the forecasting performance of the proposed model with alternative SV models which account for the heavy tails of the distribution of the observed returns.

The left panel of Figure 7 compares the predictive performance of independent, univariate SV models with and without jumps by presenting predictive log Bayes factors. The Figure is constructed by considering the log Bayes factors defined in (15) for the $\ell=30$ out-of-sample observations of our dataset. Thus, positive log Bayes factor indicates that the predictions obtained with univariate SV models with jumps are more accurate than those of SV models without jumps. According to Kass and Raftery (1995) a value of the log Bayes factor greater than $5$ provides strong evidence in favour of the model with jumps. The increasing nature of the log Bayes factor reassures that as more data are collected and as more jumps are identified in stock returns, the evidence that the model with jumps outperforms the model without jumps increases. In the right panel of Figure 7 we compare the predictive performance of our proposed SV with jumps model in which the jump intensities are modelled by using a dynamic factor model with SV models in which the jump intensities are independent over time and across stocks. The Figure compare the two models by presenting the logarithms of the approximate log Bayes factors in (20). Clearly, there is strong evidence that our proposed model outperforms the independent modelling of stock returns with a SV with jumps formulation.

Figure 8 compares the forecasts of the proposed model with those obtained from SV models without jumps and with independent jumps respectively by employing their probabilistic evaluation (Gneiting and Raftery, 2007). Joint modelling of the intensities leads to $95\%$ prediction credible intervals with the lowest interval scores at each one of the out-of-sample points for the vast majority of the stocks. The definition of the interval score in (16) rewards the forecaster that delivers narrow prediction intervals and when the observation is outside the predicted interval a penalty proportional to the significance level of the interval is incurred. It is thus indicated that the proposed SV model delivers the most accurate prediction intervals compared to those obtained from univariate SV models with and without jumps. The fact that Figure 8 displays the median instead of the mean of the $571$ interval scores emphasizes an important characteristic of our modelling approach. The efficient identification of the in-sample jumps has resulted in estimating lower volatilities for the stock returns; see Figures S$.11$ and S$.12$ in the supplementary material. This is depicted in the conducted out-of-sample exercise by more narrow prediction intervals for a given coverage probability. Moreover, since the prediction of future jumps is based on common across the stocks, autoregressive factors, we expect that jumps that occur suddenly in a small number of stocks at a given day are not easily predictable. On the other hand, the large volatility intervals delivered from SV models without jumps or models with less accurate jump identification can accidentally predict more efficiently these type of movements. The median is chosen to avoid taking into account such predictions. The superiority of the proposed approach is also emphasized by Figure 8 by the lower continuous ranked probability scores that correspond to the proposed SV model than those of the univariate SV models. This implies that the posterior predictive distribution of the proposed model are more concentrated than those obtained from univariate SV models with and without jumps while deliver predictions that are not easily distinguishable from the materialized observations. This is a desired property of probabilistic forecasts described as the main goal of probabilistic forecasting; see for example Gneiting and Katzfuss (2014). Finally, Figure S$.14$ in the supplementary material which presents root mean squared errors for the predictions obtained from the different models indicates that joint modelling of the jump intensities improves slightly the point forecasts. We conclude that the proposed modelling approach offers a more accurate identification of jumps and consequently a more precise estimation of the returns volatility compared to the one conducted by using SV models without or with independent jumps. As a result, the predictive performance of the proposed model outperforms the other models. This makes the suggested model a competitive alternative for practitioners who are interested in obtaining accurate predictions for the stocks of large portfolios.

7 Discussion

We have developed a general modelling framework together with carefully designed MCMC algorithms to perform Bayesian inference for SV models with Poisson-driven jumps. We have shown that for the data we applied our models there is evidence that, with respect to predictive Bayes factors, (i) univariate SV models with jumps outperform univariate SV models without jumps and (ii) models that jointly model jump intensities with a dynamic factor model outperform SV models with jumps that are applied independently across stocks. We feel that (ii) is an interesting result that adds considerable insight to the growing literature of modelling financial returns with jumps, adding to the observation by Aït-Sahalia et al. (2015) that there is indeed predictability in jump intensities.

There are various issues that have not been addressed in this paper. The proposed MCMC algorithms can be easily extended in order to conduct inference for the more realistic (Harvey and Shephard, 1996) class of SV models which allow for correlation between the error terms $\epsilon_{it}$ and $\eta_{it}$ in (1) and (2), a phenomenon which is commonly known as leverage effect. In this case samples from their full conditional distributions can be obtained by using the auxiliary gradient-based sampler (Titsias and Papaspiliopoulos, 2018) to draw each correlation parameter jointly with $\phi_{i}$, $\sigma_{i\eta}^{2}$ and $H_{i}$. Other modelling aspects include the relaxing of independence of $E_{t}$ in (8) or the assumption that $A$ is lower diagonal as in the models proposed in Dellaportas and Pourahmadi (2012).

Finally, a series of interesting financial questions can be addressed with our models by exploiting the fact that panel of stock returns carry additional, possibly important, information. For example, in our dataset one can explore the effects of country and sector effects or could investigate whether the joint tail risk dependence does or does not change before, during and after the financial crisis in Europe. We feel that our proposed models together with our algorithmic guidelines will serve as useful tools for such future research pursuits in financial literature.

Empirical detection of jumps

Let $R_{i}=(r_{i1},\ldots,r_{iT})$ be the vector with log-returns of the $i$th stock. The detection of jumps is based on the statistic

[TABLE]

where $MAD=1.48\times median(|r_{it}-median(R_{i})|)$ is the Mean Absolute Deviation from the median and it is a robust estimator for the scale parameter of a Gaussian distribution (Rousseeuw and Croux, 1993). A log-return is considered a jump if for this log-return the statistic in (21) has value greater than $3$; see Rousseeuw and Croux (1993) for a detailed discussion about the use of the MAD estimator in outlier detection problems.

Posterior jump probabilities

Figure $1$ in the main paper was constructed by considering a price movement as jump if the jump probability for the given stock at that particular day was estimated greater than $0.5$. Figure 9 indicates that the majority of the posterior jump probabilities have been estimated to be close to zero. In particular, $98\%$ of them are lower than $0.5$ and from those $95\%$ is less than $0.2$.

Choosing hyperparameters for the factor model

We have assumed that $b_{i}\overset{iid}{\sim}\mathcal{N}(\mu_{b},\sigma^{2}_{b})$, $w_{ik}\overset{iid}{\sim}\mathcal{N}(0,\sigma_{w}^{2})$ and $\alpha_{k}\overset{iid}{\sim}U(-1,1)$, $i=1,\ldots,p$ and $k=1,\ldots,K$. Our aim is to specify the hyperparameters $\mu_{b}$, $\sigma_{b}^{2}$ and $\sigma_{w}^{2}$ to be such that the resulting prior on $\lambda_{it}$ has mode, mean and variance close to the corresponding values of the $Gam(1,50)$ distribution which was used as prior in the case of modelling the jump intensities independently over time and across stocks.

From equations $(3.5)$ and $(3.6)$ of the main paper we have that

[TABLE]

where we have set, see Section $3.3$, $\lambda^{*}=0.15$. To specify $\mu_{y}$ and $\sigma_{y}^{2}$ we note that equations $(3.6)$–$(3.8)$ in the main paper imply

[TABLE]

Thus, by assuming that $\alpha_{k}=\mathrm{E}(\alpha_{k})=0$ we have

[TABLE]

From (22) we have that the density of $\lambda_{it}$ is

[TABLE]

Therefore, the location of its mode satisfies the equation

[TABLE]

which can be solved, for given values of $\lambda^{*}$, $\mu_{y}$ and $\sigma^{2}_{y}$, with numerical methods. Here we use the R-package rootSolve (rootSolve) to find all the roots of the equation in (25).

Since our aim is to match $\mu_{y}$, $\sigma^{2}_{y}$ and the mode of the density in (24) with the corresponding quantities implied from the $Gam(1,50)$ distribution we first note that if $\lambda_{it}\sim Gam(1,50)$, right truncated at $\lambda^{*}$, then for $\lambda^{*}=0.15$ we have that

[TABLE]

and that the mode of the $Gam(1,50)$ distribution is located at zero. By noting that in our applications we have used models with $K=2$ and $K=3$ latent factors and taking into account (23) we set $\sigma_{w}^{2}=0.5$ and $\sigma^{2}_{b}=1$. Then, from (23) we have that $\sigma^{2}_{y}$ is close to the variance that we obtain for $y_{it}$ in the univariate case; see equation (26). To specify $\mu_{b}$ except from matching the means in (23) and (26) we are also trying to set the mode of the density in (24) close to zero. Table 2 displays the mean, the variance and the mode of (24) for three different values of $\mu_{b}$. Note that by choosing $\mu_{b}=-2.4$ we have that the mean and the variance of the density in (24) are quite close to the corresponding moments of the $Gam(1,50)$ distribution while by taking smaller values for $\mu_{b}$ the mode in (24) becomes closer to zero.

Proof of Proposition 1

Proof.

A discrete distribution with pmf $p(n)$ is called log-concave when the inequality $p(n)^{2}\geq p(n-1)p(n+1)$ is satisfied for all $n$. Let $\sigma^{2}_{it}=\exp(h_{it})$, in the case of joint modelled jump intensities we have that \lambda_{it}=\lambda^{*}/\big{(}1+\exp(-b_{i}-W^{\prime}_{i}F_{t})\big{)} and in the case of independent modelling $\lambda_{it}$ is a constant. We denote with $p(n)$ the pmf $p(n_{it}=n|h_{it},\mu_{i\xi},\sigma^{2}_{i\xi},b_{i},W_{i},F_{t},r_{it})$. By integrating out the jump sizes $\Xi_{i}$ from the likelihood we have that

[TABLE]

To prove log-concavity for $n=1,2,...$, we need to show that for $n=2,3,...$

[TABLE]

or equivalently that,

[TABLE]

or that

[TABLE]

where $L=\frac{2\left(\sigma^{2}_{i\xi}(r_{it}-n\mu_{i\xi})^{2}+\mu_{i\xi}(\sigma^{2}_{it}+n\sigma^{2}_{i\xi})\right)^{2}}{(\sigma^{2}_{it}+n\sigma^{2}_{i\xi})\left((\sigma^{2}_{it}+n\sigma^{2}_{i\xi})^{2}-\sigma^{4}_{i\xi}\right)}$. The last inequality holds; $L>0$ and it is readily shown that the numerator is greater than the denominator for all $n=2,3...$. ∎

Proof of Proposition 2

Proof.

Let $\boldsymbol{1}_{n}^{\prime}$ be an $n$-dimensional column vector with ones and $I_{n}$ the $n\times n$ identity matrix. For each $i=1,\ldots,p$ and $t=1,\ldots,T$ we have that $\xi_{it}^{\kappa}\overset{iid}{\sim}N(\mu_{i\xi},\sigma^{2}_{i\xi})$, $\kappa=1,\ldots,n_{it}$. Let $\sigma^{2}_{it}=\exp(h_{it})$ then from the canonical form of the Gaussian distribution we have that

[TABLE]

where $c=(\frac{\mu_{i\xi}}{\sigma^{2}_{i\xi}}+\frac{r_{it}}{\sigma_{it}^{2}})\boldsymbol{1}_{n_{it}}$ and $Q_{i\xi}=Q_{1}+Q_{2}$, with $Q_{1}=\sigma^{-2}_{i\xi}I_{n_{it}}$ and $Q_{2}=\sigma_{it}^{-2}\boldsymbol{1}_{n_{it}}\boldsymbol{1}_{n_{it}}^{\prime}$. From miller1981inverse we have that

[TABLE]

where $\nu=tr(Q_{2}Q_{1}^{-1})$ or $(1+\nu)^{-1}=\dfrac{\sigma^{2}_{it}}{\sigma^{2}_{it}+n_{it}\sigma^{2}_{i\xi}}$. Therefore

[TABLE]

∎

Auxiliary gradient-based sampler

Here we describe how we apply the auxiliary gradient-based sampler (Titsias and Papaspiliopoulos, 2018) to draw, at each MCMC iteration, the whole latent paths of the log-volatilities $H_{i}$ jointly with the parameters $\phi_{i}$ and $\sigma_{i\eta}^{2}$ and the whole path of the latent factors $F$ jointly with the parameters $A$.

Sampling latent log-volatilities and parameters

To sample, for a fixed $i$, $H_{i}$ jointly with $\phi_{i}$ and $\sigma_{i\eta}^{2}$, we need to draw samples from the distribution with density

[TABLE]

where the the jump sizes $\Xi_{i}$ have been integrated out, $\pi(H_{i}|\phi_{i},\sigma_{i\eta}^{2},\mu_{i})$ denotes the prior of the latent autoregressive log-volatility process defined in $(3.2)$ and $\pi(\phi_{i})$ and $\pi(\sigma_{i\eta}^{2})$ denote the densities of the priors for the parameters $\phi_{i}$ and $\sigma^{2}_{i\eta}$ as described in Section $3.2$. To draw samples from (27) we employ Algorithm $3$ presented in Section $4.2$ of the main paper where $\mathrm{x}$ is the latent volatility vector $H_{i}$ and $\theta$ is consisted of the parameters $\phi_{i}$ and $\sigma^{2}_{i\eta}$; it is also easy to see that $C_{\theta}^{-1}$ has a tridiagonal form and both sampling from the proposal and calculating the acceptance ratio in $(3.6)$ can be conducted, efficiently, based on the Cholesky decomposition of tridiagonal matrices. Finally, we note that although Algorithm $3$ corresponds to a latent Gaussian model in which the Gaussian field has zero mean, each log-volatility vector $H_{i}$ has mean $\mu_{i}$; we take this into account and we present below how the steps of Algorithm $3$ of the main paper are adjusted in order to draw the desires samples $(H_{i},\phi_{i},\sigma^{2}_{i\eta})$. In particular, we iterate through the following steps:

1. (i)

Sample Z_{i}\sim\mathcal{N}_{T+1}\big{(}Z_{i}|H_{i}+(\delta/2)\nabla g(H_{i}|R_{i},N_{i},\mu_{i\xi},\sigma^{2}_{i\xi}),(\delta/2)I_{T+1}\big{)} 2. (ii)

Sample from $p(H_{i},\phi_{i},\sigma_{i\eta}^{2}|Z_{i},\mu_{i},\mu_{i\xi},\sigma^{2}_{i\xi},N_{i},R_{i})$ by using a Metropolis-Hastings in which we transform \hat{\phi}_{i}=\log\big((\phi_{i}+1)/(1-\phi_{i})\big{missing}) and $\hat{\sigma}_{i\eta}^{2}=2\log\sigma_{i\eta}$ and we propose new values $(H_{i}^{*},\hat{\phi}_{i}^{*},\hat{\sigma}_{i\eta}^{2*})$ by first drawing $(\hat{\phi}_{i}^{*},\hat{\sigma}_{i\eta}^{2*})$ from the bivariate Gaussian density centred at $(\hat{\phi}_{i},\hat{\sigma}_{i\eta}^{2})$ with covariance matrix $\kappa_{i}I_{2}$ and then draw $H_{i}^{*}$ from

[TABLE]

with $Q=(C_{\theta}^{*-1}+\frac{2}{\delta}I_{T+1})^{-1}$, M^{*}_{i}=\dfrac{\mu_{i}(1-\phi^{*}_{i})}{\sigma_{i\eta}^{2*}}\big{(}1,\,\,1-\phi^{*}_{i},\ldots,1-\phi^{*}_{i},\,\,1\big{)} and

[TABLE]

where $\phi^{*}_{i}=(e^{\hat{\phi}_{i}^{*}}-1)/(e^{\hat{\phi}_{i}^{*}}+1)$ and $\sigma_{i\eta}^{2*}=e^{\hat{\sigma}_{i\eta}^{2*}}$. The values $(H_{i}^{*},\hat{\phi}_{i}^{*},\hat{\sigma}_{i\eta}^{2*})$ are accepted or rejected according to the ratio

[TABLE]

where \nu\big{(}Z_{i},H_{i}\big{)}=\big{(}Z_{i}-H_{i}-\frac{\delta}{4}\nabla g\big{(}H_{i}\big{|}R_{i},N_{i},\mu_{i\xi},\sigma^{2}_{i\xi}\big{)}\big{)}^{\prime}\nabla g\big{(}H_{i}\big{|}N_{i},\mu_{i\xi},\sigma^{2}_{i\xi}\big{)} and $J(\cdot,\cdot)$ denotes the Jacobian of the inverse of the transformation $(\phi_{i},\sigma_{i\eta}^{2})\rightarrow(\hat{\phi}_{i},\hat{\sigma}_{i\eta}^{2})$.

After drawing $(H_{i},\phi_{i},\sigma_{i\eta}^{2})$ we sample $\mu_{i}$ from its full conditional $\mathcal{N}(m_{i},s_{i}^{2})$, where

[TABLE]

and

[TABLE]

Ancillarity-sufficiency interweaving strategy for log-volatility parameters

In order to implement the ancillarity-sufficiency interweaving strategy (Yu and Meng, 2011) for the parameters $\mu_{i}$, $\phi_{i}$ and $\sigma_{i\eta}^{2}$ we first consider the model defined, in the main paper, by equations $(3.1)$ and $(3.2)$ with the non-centred parametrization of the log-volatilities. This is specified, for each $i=1,\ldots,p$ and $t=1,\ldots,T$, by

[TABLE]

[TABLE]

where $\epsilon_{it}$ and $\eta_{it}$ are independent random variables and $\tilde{h}_{it}=(h_{it}-\mu_{i})/\sigma_{i\eta}$. Then, we utilize random walk Metropolis-Hastings steps to draw $\mu_{i}$ and $\sigma_{i\eta}^{2}$ from their full conditionals obtained from

[TABLE]

and we sample $\phi_{i}$ from

[TABLE]

by using a Metropolis-Hastings step with proposal distribution

[TABLE]

Simulation experiments for samplers with and without interweaving

Figures 10 and 11 evaluate the efficiency of the auxiliary gradient-based sampler in drawing the parameters of the latent log-volatility path of univariate SV models. We compare three sampling schemes; the auxiliary gradient-based sampler with interweaving, without interweaving, and a scheme in which instead of interweaving we iterate twice the auxiliary sampler. To evaluate the different sampling schemes we compare the effective sample size (ESS) of the samples drawn from the posterior distributions of interest. The ESS of $S$ samples drawn by using an MCMC algorithm can be estimated as $v^{2}S/\nu_{0}$ where $v^{2}$ is the sample variance of the posterior samples and $\nu_{0}$ is an estimation of the spectral density of the Markov chain at zero. We used the R-package coda (rcoda) to estimate the ESS of the posterior samples obtained from the samplers under comparison.

The comparisons of the three sampling schemes are based on $10$ replicates of two datasets consisted of $T=1,500$ log-returns. These were simulated from the univariate SV model without jumps which is defined by $(3.1)$ and $(3.2)$ of the main paper, with $p=1$, and by omitting the jump component in $(3.1)$. For the remaining of this subsection we omit the subscript $i$ from our notation; e.g. $\mu\equiv\mu_{i}$ and $H\equiv H_{i}$. We set $\mu=-0.85$, $\sigma_{\eta}=0.15$ to generate both datasets and $\phi=0.98$ in the first and $\phi=0.90$ in the second. In each replicate of the datasets we applied the sampling schemes under comparison in order to draw samples for the latent log-volatility path $H$ and its parameters $\mu$, $\sigma_{\eta}$ and $\phi$. We update $H$ jointly with $\phi$ and $\sigma^{2}_{\eta}$ by using the auxiliary gradient-based sampler to draw from (27) (with $n_{t}=0$) and then we draw $\mu$ from its full conditional. In the case of combining the sampler with the ancillarity-sufficiency interweaving strategy we follow the procedure described in the previous section in order to update the parameters $\mu$, $\sigma_{\eta}$ and $\phi$ under the non-centered parametrization of the log-volatility process. We ran all the MCMC algorithms for $30,000$ iterations and we discarded the first $10,000$ as burn-in period obtaining, thus, $20,000$ posterior samples for the log-volatility path and its parameters for each replicate of the two simulated datasets. In the Figures 10 and 11 we present, for the datasets with $\phi=0.98$ and $\phi=0.90$ respectively, the boxplots with the ESS of the posterior samples of the parameters $\mu$, $\sigma_{\eta}$ and $\phi$ obtained from the $10$ replicates. The two figures indicate that by combining the auxiliary gradient-based sampler with the ancillarity-sufficiency interweaving strategy the ESS of $\sigma_{\eta}$ and $\mu$ is clearly improved, while the three schemes have almost the same efficiency in drawing $\phi$.

Sampling the latent factors and their parameters

To draw the latent factors $F$ jointly with the parameters $A$ we utilize the auxiliary gradient-based sampler summarized by Algorithm $3$ in the main paper and we target the distribution with density

[TABLE]

where $\pi(A)$ denotes the prior density of $A$ and $C_{A}$ is the covariance matrix (which depends on the values of the parameters $A$) of the autoregressive process defined by $(3.7)$ and $(3.8)$. Since the latent factor model defined by equations $(3.6)$–$(3.8)$ corresponds to a latent Gaussian model with zero mean Gaussian field, the application of Algorithm $3$ of the main paper is straightforward. The latent path $\mathrm{x}$ in Algorithm $3$ corresponds to the path of the latent factors $F$ and the parameters $\theta$ are the parameters $\alpha_{1},\ldots,\alpha_{K}$ of the matrix $A$. We note that the parameters $\alpha_{k}$, $k=1,\ldots,K$, are defined on the interval $(-1,1)$. Therefore, we apply the transformation \hat{\alpha}_{k}=\log\big((\alpha_{k}+1)/(1-\alpha_{k})\big{missing}) in order to obtain parameters in the real line. In the acceptance ratio in $(4.6)$ of the main paper we need thus to account for the Jacobian of the inverse transformation.

Using the interweaving strategy to improve sampling of $A$

We combine the auxiliary gradient-based sampler (Algorithm $3$ in the main paper) with the ancillarity-sufficiency interweaving strategy in order to reduce the autocorrelation of the parameters in $A$. More precisely, after sampling jointly the paths $F$ of the latent factors and the parameters $A$ we set $\Gamma_{0}=DF_{0}$, where $D$ is a $K\times K$ diagonal matrix with elements $(1-\alpha_{k}^{2})^{1/2}$, $k=1,\ldots,K$, and $\Gamma_{t}=F_{t}-AF_{t-1}$, for $t=1,\ldots,T$. Then we use a random walk Metropolis-Hastings step that targets the density

[TABLE]

with $\lambda_{it}=\lambda^{*}/(1+\exp(-b_{i}-W^{\prime}_{i}F_{t}))$. We note that $F_{t}$ and $\Gamma_{t}$ are related by $F_{0}=D^{-1}\Gamma_{0}$ and $F_{t}=AF_{t-1}+\Gamma_{t}$, $t\geq 1$.

Here, we compare the efficiency of three different MCMC schemes for the update of the parameters $A$. These are the auxiliary gradient-based sampler with interweaving, without interweaving, and a scheme in which instead of interweaving we iterate twice the auxiliary sampler. To simplify our experiment we omit equations $(3.1)$ and $(3.2)$ of the main paper and we simulate observations from the model defined by equations $(3.3)$ and $(3.5)$–$(3.8)$. This a Cox model with intensities driven by latent dynamic factors and Bayesian inference for this model is conducted by switching off certain steps (see Section $4$) of Algorithm $1$. To compare the three sampling schemes we simulated $p=100$ time series consisted of $T=1,500$ observations $n_{it}$, $t=1,\ldots,T$. We used $K=2$ autoregressive latent factors for the jump intensities, we set $\alpha_{1}=0.8$ and $\alpha_{2}=0.4$ and we generated $w_{ik}\overset{iid}{\sim}N(0,1)$. We ran the MCMC algorithms, that correspond to the three different schemes for the update of $F$ and $A$, for $30,000$ iterations and we discarded the first $10,000$ as burn-in period. Figure 12 presents the estimated autocorrelation functions of the parameters in $A$ for each one of the three sampling schemes.

Sampling $\mu_{i\xi}$ and $\sigma^{2}_{i\xi}$

The full conditional distribution of $\mu_{i\xi}$ is

[TABLE]

and the full conditional distribution of $\sigma^{2}_{i\xi}$ is

[TABLE]

Particle filter and annealed importance sampling

Here we present the pseudocode for the implementation of the particle filter and annealed importance sampling algorithms that we constructed in order to estimate predictive Bayes factors. A superscript $s$ in Algorithm 4 means that the operation is performed for all $s=1,...,S$. The infinite sum in (30) is truncated at $10$ without affecting the performance of the algorithm, see also Johannes et al. (2009) for a discussion. The step in the $9$th line of Algorithm 4 is conducted by using multinomial resampling, see doucet2009tutorial for a detailed presentation of resampling methods. The for loop that begins in the $2$nd line of Algorithm 5 is performed in parallel. In the $11$th line of Algorithm 5 we use a few iterations of Algorithm $1$ initialized with the samples obtained from the $(j-1)$th step.

Approximation with a mixture of Gaussian distributions

The approximation of the model defined, in the main paper, by equations $(3.1)$–$(3.3)$ with a mixture of Gaussian distributions relies on the transformation

[TABLE]

and on the assumption that $\log(\epsilon_{it}^{2})|u_{it}\sim\mathcal{N}(m_{u_{it}},s^{2}_{u_{it}})$ where $u_{it}\in\{1,\ldots,10\}$ denotes a mixture component indicator. The means and the variances of the mixture components as well as the probabilities $p(u_{it}=j)$, $j=1,\ldots,10$, can be chosen as suggested by Omori et al. (2007). See Chib et al. (2002) and Nakajima and Omori (2009) for details.

Continuous ranked probability score and mean root squared error

Continuous ranked probability score (CRPS)

Let $F_{r}$ be the cdf of the posterior predictive distribution of $r_{it}$. The CRPS (matheson1976scoring) is defined in terms of the predictive $F_{r}$ and is given by

[TABLE]

where $\tilde{r}_{it}^{1}$ and $\tilde{r}_{it}^{2}$ are independent random variables with distribution $F_{r}$ and it is assumed that the first moment of $F_{r}$ is finite. The CRPS defined by equation (31) is negatively oriented and, thus, a small value indicates a more accurate forecast. Moreover, the negative oriented CRPS reduces to the mean absolute error in the case of a point forecast; see for example Gneiting and Raftery (2007) for more details. Since $F_{r}$ is not available in an analytic form we need to estimate (32) by utilizing the $S$ samples $\{\tilde{r}^{s}_{it}\}_{s=1}^{S}$ drawn from the posterior predictive distribution. To achieve this we first compute the empirical cdf

[TABLE]

Then, an estimation for the CRPS is given by

[TABLE]

We note that in order to compute (34) we utilize the r-package scoringRules (jordan2017evaluating) where the computation is conducted efficiently by employing a representation of (34) which is based on order statistics.

Root mean squared error (RMSE)

To assess point predictions obtained from the different models we calculate the RMSE of the draws from the posterior predictive distribution of the true log-returns $r_{it}$ as

[TABLE]

where $\hat{\tilde{r}}_{it}$ denotes the mean of the posterior predictive distribution.

Results on simulated data

Results for independent SV with jumps models

By using the approach described in Section $3.2$ of the main paper to model the jump intensities in equation $(3.3)$ we obtain $p$ independent univariate SV with jumps models. Bayesian inference for their parameters and latent states can be conducted with the MCMC algorithm derived from Algorithm $1$ by omitting the $10$th step and by substituting the $12$th step (sampling latent factors) with the step of sampling the independent jump intensities from their inverse Gamma full conditionals.

To test our methods we simulated $T=1,500$ log-returns from $p=4$ independent univariate models as defined by equations $(3.1)$–$(3.3)$ with $Gam(1,c)$ priors for the jump intensities. We used $c=50$ to simulate the first and the second and $c=100$ to simulate the third and the fourth time series in order to represent observed log-returns with more and less jumps respectively. Following the related literature (see for example Chib et al. (2002) and Eraker et al. (2003)) we set $\mu_{i}=-0.85$, $\phi_{i}=0.98$ and $\sigma_{i\eta}=0.12$ for each $i=1,2,3,4$ and $\mu_{1\xi}=\mu_{3\xi}=-3$ and $\mu_{2\xi}=\mu_{4\xi}=0$. We ran the MCMC algorithm for $80,000$ iterations from which we discarded the first $20,000$ iterations as burn-in period and we were storing the sampling outcome of every $60$th iteration. Figure 13 presents the simulated log-returns along with the times in which jumps have been generated and times with estimated posterior probability of at least one jump greater than $50\%$. Figure 14 displays the simulated paths for the corresponding log-volatilities and the $95\%$ credible intervals of their posterior distributions. By the visual inspection of the figures it is clear that with the proposed methods we have separated the jumps from the log-volatility process in all the simulated frameworks.

In order to check the ability of the proposed methods to identify jumps in time series with varying volatility of the log-volatility process we conducted the following experiment. We simulated $p=4$ time series, from the model defined by equations $(3.1)$–$(3.3)$, with $T=1,500$ observations in each one by using different values for $\sigma_{i\eta}$, $i=1,\ldots,4$; $\sigma_{1\eta}=0.1$, $\sigma_{2\eta}=0.15$, $\sigma_{3\eta}=0.25$ and $\sigma_{4\eta}=0.35$. We used $Gam(1,50)$ priors for the jump intensities, we set $\mu_{1\xi}=\mu_{2\xi}=0$ and the same values as in the previous experiment for all the other parameters. The same MCMC algorithm as in the first experiment was used to obtain posterior samples for the parameters and the latent states of the model. Figure 15 presents the simulated log-returns for each stock and indicates the times in which we estimated posterior probabilities of at least one jump greater than $50\%$, while Figure 16 displays the posterior probabilities of at least one jump for each one of the $1,500$ time points of each stock. The figures demonstrate that using the proposed methods we are able to disentangle the volatility from the jump process in all the different scenarios for the log-volatility process, since in the majority of the simulated jumps have been identified and only jumps with size that is in the range of the volatility process are not detected.

Results for joint modelled jump intensities

We present results from the application of our methods on $p=100$ time series simulated from the model defined by equations $(3.1)$–$(3.3)$ and $(3.5)$–$(3.8)$ in the main paper. We used $K=3$ latent factors with $\alpha_{1}=0.9$, $\alpha_{2}=0.65$ and $\alpha_{3}=-0.65$ to simulate the paths of the latent factors. The factor loadings in $W$ were simulated from standard normal distributions and we simulated $b_{i}\sim\mathcal{N}(-2.45,1)$, $i=1,\ldots,p$. For the parameters of the latent log-volatility processes we used $\phi_{i}=0.98$, $\sigma_{i\eta}=0.12$ and $\mu_{i}=-0.85$ and we set $\mu_{i\xi}=0$ and $\sigma_{i\xi}=3.5$ for each $i=1,\ldots,p$. Then we ran the MCMC algorithm summarized by Algorithm $1$ for $100,000$ iterations from which we discarded the first $30,000$ as burn-in period and we were saving the outcome of every $70$th iteration.

Each line with black circles in Figure 17 refers to a stock while the circles indicate times in which a jumps has been simulated in the price of the stock. The red crosses indicate the times where the posterior probability of at least one jump has estimated to be greater than $50\%$. Figure 18 presents the posterior distributions for the parameters $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$.

Additional results from the real data analysis

Bayesian analysis of log-volatilities and jumps

Figures 19 and 20 compare the Bayesian inference conducted for the parameters and the latent log-volatility processes of SV models without jumps, with independent jumps and with jointly modelled jump intensities. From the visual inspection of the Figures it is evident that the SV models with jumps, either independent or jointly modelled, deliver lower log-volatility means $\mu_{i}$ and variances $\sigma^{2}_{i\eta}$ and higher persistent parameters $\phi_{i}$ than SV models without jumps while the proposed joint modelling of the jump intensities delivers the most narrow intervals and the SV models without jumps the widest. Figure 21 depicts the posterior probabilities of the number of jumps for each weekday. This is clearly related with the parameter $\Delta_{it}$ in $(3.3)$ which takes into account the time distance of two successive observations.

Root mean squared errors of the predictions

Predictions by using different prior-hyperparameters for the jump intensities

In the case of modelling independently, over time and across stocks, the jump intensities we have assumed that $\lambda_{it}\sim Gam(\delta,c)$ and by following the related literature, see for example Chib et al. (2002), we set $\delta=1$ and $c=50$ in order the prior expectation of $\lambda_{it}$ to be $0.02$ with prior variance $0.0004$. To model jointly the jump intensities we used the dynamic factor model specified by equations $(3.5)$–$(3.8)$ in the main paper. We specified hyperparameters $\sigma_{w}^{2}$, $\sigma_{b}^{2}$ and $\mu_{b}$ such that the induced prior on $\lambda_{it}$ to be similar to the gamma distributions used in the case of independent modelling. Table 2 displays the mean, the variance and the mode of the prior induced, by the dynamic factor model, on each $\lambda_{it}$, for three different choices of the hyperparameters $\mu_{b}$, $\sigma^{2}_{b}$ and $\sigma^{2}_{w}$. Since our aim with joint modelling the jump intensities was to improve the predictions of future observations we perform a sensitivity analysis of the predictions obtained from independent and joint modelled jump intensities with different choices for the prior-hypeparameters of the jump intensities.

Figure 23 compares the predictive performance of univariate SV models without jumps with the performance of SV models with jumps in which the jump intensities are modelled independently over time and across stocks by using the gamma prior distributions presented in Table 3. The figure is constructed by considering the log Bayes factors defined by equation $(5.1)$ of the main paper for the $\ell=30$ out-of-sample observations of our real dataset presented in Section $6$. Figure 24 presents the log Bayes factors in $(5.5)$ and the logarithm of the Bayes factors in $(5.6)$ for the $30$ out-of-sample observations by using the priors on the jump intensities that described in Table 2.

Evaluation of the sequential Monte Carlo methods

To evaluate the accuracy of the estimators calculated by employing sequential Monte Carlo methods such as particle filters and the annealed importance sampling we employ the effective sample size (ESS) of the samples drawn from the posterior predictive distributions of interest. The formula of the ESS is given in the $8$th line of Algorithm 4. We note that a small ESS indicates large variance for the importance weights calculated during the sequential Monte Carlo algorithms which in turn implies that only few of the weighted samples will dominate the estimates; see for example Chopin and Papaspiliopoulos (2020b) for more details.

Figure 25 presents for each out-of-sample time point a boxplot with the $571$ ESSs that correspond to particle filters algorithms while Figure 26 shows the same quantity for the samples drawn by using the annealed importance sampling algorithm. From the visual inspection of the Figures we concluded that for the vast majority of the stocks the ESS of the $1,000$ posterior predictive indicates an accurate estimation of the quantities of interest: predictive Bayes factors, continuous probability ranked scores of the predictive cumulative distribution functions and interval scores of the prediction credible intervals.

Comparison with SV with $t$ errors

An alternative formulation that accounts for the heavy tails of the distribution of the observed returns is the SV model in which the observation errors follow the Student’s $t$ distribution Harvey et al. (1994). Bayesian inference for the parameters of SV models with Student’s $t$ errors has been examined by Chib et al. (2002). Figure 27 compares the predictive performance of univariate SV models with Student’s $t$ errors with the performance of the proposed SV model with jumps and jointly modelled jump intensities. The Figure evaluates the predictive performance of the two models with respect to their interval scores, CRPS as well as root mean squared errors. From the visual inspection of the Figure we conclude that by modelling jointly the jump intensities of SV models we obtain more accurate forecasts rather than using Student’s $t$ instead of Gaussian errors. To conduct Bayesian inference for the parameters and the latent states of the SV models with Student’s $t$ errors we utilized the r-package stochvol (kastner2014dealing) to collect $1,000$ thinned MCMC samples by storing the outcome of every $100$th iteration after the first $50,000$ iterations. To draw samples from the posterior predictive distributions of interest we employed the particle filter algorithm summarized by Algorithm 4 by setting $p(r_{it}|h_{it})$ in (30) to be the density of the Student’s $t$ distribution scaled by $e^{\tfrac{1}{2}h_{it}}$.

Computing times