My stuff
A Neural Stochastic Volatility Model

In this paper, we show that the recent integration of statistical models with deep recurrent neural networks provides a new way of formulating volatility (the degree of variation of time series) models that have been widely used in time series analysis and prediction in finance. The model comprises a pair of complementary stochastic recurrent neural networks: the generative network models the joint distribution of the stochastic volatility process; the inference network approximates the conditional distribution of the latent variables given the observables. Our focus here is on the formulation of temporal dynamics of volatility over time under a stochastic recurrent neural network framework. Experiments on real-world stock price datasets demonstrate that the proposed model generates a better volatility estimation and prediction that outperforms mainstream methods, e.g., deterministic models such as GARCH and its variants, and stochastic models namely the MCMC-based model stochvol as well as the Gaussian process volatility model GPVol, on average negative log-likelihood.

The volatility of the price movements reflects the ubiquitous uncertainty within financial markets. It is critical that the level of risk (aka, the degree of variation), indicated by volatility, is taken into consideration before investment decisions are made and portfolio are optimised (Hull 2006); volatility is substantially a key variable in the pricing of derivative securities. Hence, estimating and forecasting volatility is of great importance in branches of financial studies, including investment, risk management, security valuation and monetary policy making (Poon and Granger 2003).

Volatility is measured typically by employing the standard deviation of price change in a fixed time interval, such as a day, a month or a year. The higher the volatility is, the riskier the asset should be. One of the primary challenges in designing volatility models is to identify the existence of latent stochastic processes and to characterise the underlying dependences or interactions between variables within a certain time span. A classic approach has been to handcraft the characteristic features of volatility models by imposing assumptions and constraints, given prior knowledge and observations. Notable examples include autoregressive conditional heteroscedasticity (ARCH) model (Engle 1982) and the extension, generalised ARCH (GARCH) (Bollerslev 1986), which makes use of autoregression to capture the properties of time-varying volatility within many time series. As an alternative to the GARCH model family, the class of stochastic volatility (SV) models specify the variance to follow some latent stochastic process (Hull and White 1987). Heston (Heston 1993) proposed a continuous-time model with the volatility following an Ornstein-Uhlenbeck process and derived a closed-form solution for options pricing. Since the temporal discretisation of continuous-time dynamics sometimes leads to a deviation from the original trajectory of system, those continuous-time models are seldom applied in forecasting. For practical purposes of forecasting, the canonical model (Jacquier, Polson, and Rossi 2002; Kim, Shephard, and Chib 1998) formulated in a discrete-time fashion for regularly spaced data such as daily prices of stocks is of great interest. While theoretically sound, those approaches require strong assumptions which might involve detailed insight of the target sequences and are difficult to determine without a thorough inspection.

In this paper, we take a fully data driven approach and determine the configurations with as few exogenous input as possible, or even purely from the historical data. We propose a neural network re-formulation of stochastic volatility by leveraging stochastic models and recurrent neural networks (RNNs). In inspired by the work from Chung et al. (Chung et al. 2015) and Fraccaro et al. (Fraccaro et al. 2016), the proposed model is rooted in variational inference and equipped with the latest advances of stochastic neural networks. The model inherits the fundamentals of SV model and provides a general framework for volatility modelling; it extends previous sequential frameworks with autoregressive and bidirectional architecture and provide with a more systematic and volatility-specific formulation on stochastic volatility modelling for financial time series. We presume that the latent variables follow a Gaussian autoregressive process, which is then utilised to model the variance process. Our neural network formulation is essentially a general framework for volatility modelling, which covers two major classes of volatility models in financial study as the special cases with specific weights and activations on neurons.

Experiments with real-world stock price datasets are performed. The result shows that the proposed model produces more accurate estimation and prediction, outperforming various widely-used deterministic models in the GARCH family and several recently proposed stochastic models on average negative log-likelihood; the high flexibility and rich expressive power are validated.

A notable framework for volatility is autoregressive conditional heteroscedasticity (ARCH) model (Engle 1982): it can accurately identify the characteristics of time-varying volatility within many types of time series. Inspired by ARCH model, a large body of diverse work based on stochastic process for volatility modelling has emerged (Bollerslev, Engle, and Nelson 1994). Bollerslev (Bollerslev 1986) generalised ARCH model to the generalised autoregressive conditional heteroscedasticity (GARCH) model in a manner analogous to the extension from autoregressive (AR) model to autoregressive moving average (ARMA) model by introducing the past conditional variances in the current conditional variance estimation. Engle and Kroner (Engle and Kroner 1995) presented theoretical results on the formulation and estimation of multivariate GARCH model within simultaneous equations systems. The extension to multivariate model allows the covariance to present and depend on the historical information, which are particularly useful in multivariate financial models. An alternative to the conditionally deterministic GARCH model family is the class of stochastic volatility (SV) models, which first appeared in the theoretical finance literature on option pricing (Hull and White 1987). The SV models specify the variance to follow some latent stochastic process such that the current volatility is no longer a deterministic function even if the historical information is provided. As an example, Heston’s model (Heston 1993) characterises the variance process as a Cox-Ingersoll-Ross process driven by a latent Wiener process. While theoretically sound, those approaches require strong assumptions which might involve complex probability distributions and non-linear dynamics that drive the process. Nevertheless, empirical evidences have confirmed that volatility models provide accurate prediction (Ander- sen and Bollerslev 1998) and models such as ARCH and its descendants/variants have become indispensable tools in asset pricing and risk evaluation. Notably, several models have been recently proposed for practical forecasting tasks: Kastner et al. (Kastner and Fr¨uhwirth-Schnatter 2014) implemented the MCMC-based framework stochvol where the ancillarity-sufficiency interweaving strategy (ASIS) is applied for boosting MCMC estimation of stochastic volatility; Wu et al. (Wu, Hern´andez-Lobato, and Ghahramani 2014) designed the GP-Vol, a non-parametric model which utilises Gaussian processes to characterise the dynamics and jointly learns the process and hidden states via online inference algorithm. Despite the fact that it provides us with a practical approach towards stochastic volatility forecasting, both models require a relatively large volume of samples to ensure the accuracy, which involves very expensive sampling routine at each time step. Another drawback is that those models are incapable to handle the forecasting task for multivariate time series.

On the other hand, deep learning (LeCun, Bengio, and Hinton 2015; Schmidhuber 2015) that utilises nonlinear structures known as deep neural networks, powers various applications. It has triumph over pattern recognition challenges, such as image recognition (Krizhevsky, Sutskever, and Hinton 2012), speech recognition (Chorowski et al. 2015), machine translation (Bahdanau, Cho, and Bengio 2014) to name a few.

Time-dependent neural networks models include RNNs with neuron structures such as long short-term memory (LSTM) (Hochreiter and Schmidhuber 1997), bidirectional RNN (BRNN) (Schuster and Paliwal 1997), gated recurrent unit (GRU) (Cho et al. 2014) and attention mechanism (Xu et al. 2015). Recent results show that RNNs excel for sequence modelling and generation in various applications (van den Oord, Kalchbrenner, and Kavukcuoglu 2016; Cho et al. 2014; Xu et al. 2015). However, despite its capability as non-linear universal approximator, one of the drawbacks of neural networks is its deterministic nature. Adding latent variables and their processes into neural networks would easily make the posteriori computationally intractable. Recent work shows that efficient inference can be found by variational inference when hidden continuous variables are embedded into the neural networks structure (Kingma and Welling 2013; Rezende, Mohamed, and Wier- stra 2014). Some early work has started to explore the use of variational inference to make RNNs stochastic: Chung et al. (Chung et al. 2015) defined a sequential framework with complex interacting dynamics of coupling observable and latent variables whereas Fraccaro et al. (Fraccaro et al. 2016) utilised heterogeneous backward propagating layers in inference network according to its Markovian properties.

In this paper, we apply the stochastic neural networks to solve the volatility modelling problem. In other words, we model the dynamics and stochastic nature of the degree of variation, not only the mean itself. Our neural network treatment of volatility modelling is a general one and existing volatility models (e.g., the Heston and GARCH models) are special cases in our formulation.

Volatility models characterise the dynamics of volatility processes, and help estimate and forecast the fluctuation within time series. As it is often the case that one seeks for prediction on quantity of interest with a collection of historical information at hand, we presume the conditional variance to have dependency – either deterministic or stochastic – on history, which results in two categories of volatility models.

Deterministic Volatility Models: the GARCH Model Family

The GARCH model family comprises various linear models that formulate the conditional variance at present as a linear function of observations and variances from the past. Bollerslev’s extension (Bollerslev 1986) of Engle’s primitive ARCH model (Engle 1982), referred as generalised ARCH (GARCH) model, is one of the most well-studied

and widely-used volatility models:


where Eq. (2) represents the assumption that the observation xtfollows from the Gaussian distribution with mean 0 and variance  σ2t; the (conditional) variance  σ2tis fully de- termined by a linear function (Eq. (1)) of previous observations  {x<t}and variances  {σ2<t}. Note that if q = 0 in Eq. (1), GARCH model degenerates to ARCH model.

Various variants have been proposed ever since. Glosten, Jagannathan and Runkle (Glosten, Jagannathan, and Run- kle 1993) extended GARCH model with additional terms to account for asymmetries in the volatility and proposed GJR-GARCH model; Zakoian (Zakoian 1994) replaced the quadratic operators with absolute values, leading to threshold ARCH/GARCH (TARCH) models. The general functional form is formulated as


where I[xt−k < 0]is the indicator function: I[xt−k < 0] = 1if xt−k < 0, and 0 otherwise, which allows for asymmetric reactions of volatility in terms of the sign of previous  {x<t}.

Various variants of GARCH model can be expressed by assigning values to parameters p, o, q, d in Eq. (3):


Another fruitful specification shall be Nelson’s exponential GARCH (EGARCH) model (Nelson 1991), which instead formulates the dependencies in log-variance  log(σ2t ):


where  g(xt)(Eq. (5)) accommodates the asymmetric relation between observations and volatility changes. If setting q  ≡ 0in Eq. (4), EGARCH(p, q) then degenerates to EARCH(p).

Stochastic Volatility Models An alternative to the (conditionally) deterministic volatility models is the class of stochastic volatility (SV) models. First introduced in the theoretical finance literature, earliest SV models such as Hull and White’s (Hull and White 1987) as well as Heston model (Heston 1993) are formulated by

stochastic differential equations in a continuous-time fashion for analysis convenience. In particular, Heston model instantiates a continuous-time stochastic volatility model for univariate processes:


where x(t) = log s(t) is the logarithm of stock price stat time t,  wx(t)and  wσ(t)represent two correlated Wiener processes and the correlation between  dwx(t)and  dwσ(t)is expressed as  E[dwx(t) · dwσ(t)] = ρ dt.

For practical use, empirical versions of the SV model, typically formulated in a discrete-time fashion, are of great interest. The canonical model (Jacquier, Polson, and Rossi 2002; Kim, Shephard, and Chib 1998) for regularly spaced data is formulated as


Equation (8) indicates that the (conditional) log-variance log(σ2t )depends on not only the historical log-variances {log(σ2t )}but a latent stochastic process  {zt}. The latent process  {zt}is, according to Eq. (9), white noise process with i.i.d. Gaussian variables.

Notably, the volatility  σ2tis no longer conditionally deter- ministic (i.e. deterministic given the complete history  {σ2<t}) but to some extent stochastic in the setting of SV models: Heston model involves two correlated continuous-time Wiener processes while the canonical model is driven by a discrete-time Gaussian white-noise process.

Volatility Models in a General Form

Hereafter we denote the sequence of observations as  {xt}and the latent stochastic process as  {zt}. As seen in previous sections, the dynamics of the volatility process  {σ2t }can be abstracted in the form of


The latter equality holds when we recurrently substitute σ2τwith f  (σ2<τ, x<τ, z≤τ)for all  τ <t. For models within the GARCH family, we discard z≤tin  Σ x(x<t, z≤t)(Eq. (10)); on the other hand, for the primitive SV model, x<tis ignored instead. We can loosen the constraint that xtis zero-mean to a time-varying mean  µx(x<t, z≤t)for more flexibility.

Recall that the latent stochastic process  {zt}(Eq. (9)) in the SV model is by definition an i.i.d. Gaussian white noise process. We may extend this process to one with an inherent autoregressive dynamics and more flexibility that the mean µz(z<t)and variance  Σz(z<t)are functions of autoregressive structure on historical values. Hence, the generalised model can be formulated in the following framework:


where we have presumed that both the observation xtand the latent variable ztare normally distributed. Note that the autoregressive process degenerates to i.i.d. white noise process when  µz(z<t) ≡ 0and  Σz(z<t) ≡ σ2z. It should be emphasised that the purpose of reinforcing an autoregressive structure (11) of the latent variable ztis that we believe such formulation fits better to real scenarios from financial aspect compared with the i.i.d. convention: the price fluctua-tion of a certain stock is the consequence of not only its own history but also the influence from the environment, e.g. its competitors, up/downstream industries, relevant companies in the market, etc. Such external influence is ever-changing and may preserve memory and hence hard to characterise if restricted to i.i.d. noise. The latent variable ztwith an autoregressive structure provides a possibility of decoupling the internal influential factors from the external ones, which we believe is the essence of introducing zt.

In this section, we establish the neural stochastic volatility model (NSVM) for volatility estimation and prediction.

Generating Observable Sequence Recall that the observable variable xt(Eq. (12)) and the latent variable zt(Eq. (11)) are described by autoregressive models (as xtalso involves an exogenous input z≤t). Let pΦ(xt|x<t, z≤t)and pΦ(zt|z<t)denote the probability distributions of xtand ztat time t. The factorisation on the joint distributions of sequences  {xt}and  {zt}applies as follow:


where X  = {x1:T }and Z  = {z1:T }represents the sequences of observable and latent variables, respectively, whereas  Φstands for the collection of parameters of generative model. The unconditional generative model is defined as the joint distribution w.r.t. the latent variable Z and observable X:


It can be observed that the mean and variance are conditionally deterministic: given the historical information  {z<t}, the current mean  µzt = µzΦ(z<t)and variance Σzt = ΣzΦ(z<t)of ztis obtained and hence the distribu- tion  N(zt; µzt, Σzt )of ztis specified; after sampling ztfrom the specified distribution, we incorporate  {x<t}and calculate the current mean  µxt = µxΦ(x<t, z≤t)and variance Σ xt = Σ xΦ(x<t, z≤t)of xtand determine its distribution N(xt; µxt , Σ xt )of xt. It is natural and convenient to present such a procedure in a recurrent fashion because of its autoregressive nature. Since RNNs can essentially approximate arbitrary function of recurrent form, the means and variances, which may be driven by complex non-linear dynamics, can be efficiently computed using RNNs.

The unconditional generative model consists of two pairs of RNN and multi-layer perceptron (MLP), namely RNNzG/MLPzGfor the latent variable and  RNNxG/MLPxGfor the observable. We stack those two RNN/MLP pairs together according to the causal dependency between variables. The unconditional generative model is implemented as the generative network abstracted as follows:


where hztand hxtdenote the hidden states of the correspond- ing RNNs. The MLPs map the hidden states of RNNs into the means and deviations of variables of interest. The collection of parameters  Φis comprised of the weights of RNNs and MLPs. NSVM relaxes the conventional constraint that the latent variable ztis N(0, 1) in a way that ztis no longer i.i.d noise but a time-varying signal from external process with self-evolving nature. As discussed above, this relaxation will benefit the effectiveness in real scenarios.

One should notice that when the latent variable ztis obtained, e.g. by inference (see details in the next subsection), the conditional distribution pΦ(X|Z)(Eq. (14)) will be involved in generating the observable xtinstead of the joint distribution pΦ(X, Z)(Eq. (15)). This is essentially the scenario of predicting future values of the observable variable given its history. We will use the term “generative model” and will not discriminate the unconditional generative model or the conditional one as it can be inferred in context.

Inferencing the Latent Process

As the generative model involves the latent variable zt, of which the true values are inaccessible even we have observed xt, the marginal distribution pΦ(X)becomes the key that bridges the model and the data. However, the calcula- tion of pΦ(X)itself or its complement, the posterior distribution pΦ(Z|X), is often intractable as complex integrals are involved. We are unable to learn the parameters by differentiating the marginal log-likelihood  log pΦ(X)or to infer the latent variables through the true posterior. Therefore, we consider instead a restricted family of tractable distributions qΨ(Z|X), referred to as the approximate posterior family, as approximations to the true posterior pΦ(Z|X)such that the family is sufficiently rich and of high capacity to provide good approximations.

It is straightforward to verify that given a sequence of observations X  = {x1:T }, for any  1 ≤ t ≤T, ztis dependent on the entire observation sequences. Hence, we define the inference model with the spirit of mean-field approximation where the approximate posterior is Gaussian and the follow-

ing factorisation applies:


where  ˜µzΨ(zt−1, x1:T)and  ˜ΣzΨ(zt−1, x1:T)are functions of the given observation sequence  {x1:T }, representing the approximated mean and variance of the latent variable zt; Ψdenotes the collection of parameters of inference model.

The neural network implementation of the model, referred to as the inference network, is designed to equip a cascaded architecture with an autoregressive RNN and a bidirectional RNN, where the bidirectional RNN incorporates both the forward and backward dependencies on the entire observations whereas the autoregressive RNN models the temporal dependencies on the latent variables:


where ˜h→tand ˜h←trepresent the hidden states of the forward and backward directions of the bidirectional RNN. The autoregressive RNN with hidden state ˜hzttakes the joint state [˜h→t , ˜h←t ]of the bidirectional RNN and the previous value of zt−1as input. The inference mean  ˜µztand variance  ˜Σztis computed by an MLP from the hidden state ˜hztof the autore- gressive RNN. We use the subscript I instead of G to distinguish the architecture used in inference model in contrast to that of the generative model. It should be emphasised that the inference network will collaborates with the generative network on conditional generating procedure.


Forecasting Observations in Future

For a volatility model to be practically applicable in forecasting, the generating procedure conditioning on the history is of essential interest. We start with 1-step-ahead prediction, which serves as building block of multi-step forecasting.

Given the historical observations  {x1:T }up to time stepT, 1-step-ahead prediction of either  Σ xT+1or xT+1is fully

depicted by the conditional predictive distribution:


where the distributions on the right-hand side refer to those in the generative model with the generative parameters  Φomitted. As the true posterior p(z1:T |x1:T)involved in Eq. (28) is intractable, the exact evaluation of conditional predictive distribution p(xT+1|x1:T)is difficult.

A straightforward solution is that we substitute the true posterior p(z1:T |x1:T)with the approximation q(z1:T |x1:T)(see Eq. (22)) and leverage q(z1:T |x1:T)to inference S sample paths  {z⟨1:S⟩1:T }of the latent variables according to the historical observations  {x1:T }. The approximate posterior from a well-trained model is presumed to be a good approximation to the truth; hence the sample paths shall be mimics of the true but unobservable path. We then extend the sample paths one step further from T to T + 1 using the autoregressive generative distribution p(zT+1|z1:T)(see Eq. (13)). The conditional predictive distribution is thus approximated as


which is essentially a mixture of S Gaussians. In the case of multi-step forecasting, a common solution in practice is to perform a recursive 1-step-ahead forecasting routine with model updated as new observation comes in; the very same procedure can be applied except that more sample paths should be evaluated due to the accumulation of uncertainty. Algorithm 1 gives the detailed rolling scheme.

In this section, we present the experiment on real-world stock price time series to validate the effectiveness and to evaluate the performance of the prosed model.

Dataset and Pre-processing The raw dataset comprises 162 univariate time series of the daily closing stock price, chosen from China’s A-shares and collected from 3 institutions. The choice is made by selecting those with earlier listing date of trading (from 2006 or earlier) and fewer suspension days (at most 50 suspension days within the entire period of observation), such that the undesired noises introduced by insufficient observation or missing values – highly influential on the performance but essentially irrelevant to the purpose of volatility modelling – can be reduced to the minimum. The raw price series is cleaned by aligning and removing abnormalities: we manually aligned the mismatched part and interpolated the missing value by stochastic regression imputation (Little and Ru- bin 2014) where the imputed value is drawn from a Gaussian distribution with mean and variance calculated by regression on the empirical value within a short interval of 20 recent days. The series is then transformed from actual prices stinto log-returns xt = log(st/st−1)and normalised. Moreover, we combinatorically choose a predefined number d out

of 162 univariate log-return series and aggregate the selected series at each time step to form a d-dimensional multivariate time series, the choice of d is in accordance with the rank of correlation, e.g. d = 6 in our experiments. Theoretically, it leads to a much larger volume of data as �1626� > 2 × 1010. Specifically, the actual dataset for training and evaluation comprises a collection of 2000 series of d-dimensional normalised log-return vectors of length  2570 (∼7 years) with no missing values. We divide the whole dataset into two subsets for training and testing along the time axis: the first 2000 time steps of each series have been used as training samples whereas the rest 570 steps of each series as the test samples.

Baselines We select several deterministic volatility models from the GARCH family as baselines: 1. Quadratic models

ARCH(1); GARCH(1,1); GJR-GARCH(1,1,1);

2. Absolute value models

AVARCH(1); AVGARCH(1,1); TARCH(1,1,1);

3. Exponential models.

EARCH(1); EGARCH(1,1);

Moreover, two stochastic volatility models are compared:

1. MCMC volatility model: stochvol;

2. Gaussian process volatility model GP-Vol.

For the listed models, we retrieve the authors’ implementations or tools: stochvol1, GP-Vol2 (the hyperparameters are chosen as suggested in (Wu, Hern´andez-Lobato, and Ghahramani 2014)) and implement the models, such as GARCH, EGARCH, GJR-GARCH, etc., based on several widely-used packages345 for time series analysis. All baselines are evaluated in terms of the negative log-likelihood on the test samples, where 1-step-ahead forecasting is carried out in a recursive fashion similar to Algorithm 1.

Model Implementation In our experiments, we predefine the dimensions of observable variables to be  dim xt = 6and the latent variables dim zt = 4. Note that the dimension of the latent variable is smaller than that of the observable, which allows us to extract a compact representation. The NSVM implementation in our experiments is composed of two neural networks, namely the generative network (see Eq. (16)-(21)) and inference network (see Eq. (23)-(27)). Each RNN module contains one hidden layer of size 10 with GRU cells; MLP modules are 2-layered fully-connected feedforward networks, where the hidden layer is also of size 10 whereas the output layer splits into two equal-sized sublayers with different activation functions: one applies exponential function to ensure the non-negativity for variance while the

other uses linear function to calculate mean estimates. Thus MLPzI’s output layer is of size 4 + 4 for  { ˜µz, ˜Σz}whereas the size of  MLPxG’s output layer is 6 + 6 for  {µx, Σ x}. Dur- ing the training phase, the inference network is connected with the conditional generative network (see, Eq. (16)-(18)) to establish a bottleneck structure, the latent variable ztinferred by variational inference (Kingma and Welling 2013; Rezende, Mohamed, and Wierstra 2014) follows a Gaussian approximate posterior; the size of sample paths is set to S = 100. The parameters of both networks are jointly learned, including those for the prior. We introduce Dropout (Srivastava et al. 2014) into each RNN modules and impose L2-norm on the weights of MLP modules as regularistion to prevent overshooting; Adam optimiser (Kingma and Ba 2014) is exploited for fast convergence; exponential learning rate decay is adopted to anneal the variations of convergence as time goes. Two covariance configurations are adopted: 1. we stick with diagonal covariance matrices con-figurations; 2. we start with diagonal covariance and then apply rank-1 perturbation (Rezende, Mohamed, and Wierstra 2014) during fine-tuning until training is finished. The recursive 1-step-ahead forecasting routine illustrated as Algorithm 1 is applied in the experiment for both training and test phase: during the training phase, a single NSVM is trained, at each time step, on the entire training samples to learn a holistic dynamics, where the latent shall reflect the evolution of environment; in the test phase, on the other hand, the model is optionally retrained, at every 20 time steps, on each particular input series of the test samples to keep track on the specific trend of that series. In other words, the trained NSVM predicts 20 consecutive steps before it is retrained using all historical time steps of the input series at present. Correspondingly, all baselines are trained and tested at every time step of each univariate series using standard calibration procedures. The negative log-likelihood on test samples has been collected for performance assessment. We train the model on a single-GPU (Titan X Pascal) server for roughly two hours before it converges to a certain degree of accuracy on the training samples. Empirically, the training phase can be processed on CPU in reasonable time, as the complexity of the model as well as the size of parameters is moderate.

Result and Discussion

The performance of NSVM and baselines is listed for comparison in Table 1: the performance on the first 10 individual stocks (chosen in alphacetical order but anonymised here) and the average score on all 162 stocks are reported in terms of negative log-likelihood (NLL) measure. The result shows that NSVM has achieved higher accuracy over the baselines on the task of volatility modelling and forecasting on NLL, which validates the high flexibility and rich expressive power of NSVM for volatility modelling and forecasting. In particular, NSVM with rank-1 perturbation (referred to as NSVM-corr in Table 1) beats all other models in terms of NLL, while NSVM with diagonal covariance matrix (i.e. NSVM-diag) outperforms GARCH(1,1) on 142 out of 162 stocks. Although the improvement comes at the cost of longer training time before convergence, it can be mitigated by applying parallel computing techniques as well as more advanced network architecture or training methods.

Table 1: The performance of the proposed model and the baselines in terms of negative log-likelihood (NLL) evaluated on the test samples of real-world stock price time series: each row from 1 to 10 lists the average NLL for a specific individual stock; the last row summarises the average NLL of the entire test samples of all 162 stocks.


Figure 1: Case studies of volatility forecasting.

Apart from the higher accuracy NSVM obtained, it provides us with a rather general framework to generalise univariate time series models of any specific functional form to the corresponding multivariate cases by extending network dimensions and manipulating the covariance matrices. A case study on real-world financial datasets is illustrated in Fig. 1.

NSVM shows higher sensibility on drastic changes and better stability on moderate fluctuations: the response of NSVM in Fig. 1a is more stable in t  ∈ [1600, 2250], the period of moderate price fluctuation; while for drastic price change at t = 2250, the model responds with a sharper


spike compared with the quadratic GARCH model. Furthermore, NSVM demonstrates the inherent non-linearity in both Fig. 1a and 1b: at each time step within t ∈[1000, 2000], the model quickly adapts to the current fluctu-ation level whereas GARCH suffers from a relatively slower decay from the previous influences. The cyan vertical line at t = 2000 splits the training samples and test samples. We show only one instance within our dataset due to the limitation of pages, the performance of other instances are similar.

In this paper, we proposed a new volatility model, referred to as NSVM, for volatility estimation and prediction. We integrated statistical models with deep neural networks, leveraged the characteristics of each model, organised the dependences between random variables in the form of graphical models, implemented the mappings among variables and parameters through RNNs and MLPs, and finally established a powerful stochastic recurrent model with universal approximation capability. The proposed architecture comprises a pair of complementary stochastic neural networks: the generative network and inference network. The former models the joint distribution of the stochastic volatility process with both observable and latent variables of interest; the latter provides with the approximate posterior i.e. an analytical approximation to the (intractable) conditional distribution of the latent variables given the observable ones. The parameters (and consequently the underlying distributions) are learned (and inferred) via variational inference, which maximises the lower bound for the marginal log-likelihood of the observable variables. NSVM has presented higher accuracy on the task of volatility modelling and forecasting on real-world financial datasets, compared with various widely-used models, such as GARCH, EGARCH, GJRGARCH, TARCH in the GARCH family, MCMC volatility model stochvol as well as Gaussian process volatility model GP-Vol. Future work on NSVM would be to investigate the modelling of time series with non-Gaussian residual distributions, in particular the heavy-tailed distributions e.g. LogNormal log Nand Student’s t-distribution.

[Andersen and Bollerslev 1998] Andersen, T. G., and Bollerslev, T. 1998. Answering the skeptics: Yes, standard volatility models do provide accurate forecasts. International economic review 885– 905.

[Bahdanau, Cho, and Bengio 2014] Bahdanau, D.; Cho, K.; and Bengio, Y. 2014. Neural machine translation by jointly learning to align and translate. CoRR abs/1409.0473.

[Bollerslev, Engle, and Nelson 1994] Bollerslev, T.; Engle, R. F.; and Nelson, D. B. 1994. Arch models. Handbook of econometrics 4:2959–3038.

[Bollerslev 1986] Bollerslev, T. 1986. Generalized autoregressive conditional heteroskedasticity. Journal of econometrics 31(3):307– 327.

[Cho et al. 2014] Cho, K.; van Merrienboer, B.; G¨ulc¸ehre, C¸ .; Bah- danau, D.; Bougares, F.; Schwenk, H.; and Bengio, Y. 2014. Learning phrase representations using RNN encoder-decoder for statistical machine translation. In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing, 1724–1734.

[Chorowski et al. 2015] Chorowski, J.; Bahdanau, D.; Serdyuk, D.; Cho, K.; and Bengio, Y. 2015. Attention-based models for speech recognition. In Advances in Neural Information Processing Systems 28, 577–585.

[Chung et al. 2015] Chung, J.; Kastner, K.; Dinh, L.; Goel, K.; Courville, A. C.; and Bengio, Y. 2015. A recurrent latent variable model for sequential data. In Advances in Neural Information Processing Systems 28, 2980–2988.

[Engle and Kroner 1995] Engle, R. F., and Kroner, K. F. 1995. Multivariate simultaneous generalized arch. Econometric theory 11(01):122–150.

[Engle 1982] Engle, R. F. 1982. Autoregressive conditional het- eroscedasticity with estimates of the variance of united kingdom inflation. Econometrica: Journal of the Econometric Society 987– 1007.

[Fraccaro et al. 2016] Fraccaro, M.; Sønderby, S. K.; Paquet, U.; and Winther, O. 2016. Sequential neural models with stochastic layers. In Advances in Neural Information Processing Systems 29, 2199–2207.

[Glosten, Jagannathan, and Runkle 1993] Glosten, L. R.; Jagan- nathan, R.; and Runkle, D. E. 1993. On the relation between the expected value and the volatility of the nominal excess return on stocks. The journal of finance 48(5):1779–1801.

[Heston 1993] Heston, S. L. 1993. A closed-form solution for op- tions with stochastic volatility with applications to bond and currency options. Review of financial studies 6(2):327–343.

[Hochreiter and Schmidhuber 1997] Hochreiter, S., and Schmidhu- ber, J. 1997. Long short-term memory. Neural Computation 9(8):1735–1780.

[Hull and White 1987] Hull, J., and White, A. 1987. The pricing of options on assets with stochastic volatilities. The journal of finance 42(2):281–300.

[Hull 2006] Hull, J. C. 2006. Options, futures, and other derivatives. Pearson Education India.

[Jacquier, Polson, and Rossi 2002] Jacquier, E.; Polson, N. G.; and Rossi, P. E. 2002. Bayesian analysis of stochastic volatility models. Journal of Business & Economic Statistics 20(1):69–87.

[Kastner and Fr¨uhwirth-Schnatter 2014] Kastner, G., and Fr¨uhwirth-Schnatter, S. 2014. Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models. Computational Statistics & Data Analysis 76:408–423.

[Kim, Shephard, and Chib 1998] Kim, S.; Shephard, N.; and Chib, S. 1998. Stochastic volatility: likelihood inference and comparison with arch models. The review of economic studies 65(3):361–393.

[Kingma and Ba 2014] Kingma, D. P., and Ba, J. 2014. Adam: A method for stochastic optimization. CoRR abs/1412.6980.

[Kingma and Welling 2013] Kingma, D. P., and Welling, M. 2013. Auto-encoding variational bayes. CoRR abs/1312.6114.

[Krizhevsky, Sutskever, and Hinton 2012] Krizhevsky, A.; Sutskever, I.; and Hinton, G. E. 2012. Imagenet classifica-tion with deep convolutional neural networks. In Advances in Neural Information Processing Systems 25, 1106–1114.

[LeCun, Bengio, and Hinton 2015] LeCun, Y.; Bengio, Y.; and Hin- ton, G. E. 2015. Deep learning. Nature 521(7553):436–444.

[Little and Rubin 2014] Little, R. J., and Rubin, D. B. 2014. Statistical analysis with missing data. John Wiley & Sons.

[Nelson 1991] Nelson, D. B. 1991. Conditional heteroskedasticity in asset returns: A new approach. Econometrica: Journal of the Econometric Society 347–370.

[Poon and Granger 2003] Poon, S.-H., and Granger, C. W. 2003. Forecasting volatility in financial markets: A review. Journal of economic literature 41(2):478–539.

[Rezende, Mohamed, and Wierstra 2014] Rezende, D. J.; Mohamed, S.; and Wierstra, D. 2014. Stochastic backpropagation and approximate inference in deep generative models. In Proceedings of the 31th International Conference on Machine Learning, 1278–1286.

[Schmidhuber 2015] Schmidhuber, J. 2015. Deep learning in neural networks: An overview. Neural Networks 61:85–117.

[Schuster and Paliwal 1997] Schuster, M., and Paliwal, K. K. 1997. Bidirectional recurrent neural networks. IEEE Trans. Signal Processing 45(11):2673–2681.

[Srivastava et al. 2014] Srivastava, N.; Hinton, G. E.; Krizhevsky, A.; Sutskever, I.; and Salakhutdinov, R. 2014. Dropout: a simple way to prevent neural networks from overfitting. Journal of Machine Learning Research 15(1):1929–1958.

[van den Oord, Kalchbrenner, and Kavukcuoglu 2016] van den Oord, A.; Kalchbrenner, N.; and Kavukcuoglu, K. 2016. Pixel recurrent neural networks. In Proceedings of the 33nd International Conference on Machine Learning, 1747–1756.

[Wu, Hern´andez-Lobato, and Ghahramani 2014] Wu, Y.; Hern´andez-Lobato, J. M.; and Ghahramani, Z. 2014. Gaussian process volatility model. In Advances in Neural Information Processing Systems 27, 1044–1052.

[Xu et al. 2015] Xu, K.; Ba, J.; Kiros, R.; Cho, K.; Courville, A. C.; Salakhutdinov, R.; Zemel, R. S.; and Bengio, Y. 2015. Show, attend and tell: Neural image caption generation with visual attention. In Proceedings of the 32nd International Conference on Machine Learning, 2048–2057.

[Zakoian 1994] Zakoian, J.-M. 1994. Threshold heteroskedastic models. Journal of Economic Dynamics and control 18(5):931– 955.

Given the observations X  = {x1:T }, we target at maximising the marginal log-likelihood pΦ(X)w.r.t.  Φ, where the actual posterior pΦ(Z|X)is involved. Because of the intractability of pΦ(Z|X), the exact inference is not applicable; we have to seek for approximate solutions instead.

We factorise the marginal log-likelihood  log pΦ(X)as


Note that we have introduced a tractable,  Ψ-parameterised distribution qΨ(Z|X)from a flexible family of distributions to approximate the actual posterior pΦ(Z|X). The evidence lower bound (ELBO)  L[q; X, Φ, Ψ]in Eq. (30) is essentially a functional w.r.t. q, conditioning on the observations X and parameterised by the parameter sets  Φ, Ψof both generative and inference models. Theoretically, ELBO ensures a lower bound on the marginal log-likelihood, and can be maximised via gradient-based optimisers.

It is usually the case that Eq. (30) lacks a closed-form expression. We have to estimate the ELBO using Monte Carlo integration on the latent variable zt. Provided S sample paths drawn by the inference model defined in Eq. (22), the estimator of ELBO can be calculated as the path average:


where  {z⟨1:S⟩1:T }denotes the collection of S sample paths.

By assuming the latent variable ztbeing Gaussian, we can readily apply the reparameterisation technique (Kingma and Welling 2013) to ztto form an unbiased gradient estimator:


where  ϵ ⟨s⟩t ∼ N(0, Iz)is the standard Gaussian variable.


Figure 2: Generalised stochastic volatility model.

The reparameterisation extracts the randomness out of the latent variable z⟨s⟩tvia  ϵ ⟨s⟩t, leaving  ˜µztand  ˜Azt, i.e. the mean and standard deviation of z⟨s⟩t, being deterministic functions. It guarantees that the gradient-based optimisation techniques are applicable by isolating the model parameters (  ˜µztand  ˜Azt) from the sampling procedure (involving  ϵ ⟨s⟩t).


Figure 3: Two components of the generative network.


By introducing hidden state hztand hxtas memory for his- torical information integration, the formulation is essentially equivalent to the recurrent model illustrated as Fig. 2.

We decompose the recurrent model into two components in a similar way as one would apply in factorising pΦ(X, Z)into the marginal distribution pΦ(Z)in Eq. (13) and the conditional distribution pΦ(X|Z)in Eq. (14).

The marginal pΦ(Z)is implemented by


which represents an autoregressive network for the latent ztas illustrated in Fig. 3a. The conditional pΦ(X|Z)is built as


which corresponds to a conditional generative network for the observable xtas in Fig. 3b.


Figure 4: Illustration of the training setup.


Figure 5: Architecture of the inference network.

On the other hand, the inference network is implemented in a similar recurrent fashion, as an autoregressive network with bidirectional dependencies:


The architecture of inference network is illustrated in Fig. 5.

The training procedure involves the inference network (in Fig. 5) and the conditional generative network (in Fig. 3b); the autoregressive network (in Fig. 3a) will not be utilised.

The historical observations  {x<t}is fed into the inference network, which outputs the inferred sequence  {z<t}of the causing latent variable. The latent sequence  {z<t}is then put into the conditional generative network (in Fig. 3b) to generate the predictions  { ˆx<t}. The likelihood of the predictions { ˆx<t}regarding the actual observations  {x<t}is calculated and the networks are optimised via gradient-based methods. Figure 4 illustrates the setting of networks during training.

The procedure of forecasting contains three steps: 1. feeding the historical data as the observations  {x<t}into the inference network (in Fig. 5) to infer the latent variables  {z<t}that might have caused that observations, 2. evolving the latent dynamics using the autoregressive network (in Fig. 3b) with the inferred sequence  {z<t}to produce the latent variable ztfor the next time step, and 3. invoking the conditional generative network (in Fig. 3b) with  {z1:t}and  {x<t}to predict the next possible xt. The procedure is shown in Fig. 6.

So far we have kept the covariance matrix  Σin our formulas to allow for multivariate forecasting. It entails a complexity of computation of order  O(M3)to maintain and update the full-size covariance matrix  Σof M dimensions; for cases in high dimensions, the computational cost for the full-size  Σbecomes unaffordable. Thus, we have to seek for alternatives that are more economic in terms of computation. A practical


Figure 6: Illustration of the forecasting setup.


approach is to leverage low-rank perturbations on diagonal matrices D such that  Σ−1 = D +VV⊤, where V  = [v1:K]is the rank-K perturbation with each  vkbeing independent M-dimensional column vector. The corresponding covariance matrix and its determinant can be readily calculated using Woodbury identity and matrix determinant lemma:


To solve the standard deviation A in the factorisation  Σ =AA⊤, we start with rank-1 perturbation. Given K = 1, matrix V = [v] is basically a column vector, hence I  + V⊤D−1V =1 + v⊤D−1vreturns a real number. A solution of A reads


where  γ = v⊤D−1vand  η = (1 + γ)−1. The complexity of computation is just of order O(M). Observe that VV⊤ =�Kk=1 vkv⊤k, the perturbation of rank K is merely the super- position of K rank-1 perturbations. Thus, we can calculate A in a recurrent fashion, where the complexity of computation for rank-K perturbation remains O(M) when K  ≪M holds. Algorithm 2 describes the detailed scheme of computation.

Designed for Accessibility and to further Open Science