Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
96,719 characters · 15 sections · 83 citation commands
A Statistical Recurrent Stochastic Volatility Model for Stock Markets
The volatility of a financial time series, such as stock returns, is defined as the variance of the returns and serves as a measure of the uncertainty about the returns. The volatility, which is of great interest to financial econometricians, is unobserved and therefore often modelled statistically in order to estimate it. The two model classes most frequently used in volatility modelling are the Generalized Autoregressive Conditional Heteroscedastic (GARCH) models and the Stochastic Volatility (SV) models. The GARCH model Bollerslev1986 expresses the current volatility, conditional on the previous returns and volatilities, as a {\it deterministic} and linear function of the squared returns and the conditional volatilities in the previous time period. The SV model Taylor:1982,Taylor:1986, on the other hand, uses a latent stochastic process to model the volatility, which is usually taken as a first order autoregressive process. It is well documented that the GARCH and SV models are able to capture important effects exhibited in the variance of financial returns. For example, the volatilities in financial returns are observed to be highly autocorrelated in certain time periods and exhibit periods of both low and high volatility Benoit1967. This so-called volatility clustering phenomenon can be modeled by the volatility processes introduced in the GARCH and SV models, making these volatility models widely employed in financial time series modelling.
Although the GARCH and SV models were independently and almost concurrently introduced, the GARCH models were initially more widely adopted as it is much easier to estimate GARCH models than SV models. This is because the likelihood of a GARCH model can be obtained explicitly, while the likelihood of a SV model is intractable as it is an integral over the latent process. However, the conditional variance process of GARCH models is deterministic and hence GARCH models might not capture efficiently the random oscillatory behavior of financial volatility Nelson:1991. SV models are considered as an attractive alternative to GARCH models because they overcome this limitation Kim:1998,Yu:2002. Recent advances in Bayesian computation such as particle Markov chain Monte Carlo (PMCMC) Andrieu:2010 allow straightforward estimation and inference for SV models.
Standard SV models still cannot appropriately capture some important features arising in financial volatility. For example, a large amount of both theoretical and empirical evidence indicates that there exists long-range persistence in the volatility process of many financial returns, see, e.g, Lo:1991, Ding:1993, Crato:1994 and Bollerslev:1996. The long-memory property of a time series implies that the decay of the autocorrelations of the series is slower than exponential. The standard SV model of Taylor:1982 uses an AR(1) process to model the log of the volatility and hence might fail to capture this type of persistence Breidt:1998. Another line of the literature shows strong evidence of non-linear auto-dependence in the volatility process of some stock and currency exchange returns KILIC2011368 and that the simple linear AR(1) process cannot effectively capture the underlying non-linear volatility dynamics.
Breidt:1998 proposed the Long Memory Stochastic Volatility (LMSV) model to overcome the short-memory limitation of the standard SV model. LMSV uses an ARFIMA process Granger:1980 as an alternative to the AR(1) process to capture the long-memory dependence in the volatility. The empirical evidence in Breidt:1998 suggests that the LMSV model is able to capture the long-memory volatility behaviour in some stock return datasets. However, the literature is unclear about whether the LMSV model can capture non-linear dynamics within the volatility process, because the ARFIMA model is linear. Additionally, it is challenging to estimate the LMSV model as its likelihood is intractable. We are unaware of any available software package that implements the LMSV methodology. In another approach, YuJun:2006 introduced a family of non-linear SV (N-SV) models to capture the possible departure from the log transform commonly used in SV models. In the standard SV model, the logarithm of volatility is assumed to follow an AR(1) process; N-SV uses other non-linear transformations, such as the Box-Cox power function, rather than the logarithm. The simulation studies and empirical results on currency exchange and option pricing data in YuJun:2006 show that the N-SV model using the Box-Cox transformation is able to detect some interesting effects in the underlying volatility process. The general use of N-SV models requires the user to select an appropriate non-linear transformation for the dataset under consideration, and this might lead to a challenging model selection problem. Neither Breidt:1998 nor YuJun:2006 clearly discussed the out-of-sample forecast performance of their LMSV and N-SV models.
Recurrent neural networks (RNN) in the Deep Learning literature have impressive prediction performance and have been successfully deployed in a large number of industrial-level applications (language translation, image captioning, speech synthesis, etc.). The RNN models are well-known for their ability to efficiently capture the long-range memory and non-linear dependence existing within various types of sequential data, and are considered as the state-of-the-art models for many sequence learning problems Lipton:2018. Many researchers and practitioners have used RNN for mean modelling in financial time series analysis, but the general consensus is that these machine learning models do not clearly outperform the traditional time series models such as ARMA and ARIMA (see, e.g., Makridakis:2018 and ZHANG2003159). Makridakis:2018 note that without careful modifications, Machine Learning models are usually less accurate than the statistical approaches that have been extensively investigated in the financial time series literature. Recently, the idea of using the RNN models to improve the predictive performance of GARCH-type models has also been proposed for volatility modelling. For example, Kim:2018 use the volatility estimates from several GARCH-type models as inputs to a RNN model, which then non-linearly transforms these inputs to output the final estimate of the volatility. The empirical results on the Korean stock market KOPSI 200 index show a significant improvement of forecast performance of the proposed hybrid model over several GARCH-type benchmark models. However, similar to many engineering-oriented Machine Learning models, Kim:2018's model overlooks the interpretation aspect in volatility modelling, which is often of main interest to econometricians. One of the main motivations of our article is to develop deep learning based volatility models that are not only able to produce accurate prediction, but also interpretable and have meaningful in-sample analysis. These models should not overlook the well-established features of traditional econometric models, that are motivated by the well-known stylized facts in financial time series such as volatility clustering and fat tails.
In the SV literature, there is still lack of research using RNN structures to model the stochastic volatility dynamics of financial time series, perhaps because of two reasons. First, it is non-trivial to sensibly incorporate RNN into the statistical volatility models. Simple adaptations of RNN to volatility models easily overlook the important stylized facts exhibited in financial volatility, which are well captured by the AR(1) process in the SV model. It is important to select appropriate RNN structures that are not only able to produce accurate out-of-sample volatility forecast, but also explain well the volatility dynamics. Second, a stochastic volatility model that incorporates a RNN structure into its latent stochastic process is highly sophisticated and thus challenging to estimate.
This paper combines the SV and RNN models in a non-trivial way, and proposes a new model, called the Statistical Recurrent Stochastic Volatility (SR-SV) model. In particular, we use the Statistical Recurrent Unit (SRU) structure of Oliva:2017, which is a special type of RNN models, to capture complex volatility effects overlooked by an AR(1) process in the standard SV model but still retain the essential components of the SV model. This combination allows the SR-SV model to enjoy much of advances from both worlds of deep learning (e.g., flexibility and excellent predictive performance) and econometric volatility modelling (e.g., excellent interpretability of volatility effects). The SR-SV model belongs to the class of parametric state space models whose Bayesian inference can be performed using recent advances in the Sequential Monte Carlo (SMC) and particle MCMC literature Andrieu:2009,Andrieu:2010,Duan:2015,Deligiannidis:2018. The simulation studies and empirical results on the five stock index datasets demonstrate that the SR-SV model can efficiently capture the potential non-linear and long-memory effects in the underlying volatility dynamics, and provide better out-of-sample forecasts than the standard SV, N-SV and LMSV models. We note that we have tested SR-SV on a wider range of stock returns but only report in the paper the results for five of them, as we constantly observed a similar improvement of the model compared to the other three counterpart. A Matlab software package implementing Bayesian estimation and inference for SR-SV together with the examples reported in this paper are available on github\footnote{The link is provided in the unblinded version.}.
The article is organized as follows. Section (ref) briefly reviews the SV and SRU models, and presents the SR-SV model. Section (ref) discusses in detail Bayesian estimation and inference for the SR-SV model. Section (ref) presents the simulation study and applies the SR-SV model to analyze the five stock index datasets. Section (ref) concludes. The Appendix gives details of the implementation and further empirical results.
Let $y=\{y_t,\ t=1,...,T\}$ be a series of financial returns. We consider a basic version of SV models Taylor:1982 \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu The persistence parameter $\phi$ is assumed to be in $(-1,1)$ to enforce stationarity of both the $z$ and $y$ processes. In this SV model, the log volatility process $z$ is assumed to follow an AR(1) model. It is well documented in the financial econometrics literature that financial time series data often exhibit a long-term auto-dependence, which forces the persistence parameter $\phi$ to be close to 1 Jacquier:1994,Kim:1998. Write $p(z|\theta)$ for the density of $z$ given the model parameters $\theta=(\mu,\phi,\sigma^2)$ and $p(y|z)$ for the density of the data $y$ conditional on $z$. We can view $p(z|\theta)$ as the prior with $\theta$ being the hyper-parameters and $p(y|z)$ as the likelihood Jacquier:1994. Under this perspective, the SV model (ref)-(ref) puts non-zero prior mass on AR(1) stochastic processes, and zero or almost-zero mass on stochastic processes that are far from being well approximated by an AR(1). This means that the SV model in (ref)-(ref) might not be able to capture more complex dynamics in the posterior behavior of the log volatility process $z$, such as long-term memory or non-linear auto-dependence, and that a more flexible prior distribution should be put on $z$. We will design such a flexible prior by combining the attractive features from both SV and RNN time series modeling techniques.
YuJun:2006 propose a class of non-linearity N-SV models as a variant of SV which allows a more flexible link between the variance $\text{\rm Var}(y_t|z_t)$ and the AR(1) process $z$. Their N-SV model, using the Box-Cox transformation for $\text{\rm Var}(y_t|z_t)$, is written as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu where $\delta$ is the auxiliary parameter that measures the degree of non-linearity rather than the log transform. As $\delta \rightarrow 0$, $(1+\delta z_t)^{1/2\delta} \rightarrow e^{\frac12z_t}$ and hence the N-SV model includes the SV model as a special case. The term non-linearity here might cause some confusion, as it does not refer to the non-linear auto-dependence within the log volatility process $z$, but the non-linearity between $\text{\rm Var}(y_t|z_t)$ and $z_t$.
Breidt:1998 suggest to use an ARFIMA$(p,d,q)$ process Granger:1980,Hosking:1981 for the log volatility $z_t$ to capture the long-memory auto-dependence exhibited in financial time series. Their LMSV model is written as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu where $\Phi(B) = 1-\phi_1 B - \phi_2 B^2-...- \phi_p B^p$, $\Theta(B) = 1+\theta_1 B +\theta_2 B^2+ ... + \theta_q B^q$, and $B$ is a backshift operator, i.e., $B^s X_t = X_{t-s}$. To ensure the stationarity and invertability of the log volatility process $z_t$, the fractional integration parameter $d$ is assumed to be in $(-0.5,0.5)$ and the roots of $\Phi(B)$ and $\Theta(B)$ have to lie outside the unit circle.
Another notable line of research in the SV literature is the class of semi-parametric stochastic volatility models that incorporate non-parametric techniques into modelling the conditional distribution of financial returns. For example, the stochastic volatility, Dirichlet process mixture (SV-DPM) model of Jensen:2010 uses a Dirichlet process prior Ferguson:1973 to characterize the conditional distribution of $y_t$. Semi-parametric models are different from the SR-SV model in two important aspects. First, the SV-DPM model is proposed to capture the asymmetries and leptokurtotic behaviors of financial returns, while the SR-SV model focuses on modeling the non-linearity and long-memory auto-dependence in the log-volatility dynamics. Second, the SV-DPM model is a semi-parametric model in the sense that the model cannot be described using a finite number of parameters as it uses a non-parametric prior, e.g. Dirichlet process, to simulate the conditional return and retains the parametric structure, e.g. AR(1), of the log-volatility in the standard SV model. The SR-SV model, on the other hand, is a parametric model with eleven parameters whose mathematical representation will be discussed in Section (ref). Our article therefore uses parametric models including the standard SV, N-SV and LMSV models as the benchmarks to evaluate the SR-SV model.
There are at least two approaches to modeling time series data. One approach is to represent time effects {\it explicitly} via some simple function, often a linear function, of the lagged values of the time series. This is the mainstream time series data analysis approach in the statistics literature with the well-known models such as AR or ARMA. The alternative approach is to represent time effects {\it implicitly} via latent variables, which are designed to store the memory of the dynamics in the data. These latent variables, also called hidden states, are updated in a recurrent and deterministic manner using the information carried over by their values from the previous time steps and the information from the data at the current time step. Recurrent neural networks (RNN), belong to the second category, were first developed in cognitive science and successfully used in computer science and other fields. Another class of models that represent time implicitly is state space models, albeit the recurrent update is stochastic, which are widely used in econometrics and statistics. The SV model discussed in Section (ref) is an example of state space models.
For the purpose of this section, we denote the time series data as $\{D_t=(x_t,z_t),\ t=1,2,...\}$ where $x_t$ is the vector of inputs and $z_t$ the scalar output. In our article, it is useful to think of $x_t$ as scalar; however, the RNN approach is often efficiently used to model multivariate time series. If the time series of interest has the form $\{z_t,\ t=1,2,...\}$, it can be written as $\{(x_t,z_t),\ t=2,...\}$ with $x_t=z_{t-1}$. Our goal is to model the conditional distribution $p(z_t|x_t,D_{1:t-1})$. If the serial dependence structure is ignored, then a feedforward neural network (FNN) can be used to transform the raw input data $x_t$ into a set of hidden units $h_t$, also called {\it learned features} or {\it summary statistics}, for the purpose of explaining or predicting $z_t$. However, this approach is unsuitable for time series data as the time effects or the serial dependence are totally ignored. The main idea behind RNN is to let the set of hidden units $h_t$ to feed itself using its lagged value $h_{t-1}$ from the previous time step $t-1$. Hence, RNN can be best thought of as a FNN that allows a connection of the hidden units to their value from the previous time step, enabling the network to possess memory. Mathematically, this RNN model Elman:1990 is written as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu The model parameters include $w_x$, $w_h$, $b$, $\beta_0$ and $\beta_1$, $\Psi(\cdot)$ is a non-linear activation function, e.g., common choices are the {\it sigmoid} $\Psi(z)=1/(1+e^{-z})$ and the {\it tanh} $\Psi(z) = (e^{z} - e^{-z})/(e^{z} +e^{-z})$, and $p(z_t|\eta_t)$ is a probability density depending on the learning task. For example, if $z_t$ is continuous, then typically $p(z_t|\eta_t)$ is a Gaussian density with mean $\eta_t$; if $z_t$ is binary, then $z_t|\eta_t$ follows a Bernoulli distribution with probability $\Psi(\eta_t)=1/(1+e^{-\eta_t})$. Usually one sets $h_1=0$, i.e. the neural network initially does not have any memory.
Figure (ref) illustrates graphically the RNN model (ref)-(ref). We follow Goodfellow:2016 and use a black square to indicate the delay of a single time step in the circuit diagram (left). The circuit diagram can be interpreted as an unfolded computational graph (right), where each node is associated with a particular time step.
The unfolded graph in Figure (ref) suggests that the hidden state at time $t$ is the output of a composite function
which somewhat resembles a multiplication structure in terms of the weight $w_h$. Consequently, the gradient of $h_t$ with respect to the model parameters might either explode or vanish if $t$ is sufficiently large and $w_h$ is not equal to 1, and hence making it inefficient for the Simple RNN model to learn in long time series. See Goodfellow:2016 for further explanation.
Many sophisticated RNN structures have been proposed to overcome the aforementioned problem in the Simple RNN model; for example, the Long Short-term Memory model of Hochreiter1997, the Gated Recurrent Unit of cho-etal-2014-learning and the Statistical Recurrent Unit (SRU) of Oliva:2017. The SRU allows the vector of summary statistics $h_t$ to traverse through the network using a moving average. We will use the SRU in this paper as its structure and some of its main parameters carry statistical meaning; see Section (ref). A general SRU structure is mathematically written as
where $\alpha=(\alpha_1,...,\alpha_m) \in (0,1)$ is a vector of moving average weights, and $W_h$, $b_r$, $W_{r}$, $W_x$ and $b_{\varphi}$ are the model parameters. We denote the functional learning structure in (ref)-(ref) as $h_t=\text{SRU}(x_t,h_{t-1})$, which takes $x_t$ - the input data at current time $t$ - and $h_{t-1}$ - the previous output of the SRU - as the input arguments. See Figure (ref)(a) for the graphical representation of this SRU structure. The moving average structure of the state $h_t$ allows the RNN network with SRU units to enjoy some advantages compared to other RNN models. The current state $h_t$ is related to the previous state $h_{t-1}$ both {\it directly} and {\it indirectly} and hence mitigate the problem of multiplying the same quantities multiple times as in the Simple RNN model. The novel architecture of the SRU allows the model to capture long term dependencies in data via simple moving averages.
This section proposes the SR-SV model that combines SV and SRU for financial volatility modelling. The key idea is that we use the SRU structure to capture the complicated effects such as long-term memory and non-linear auto-dependence, in the volatility dynamics that are overlooked by the basic SV models. This leads to a prior distribution for the log volatility process $z$ that is much more flexible than the AR(1) prior (c.f. Section (ref)). Our proposed SR-SV model is as follows \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu that is, we use a SRU to model the dynamics of the hidden states $h_t$. Here, $z_0$ is the initial value of the log volatility process and a convenient choice of $z_0$ is the log of the unconditional variance of the observed series $y$, i.e., $z_0 = \log(\text{var}(y))$. We follow the literature to initialize $h_1 = 0$ as the recurrent units initially have no memory. Figure (ref) plots the graphical representation of the SR-SV model. See Appendix (ref) for the fully-written version of the SR-SV model.
We note the following important properties of the SR-SV model. First, the SR-SV model in (ref)-(ref) retains the measurement equation (ref) and the linear part $\phi z_{t-1}$ of the AR(1) process from the standard SV model, and captures the volatility effects not captured by the AR(1) process, e.g. non-linear and long-memory auto-dependence, via the latent state $h_t$ of the SRU structure. The log volatility at time $t$ in (ref) can be written as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu Therefore, the parameter $\beta_1$ characterizes all the effects in the underlying log volatility process $z$ rather than the short-term linear effect captured by the AR(1) process. We refer to $\beta_1$ as the {\it non-linearity long-memory} coefficient. If $\beta_1=0$ and $\epsilon_1^\eta{\sim}{\cal N}(\beta_0/(1-\phi),\sigma^2/(1-\phi^2))$, the SR-SV model becomes the SV model (ref)-(ref) and hence the SV model is a special case of the SR-SV model. We therefore follow the SV literature and assume that $|\phi|<1$. The $z$ process, and thus the $y$ process of the SR-SV model, is not guaranteed to be stationary unless $\beta_1=0$ and $\epsilon_1^\eta{\sim}{\cal N}(\beta_0/(1-\phi),\sigma^2/(1-\phi^2))$. Non-stationarity for volatility is often argued to be more realistic in practice (e.g. Bellegem:2012), although it might be mathematically less appealing. The equation in (ref) can be further written out as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu where $\mathfrak{N}(\cdot)$ is a non-linear function and $w_{z}$ is the weight corresponding to $z_{t-1}$; see the full version of the SR-SV model in the Appendix. If $w_{z}=0$ in (ref), then $z_t$ only depends linearly on $z_{t-1}$, therefore this equation indicates that the parameter $w_{z}$ characterizes the serial dependence rather than linearity that the previous log volatility $z_{t-1}$ has on $z_t$. We will analyse $w_{z}$ in more details in Section (ref).
Second, Oliva:2017 set the scales $\alpha$ of the SRU model to several pre-specified values to obtain a vector of summary statistics $h_t$ at different moving average weights. We, however, treat $\alpha$ as a model parameter and learn it from the data. We note that a higher $\alpha$ weighs more on the historical information while a smaller $\alpha$ puts a more weight on the current information. We show later in the empirical study that this parameter $\alpha$ is able to quantify the existence of the long-memory auto-dependence commonly exhibited in the volatility dynamics of the financial time series.
Third, neural networks are highly flexible but often subject to overfitting, i.e., they have over-confident in-sample fit and bad out-of-sample forecasts. Regularization is often needed to avoid overfitting. Injecting noise into the layers of the network has been found an effective regularization approach in the Machine Learning literature, and seen as a form of data augmentation at multiple levels of abstraction Sietsma:1991,poole:2014,Goodfellow:2016,dieng:2018. In the SR-SV model, by allowing $z_{t-1}$ and $\eta_{t-1}$ to be the inputs of the SRU structure at time $t$, we inject the noise $\epsilon^\eta_{t-1}$ of the volatility process to the input and hidden layers of the SRU. This noise-injecting regularization approach makes the SR-SV model perform well on both in-sample fitting and out-of-sample forecast, even with the simplest specification of the SRU structure where all the $r_t$, $\varphi_t$ and $h_t$ are scalars. Our SR-SV model can be categorized as a {\it parametric} model with the vector of model parameters $\theta$ consisting of eleven parameters: four main parameters $\beta_0$, $\beta_1$, $\phi$, $\sigma^2$ and the parameters in the SRU including $\alpha, w_h,b_r,w_r,w_\eta,w_z$ and $b_\varphi$.
Finally, $\beta_0$ plays the role of the scale factor $\tau=e^{\beta_0/2}$ for the variance of $y_t$. One could set $\beta_0=0$ and modify (ref) to $y_t = \tau e^{\frac12z_t}\epsilon^y_t$; however, this parameterization might be less statistically efficient in terms of Bayesian estimation, especially for the parameter $\tau$ Kim:1998.
It is straightforward to extend the SR-SV model in (ref)-(ref) by incorporating other advances in the SV literature. For example, we can use a Student's $t$ distribution instead of a Gaussian for the measurement shock $\epsilon^y_t$ and take into account the leverage effect by correlating $\epsilon^y_t$ with the volatility shock $\epsilon_t^\eta$. We do not consider these extensions here, however, because using the most basic version makes it easier to understand the strengths and weaknesses of the new model.
This section discusses Bayesian estimation and inference for the SR-SV model. For a generic sequence $\{x_t\}$ we use $x_{i:j}$ to denote the series $(x_i,...,x_j)$. The SR-SV model is a state-space model with the measurement equation \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu and the state transition equation \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu We are interested in sampling from the posterior distribution of $\theta$ \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu where $p(y_{1:T}|\theta)$ is the likelihood function, $p(\theta)$ is the prior and $p(y_{1:T}) = \int_{\Theta}p(y_{1:T}|\theta)p(\theta)d\theta$ is the marginal likelihood. Recall that the vector of model parameters $\theta$ consists of $\beta_0$, $\beta_1$, $\phi$, $\sigma^2$ and the 7 parameters within the SRU model (ref)-(ref).
The likelihood function in ((ref)) is \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu which is computationally intractable for non-linear non-Gaussian state space models like the SV and SR-SV models, but can be estimated unbiasedly by a particle filter DelMoral:2004. Bayesian inference for SR-SV can be performed using recent advances in the Sequential Monte Carlo literature that we present next.
Duan:2015 propose the Density Tempered Sequential Monte Carlo (DT-SMC) approach to Bayesian inference for state space models where the likelihood is intractable. The DT-SMC sampler generalizes the SMC method of Neal:2001 and DelMoral:2006 when the likelihood can be computed analytically. In order to sample from the posterior $\pi(\theta)$, the DT-SMC method first samples a set of $M$ weighted particles $\{W^j_0,\theta_0^j\}^M_{j=1}$ from an easy-to-sample distribution $\pi_0(\theta)$, such as the prior $p(\theta)$, and then traverses these particles through intermediate distributions $\pi_t(\theta), \;\; t=1,...,K$, which target the posterior distribution $\pi(\theta)$ eventually, i.e. $\pi_K(\theta)=\pi(\theta)$. The DT-SMC method uses the following intermediate distributions \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu where the $\gamma_t$ is referred to as the level temperature and $0 = \gamma_0 < \gamma_1 < \gamma_2 < ... < \gamma_K=1$, $\widehat{p}(y_{1:T}|\theta,u)$ is the unbiased estimator of the likelihood $p(y_{1:T}|\theta)$ and $u$ is the set of pseudo random numbers used within a particle filter to estimate the likelihood $p(y_{1:T}|\theta)$. For the purpose of this paper where it is possible to sample from the prior $p(\theta)$, we set $\pi_0(\theta)=p(\theta)$. Algorithm (ref) summarizes the DT-SMC method for the SR-SV model.
The DT-SMC method consists of three main steps: reweighting, resampling and Markov move. At the begining of SMC iteration $t$, the set of weighted particles $\{W_{t-1}^j,\theta_{t-1}^j\}^M_{j=1}$ that approximate the intermediate distribution $\pi_{t-1}(\theta)$ is reweighted to approximate the target $\pi_{t}(\theta)$. The efficiency of these weighted particles is often measured by the effective sample size (ESS) Robert:1998,LiuChen:1998 defined in (ref). If the ESS is below a prespecified threshold, the particles are resampled; the resulting equally-weighted resamples, which are now approximate samples from $\pi_{t}(\theta)$, are then refreshed by a Markov kernel whose invariant distribution is $\pi_{t}(\theta)$. For example, Duan:2015 uses the pseudo marginal Metropolis-Hastings (PMMH) kernal of Andrieu:2010 with the likelihood estimated unbiasedly by the particle filter in the Markov move step. However, Pitt:2012 suggest that the PMMH approach works efficiently when the variance of the log of the estimated likelihood is around 1. For some state space models like the SR-SV model, a large number of particles might be required to obtain a likelihood estimator with log variance to be around 1, which is computationally inefficient. To tackle this problem, we incorporate the Correlated Pseudo Marginal (CPM) approach of Deligiannidis:2018 into the Markov move step. The CPM method makes the current set of random numbers $u$ and proposal $u^\prime$ correlated, and helps reduce the variance of the ratio $\widehat{p}(y_{1:T}|\theta^{\prime},u^{\prime}) / \widehat{p}(y_{1:T}|\theta,u)$ in (ref), thus leading to a better mixing Markov chain while using less number of particles in the particle filter. Similar to the SMC methods of DelMoral:2006 and Neal:2001, the DT-SMC method is parallelizable as the particles move independently in the Markov move step, and provides an estimate of the marginal likelihood as a by-product.
In Algorithm (ref), we use a random walk proposal for $q(\theta^\prime|\theta)$. We follow Dang:2018 and choose the tempering sequence $\gamma_t$ adaptively to ensure a sufficient level of particle efficiency by selecting the next value of $\gamma_t$ such that ESS stays above a threshold.
The marginal likelihood is often used to choose between models via the Bayes factor jeffreys:1935,Kass:1995. In order to compare the relative performance between two models $M_1$ and $M_2$ on a given data $y_{1:T}$, we can use the Bayes factor defined by \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu providing a Bayesian alternative to hypothesis testing. The larger the Bayes factor $F_{M_1,M_2}$, the stronger $M_1$ is supported by the data than $M_2$. Jeffreys:1961 suggests a scale of interpretation of the Bayes factor $F_{M_1,M_2}$ as listed in Table (ref). We note that the DT-SMC sampler in the previous section provides an efficient way to compute the marginal likelihood.
This section evaluates the performance of the SR-SV model relative to the SV, N-SV and LMSV models using a simulation study and real data applications. We do not report the results for GARCH as it performs similarly to SV. We use the DT-SMC sampler for the Bayesian inference in the SV, N-SV and SR-SV models. As the LMSV model does not have an explicit state-space representation and its likelihood function is analytically intractable, we follow Breidt:1998 and estimate the LMSV model on the frequency domain. The implementation of Bayesian inference for the LMSV model is presented in Appendix (ref). Table (ref) in Appendix (ref) lists our implementation details of the DT-SMC sampler. All the computations are run on a High-Performance Computing (HPC) machine with 16-core CPU and 16 Gigabytes of RAM. The DT-SMC sampler was initialized by sampling from the priors in Table (ref).
We now motivate the choice of the priors in Table (ref). We follow YuJun:2006 and Kim:1998 to set the same prior, which is a Beta distribution, for the persistence parameters $\phi$ of the three models SV, N-SV and SR-SV. We also use an inverse-Gamma prior for the parameters $\sigma^2$ in all models but make it more flat than the priors used in YuJun:2006 and Kim:1998. We follow YuJun:2006 to use an informative but reasonably flat prior distribution for the intercept $\mu$ in the SV and N-SV models. For the SR-SV model, we found that the posterior distributions of $\beta_1$ and $w_{z}$ to be unimodal under inverse-Gamma priors. We use a normal prior with a zero mean and a small variance for the SRU parameters, except $w_{z}$, because empirical results from the SRU literature show that the values of the SRU parameters are often small. Finally, we set a normal prior with a zero mean and a small variance for the intercept $\beta_0$ in the SR-SV model as the empirical results often show small values of $\beta_0$.
Table (ref) lists the predictive scores used to measure the out-of-sample performance. The smaller the predictive scores, the better.
We consider three volatility models specified in Table (ref). Model 1 is a GARCH(1,1). Model 2 is an extension of the GARCH(1,1) model by applying a Box-Cox power transformation Box:1964 to both the conditional variance equation and the volatility dynamics. Model 2 is similar to the non-linear ARCH models of Higgins:1992 but they use lagged innovations to construct the conditional variance. Model 3 is a FIGARCH(1,$d$,1) model of Baillie:1996 which uses a long-memory process AFRIMA(1,$d$,1) to simulate the long-memory auto-dependence.
We generate time series of $T=3000$ observations from these three models and name the simulation datasets as SIM I, SIM II, SIM III, accordingly. The parameters are set so that $y_t$ somewhat resembles real financial time series data exhibiting volatility clustering with non-linearity (SIM II) and long-memory (SIM III) auto-dependence in the underlying volatility dynamics. For each dataset, the first $T_{in} = 2000$ observations are used for model estimation and the last $T_{out} = 1000$ are for out-of-sample analysis. Table (ref) shows the posterior mean estimates for the parameters of the SV and SR-SV models, with the posterior standard deviations in brackets; for the SR-SV model we only show the results for the main parameters. The last column in the table shows the marginal likelihood estimates, averaged over 10 different runs of the DT-SMC sampler, together with the Monte Carlo standard errors in the brackets. Figures (ref), (ref) and (ref) plot the filtered values of the $\eta_t$ and $h_t$ components of the SR-SV model in SIM I, SIM II and SIM III datasets, respectively. Figure (ref) in Appendix (ref) plots the true volatility together with the filtered volatility produced by the SV and SR-SV models, for all the three simulation datasets.
The estimation results suggest the following conclusions. First, for the SIM I data, the difference of marginal likelihood estimates in Table (ref) between the SV and SR-SV models is insignificant, the coefficient $\beta_1$ is insignificant, and the filtered volatilities from these two models in Figure (ref) are identical and close to the true volatility. This implies that the SV and SR-SV models fit equally well to the SIM I data and that the SR-SV model is close to the SV model if the true data generating process, which is GARCH(1,1) in this example, exhibits no other effects rather than short-memory linear auto-dependence within the volatility dynamics.
Second, the estimation results on the SIM II and SIM III data show that the additional neural network structure of the SR-SV model is able to efficiently capture the volatility effects overlooked by the basic SV model. This is supported by the Bayes factors of the SR-SV model compared to the SV model of more than $e^{2.3}$, which, according to the interpretation in Table (ref), strongly support the SR-SV model. The plots of the $h_t$ and $\eta_t$ components of the SR-SV model in Figures (ref), (ref) and (ref) clearly show that $h_t$, and hence $\eta_t$, is well responsive to volatility effects rather than the linear short-memory effects. For example, in the SIM I data when the volatility exhibits no non-linear and long-memory effects, the $h_t$ shown in Figure (ref) is significantly small at all time steps and hence the $\eta_t$ component simply fluctuates around $\beta_0$ during both low and high volatility periods. Figures (ref) and (ref), on the other hand, show that the $h_t$ response adaptively to the changes in the volatility dynamics. As the result, $\eta_t$ is small during the low volatility periods and large in the high volatility periods. The non-linear (SIM II) and long-memory (SIM III) auto-dependence of the simulated volatility are well captured by the SRU structure of the SR-SV model. The plots of filtered volatility in Figure (ref) show that the filtered volatility of the SR-SV model are generally closer to the true volatility than those of the SV model.
Third, the parameters of the SR-SV model are able to characterize well the existence of the various volatility effects in these simulation data. The estimated posterior means of parameter $\beta_1$ are more than two standard deviations from zero in the SIM II and SIM III data, suggesting the existence of volatility effects rather than linearity in the volatility dynamics of these two datasets; while $\beta_1$ is less than two standard deviations from zero in the SIM I data, suggesting that only simple linear effects are detectable in the volatility dynamics of this dataset. Similarly, the non-linear coefficient $w_{z}$ of the SR-SV model (c.f. equation (ref)) is more than two standard deviations from zero in the SIM II and SIM III data but not in the SIM I data, indicating that the SR-SV model is able to detect the serial dependence rather than the linear dependence that the past log volatility $z_{t-1}$ has on the current log volatility $z_t$. The estimated posterior mean of the moving average weight parameter $\alpha$ in the SIM III data is higher than those in the SIM I and SIM II data, supporting further the evidence of the long-memory auto-dependence exhibited in the volatility dynamic of the SIM III data, which is generated from a FIGARCH(1,$d$,1) model. Figure (ref) shows that the posterior mode of $\alpha$ in the SIM III data is much closer to 1 than those in the SIM I and SIM II data. Finally, it is worth noting that the persistence parameter $\phi$ of the SR-SV is close to the persistence parameter $\phi$ of the SV model in the SIM I data but much smaller than that of the SV model in the SIM II and SIM III model. This is probably because the non-linear coefficient $w_z$ with respect to the past log volatility $z_{t-1}$ is significant in the SIM II and SIM III data, and hence the historical information has been also well stored in the $\eta_t$ process.
Table (ref) reports the predictive performance scores of the SV and SR-SV models with the Monte Carlo standard errors in brackets. For the SIM II and SIM III data, the SR-SV model outperforms the SV model for all the predictive scores, which is consistent with the in-sample analysis showing that the SR-SV model fits these simulation datasets better than the SV model. For the SIM I data, the SR-SV model also outperforms the SV model for all scores except the PPS score. These results illustrate the impressive out-of-sample forecast ability of the SR-SV model. The results for the real data applications in the next section further support this claim.
This section evaluates the SR-SV model using five popular daily stock indexes from different international markets: The German stock index DAX30 (DAX), the Hong Kong stock index HS50 (HSI), the France market index CAC40 (FCHI), the US stock market index SP500 (SPX) and the Canada market index TSX250 (TSX).
The datasets were downloaded from the Realized Library of The Oxford-Man Institute\footnote{https://realized.oxford-man.ox.ac.uk/}. We used the adjusted closing prices $\{P_t, \ t=1,...,T_P\}$ and calculated the demeaned return process as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu and using the first $T_{\text{in}}=2000$ returns for in-sample analysis and the rest $T_{\text{out}}=1000$ for out-of-sample analysis. Table (ref) describes the relevant aspects of the datasets.
Figure (ref) in Appendix (ref) plots the time series data and shows the existence of the volatility clustering effect commonly seen in financial data. Table (ref) reports some descriptive statistics together with Lo's modified R/S test Lo:1991 for long-range memory in the absolute and squared returns. Lo's modified R/S test is widely used in the financial time series literature; see, e.g., Lo:1991, Giraitis:2003, Breidt:1998. All the index data exhibit some negative skewness, high excess kurtosis and high variation. The result of Lo's modified R/S test for long-memory dependence with several different lags $q$ indicates that there is significant evidence of long-memory dependence in the stock indices.
The Realized Library provides different realized measures\footnote{See https://realized.oxford-man.ox.ac.uk/documentation/estimators for the list of the available realized measures} that can be used in financial econometrics as a proxy to the latent $\sigma_t^2$. We use the following four common realized measures including Realized Variance (RV) Andersen:1998, Bipower Variation (BV) Nielsen:2004, Median Realized Volatility (MedRV) Andersen:2012, Realized Kernel Variance Barndorff&Nielsen:2008 with the Non-Flat Parzen kernal ($\text{RKV}$) to evaluate the forecast performance of the volatility models using the predictive scores in Table (ref). See Shephard:2010 for more details about the Realized Library.
Denote by $RV_t$ the realized measure of $\sigma_t^2$ at time $t$. As the realized measures ignore the variation of the prices overnight and sometimes the variation in the first few minutes of the trading day when recorded prices may contain large errors Shephard:2010, we follow Hansen:2005 to scale the realized measure $RV_t$ as \thinmuskip= 3mu plus 0mu minus 2mu \medmuskip= 4mu plus 2mu minus 2mu \thickmuskip=5mu plus 5mu minus 2mu
\thinmuskip= 0mu \medmuskip= 1mu plus 1mu minus 1mu \thickmuskip=2mu plus 3mu minus 1mu and use $\widetilde\sigma^2_t$ as the estimate of the latent conditional variance $\sigma_t^2$. See Table (ref) for $T_\text{in}$ and $T_\text{out}$ used in our datasets. See Martens:2002 and Fleming:2003 for the similar scaling estimator of the daily volatility.
Table (ref) summarizes the estimation results of fitting the SV, N-SV and SR-SV models to the five datasets. Table (ref) in Appendix (ref) shows the estimation results for the LMSV model. For the SR-SV model, we only show the results of the key parameters. Figure (ref) shows the posterior densities of the moving average weight parameter $\alpha$ of the SR-SV model in all simulation and real datasets. We draw some conclusions from Table (ref) and the listed figures.
First, the marginal likelihood estimates show that the SR-SV model fits the five index data better than the SV and N-SV models. The Bayes factors of the SR-SV model compared to the SV and N-SV models are more than $e^{2.3}$, which strongly support the SR-SV in all cases. There are no significant differences between the SV and N-SV models in terms of marginal likelihood estimates across the five panels.
Second, the evidence of the volatility effects rather than linearity, e.g. non-linearity and long-memory auto-dependence, in the volatility dynamics of the index datasets is clear as the posterior means of the non-linearity long-memory parameter $\beta_1$ of the SR-SV model are more than two standard deviations from zero in all cases. The estimation results of the LMSV in Table (ref) show that the posterior means of the fractional integration parameter $d$ are also more than two standard deviations from zero and close to $0.5$ in all cases, suggesting a strong evidence of the long-memory dependence in the volatility process of these five index datasets. The posterior means of the non-linear parameter $w_{z}$ with respect to the past log volatility $z_{t-1}$ are more than two standard deviations from zero, indicating the existence of the serial dependence rather than linearity that the past log volatility $z_{t-1}$ has on the current log volatility $z_t$, and that the SR-SV model is able to detect this serial dependence. The posterior density plots of parameter $\alpha$ in Figure (ref) suggest the existence of the long-memory auto-dependence in the volatility processes of the index datasets as the posterior densities of $\alpha$ are highly skewed for all cases and the posterior modes are close to 1, which is similar to the results in the SIM III data.
Third, it is worth noting that, in all cases, the persistence parameter $\phi$ in the SR-SV model is smaller than the persistence parameters in the SV and N-SV models as the parameter $w_{z}$ is significant and hence the linear effect that $z_{t-1}$ has on $z_t$ is reduced. As illustrated in the SIM I data, if the volatility process exhibits no volatility effects rather than a short-memory linear auto-dependence, the persistence parameters $\phi$ in the SV and SR-SV model are similar and the parameter $w_{z}$ is insignificant.
Using the posterior mean estimates in Table (ref), the {\it filtered} values of $z_t$ of the SR-SV, SV and N-SV models can be computed using the particle filter. We discuss how to obtain the filtered values of $z_t$ of the LMSV model in Appendix (ref). Figure (ref) plots the filtered log volatility of the SV and SR-SV model, together with the filtered values of the components $h_t$ and $\eta_t$ of the SR-SV model in all time steps, for the SPX data. Figure (ref), (ref), (ref) and (ref) in Appendix (ref) show the similar plots for the HSI, DAX, FCHI and TSX data, respectively. Figure (ref) shows that the component $h_t$, and hence $\eta_t$, of the SR-SV model is well responsive to changes in the volatility dynamics, e.g. being small during the low volatility periods and large in the high volatility periods of the SPX data. The SRU structure of the SR-SV model is able to capture these distinct behaviors of financial time series. We observe the similar behaviors of the $h_t$ and $\eta_t$ components for the other datasets as shown in the Figures (ref)-(ref).
Table (ref) provides summary statistics on the in-sample filtered volatilities and residuals $\widehat\epsilon_t^y$ of the SR-SV, SV, N-SV and LMSV models. We note that if the assumption about the normality of returns of the three models is justified then the residuals $\widehat\epsilon_t^y$ should have a zero skewness and a kurtosis of 3. Table (ref) shows that all the residual distributions produced from the four models are close to the standard normal distribution, but still slightly skewed and leptokurtic. The $p$-values of the Ljung-Box (LB) autocorrelation test of the residuals are high in all index datasets, suggesting that there is no evidence of autocorrelation in the residuals. We conjecture that extending the SR-SV model, for example, by using a Student's $t$ distribution instead of a Gaussian for the measurement shock $\epsilon^y_t$ and taking into account the leverage effect by correlating $\epsilon^y_t$ with the volatility shock $\epsilon_t^\eta$, is likely to lead to better diagnostics for the residuals. However, we do not consider these extensions here.
The marginal likelihood estimates in Table (ref) suggest that the SR-SV model fits the in-sample data of the five index datasets better than the SV and N-SV models. We now examine if this in-sample performance is consistent with the out-of-sample performance. Table (ref) provides summary statistics on the one-step-ahead forecasts of volatility and standardized residuals of the SR-SV, SV, N-SV and LMSV models. We note two conclusions from Table (ref).
First, the SR-SV model does not suffer from overfiting as often observed in neural network based volatility models Pagan:1990, Donaldson:1997, as the one-step-ahead forecast volatilities and the forecast residuals appear to be well behaved, compared to those from the more parsimonious SV and N-SV models. We emphasize that, as discussed in Section (ref), the use of noise-injecting regularization in the novel structure of the SR-SV model helps prevent it from the well-known overfitting problem.
Second, the means and standard deviations of one-step-ahead forecast volatility of the SR-SV model are smaller than those of the SV, N-SV and LMSV in all five index datasets. The SR-SV forecasts are generally more conservative in low volatility periods, in the sense that the forecast intervals often have a smaller band compared to the forecasts produced by the other models. The comparison of 99% one-step-ahead forecast intervals during the period Sep 2014 - May 2015 of the SPX data in Figure (ref) shows that the SR-SV model gives a safe buffer against abrupt changes in low volatility regions, e.g. Nov-Dec 2014, because it maintains a wider forecast band, while it does not produce overly large forecast intervals in high volatility regions, e.g. Oct 2014, Dec 2014 - Jan 2015. Therefore, the SR-SV model is less sensitive to the data values in the shorter time periods, and maintains a good trade-off between the information in recent observations and the information in the long-term memory. The SV, N-SV and LMSV models, compared to the SR-SV model, produce a smaller forecast volatility in low volatility regions and a higher volatility forecast in high volatility regions. The figure also shows that the SV and N-SV forecasts depend mainly on the return at the previous step, as the persistence parameters $\phi$ in the SV and N-SV models are larger than the persistence parameter of the SR-SV model. The SR-SV intervals are closer to the intervals made by the realized variance and hence seem to track the out-of-sample returns better than the SV, N-SV and LMSV models.
Table (ref) shows the out-of-sample performances of the four models on the SPX data, with the Monte Carlo standard errors in brackets. Tables (ref), (ref), (ref) and (ref) in Appendix (ref) show the results for the other four datasets. Each table provides predictive scores separately in four panels, corresponding to the four realized measures BV, MedRV, RKV and RV, as discussed Section (ref). In each panel, we count the number of times a particular model has lowest (best) predictive scores and list these numbers in the last column; the model with the highest count is preferred. We note that the PPS predictive score is independent of the realized measures of volatility, and that the PPS predictive score is not available for the LMSV model.
As shown in these tables, the SR-SV model consistently has the best out-of-sample performance in all the five index datasets. The superior predictive performance of the SR-SV model is consistent with its in-sample performance discussed earlier and provides further evidence to support the conclusion that the SR-SV model does not overfit the five index datasets. We note that the forecast performances of the SV and N-SV models are mixed with no model consistently outperforming the other across the predictive scores and the datasets. The LMSV model consistently makes the least accurate forecasts in all cases .
This paper proposes a statistical recurrent stochastic volatility (SR-SV) model, by combining the statistical recurrent unit architecture from Machine Learning and the stochastic volatility model from Financial Econometrics. These two techniques are combined in a principled and non trivial way to form a new approach that is, as carefully illustrated through the extensive simulation and empirical studies, highly efficient for volatility modelling and forecasting. It is easy to carry Bayesian inference in the SR-SV model using standard Bayesian computation methods, such as the Density Tempered Sequential Monte Carlo method as used in this paper. The simulation and empirical studies suggest that the SR-SV model is able to capture various volatility effects overlooked by the SV benchmark models, and is able to produce highly accurate forecast volatilities.
Extending the SR-SV model by incorporating features such as the leverage effect is an interesting research question. Another interesting research question is extending the present SR-SV model to multivariate financial time series. We conjecture that the RNN architectures will be even more powerful in this case as they can naturally capture the interaction between the inputs. This research is in progress.