Title: A NEURAL STOCHASTIC VOLATILITY MODEL

Abstract: In this paper, we show that the recent integration of statistical models with recurrent neural networks provides a new way of formulating volatility models that have been popular in time series analysis and prediction. 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 observable ones. Our focus in this paper is on the formulation of temporal dynamics of volatility over time under a stochastic recurrent neural network framework. Our derivations show that some popular volatility models are a special case of our proposed neural stochastic volatility model. Experiments demonstrate that the proposed model generates a smoother volatility estimation, and outperforms standard econometric models GARCH, EGARCH, GJR-GARCH and some other GARCH variants as well as MCMC-based model stochvol and a recent Gaussian processes based volatility model GPVOL on several metrics about the fitness of the volatility modelling and the accuracy of the prediction.

1. INTRODUCTION:

The volatility of the price movements reflects the ubiquitous uncertainty within financial markets. It is critical that the level of risk, 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 & Granger, 2003) .

Volatility is measured typically by using the standard deviation of price change in a fixed time interval, such as a day, a month or a year. The higher the volatility, the riskier the asset. One of the primary challenges in designing volatility models is to identify the existence of latent (stochastic) variables or 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 heteroskedasticity (ARCH) model (Engle, 1982) and its generalisation GARCH (Bollerslev, 1986) , which makes use of autoregression to capture the properties of time-variant volatility within many time series. Heston (1993) assumed that the volatility follows a Cox-Ingersoll-Ross (CIR) process (Cox et al., 1985) and derived a closed-form solution for options pricing. While theoretically sound, those approaches require strong assumptions which might involve complex probability distributions and non-linear dynamics that drive the process, and in practice, one may have to impose less prior knowledge and rectify a solution under the worst-case volatility case (Avellaneda & Paras, 1996) .

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). We are inspired by the recent development on variational approaches of stochastic (deep) neural networks (Kingma & Welling, 2013; Rezende et al., 2014) to a recurrent case (Chung et al., 2015; Fabius & van Amersfoort, 2014; Bayer & Osendorfer, 2014) , and our formulation shows that existing volatility models such as the GARCH (Bollerslev, 1986 ) and the Heston model (Heston, 1993) are the special cases of our neural stochastic volatility formulation. With the hidden latent variables in the neural networks we naturally uncover the underlying stochastic process formulated from the models.

Experiments with synthetic data and real-world financial data are performed, showing that the proposed model outperforms the widely-used GARCH model on several metrics of the fitness and the accuracy of time series modelling and prediction: it verifies our model's high flexibility and rich expressive power.

2. RELATED WORK:

A notable volatility method is autoregressive conditional heteroskedasticity (ARCH) model (Engle, 1982) : it can accurately capture the properties of time-variant 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 (1986) generalised ARCH model to the generalised autoregressive conditional heteroskedasticity (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 & 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 covariances to present and depend on the historical information, which are particularly useful in multivariate financial models. Heston (1993) derived a closed-form solution for option pricing with stochastic volatility where the volatility process is a CIR process driven by a latent Wiener process such that the current volatility is no longer a deterministic function even if the historical information is provided. Notably, empirical evidences have confirmed that volatility models provide accurate forecasts (Andersen & Bollerslev, 1998) and models such as ARCH and its descendants/variants have become indispensable tools in asset pricing and risk evaluation.

On the other hand, deep learning (LeCun et al., 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 et al., 2012; He et al., 2015; van den Oord et al., 2016) , speech recognition (Hinton et al., 2012; Graves et al., 2013; Chorowski et al., 2015) , machine translation (Sutskever et al., 2014; Cho et al., 2014; Bahdanau et al., 2014; Luong et al., 2015) to name a few. Time-dependent neural networks models include RNNs with advanced neuron structure such as long short-term memory (LSTM) (Hochreiter & Schmidhuber, 1997) , gated recurrent unit (GRU) (Cho et al., 2014) , and bidirectional RNN (BRNN) (Schuster & Paliwal, 1997) . Recent results show that RNNs excel for sequence modelling and generation in various applications (Graves, 2013; Gregor 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 posterori 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 & Welling, 2013; Rezende et al., 2014) . Some early work has started to explore the use of variational inference to make RNNs stochastic (Chung et al., 2015; Bayer & Osendorfer, 2014; Fabius & van Amersfoort, 2014) . Bayer & Osendorfer (2014) and Fabius & van Amersfoort (2014) considered the hidden variables are independent between times, whereas (Fraccaro et al., 2016) utilised a backward propagating inference network according to its Markovian properties. Our work in this paper extends the work (Chung et al., 2015) with a focus on volatility modelling for time series. We assume that the hidden stochastic variables follow a Gaussian autoregression process, which is then used to model both the variance and the mean. We show that the neural network formulation is a general one, which covers two major financial stochastic volatility models as the special cases by defining the specific hidden variables and non-linear transforms.

3. PRELIMINARY: VOLATILITY MODELS:

Stochastic processes are often defined by stochastic differential equations (SDEs), e.g. a (univariate) generalised Wiener process is d x t = µ d t + σ d w t , where µ and σ denote the time-invariant rates of drift and standard deviation (square root of variance) while d w t ∼ N (0, d t) is the increment of standard Wiener process at time t. In a small time interval between t and t + ∆t, the change in the variable is ∆x t = µ∆t + σ∆w t . Let ∆t = 1, we obtain the discrete-time version of basic volatility model:

EQUATION

where t ∼ N (0, 1) is a sample drawn from standard normal distribution. In the multivariate case, Σ represents the covariance matrix in place of σ 2 . As presumed that the variables are multidimensional, we will use Σ to represent variance in general case except explicitly noted.

3.1. DETERMINISTIC VOLATILITY:

The time-invariant variance Σ can be extended to be a function Σ t = Σ(x <t ) relying on history of the (observable) underlying stochastic process {x <t }. The current variance Σ t is therefore determined given the history {x <t } up to time t. An example of such extensions is the univariate GARCH(1,1) model (Bollerslev, 1986) :

EQUATION

where x t-1 is the observation from N (µ t-1 , σ 2 t-1 ) at time t -1. Note that the determinism is in a conditional sense, which means that it only holds under the condition that the complete history {x <t } is presented, such as the case of 1-step-ahead forecast. otherwise the current volatility would still be stochastic as it is built on stochastic process {x t }. However, for multi-step-ahead forecast, we usually exploit the relation

E t-1 [(x t -µ t ) 2 ] = σ 2

t to substitute the corresponding terms and calculate the forecasts with longer horizon in a recursive fashion, for example,

σ 2 t+1 = α 0 + α 1 E t-1 [(x t - µ t ) 2 ] + β 1 σ 2 t = α 0 + (α 1 + β 1 )σ 2 t .

For n-step-ahead forecast, there will be n iterations and the procedure is hence also deterministic.

3.2. STOCHASTIC VOLATILITY:

Another extension is applicable for Σ t from being conditionally deterministic (i.e. deterministic given the complete history {x <t }) to fully stochastic: Σ t = Σ(z ≤t ) is driven by another latent stochastic process {z t } instead of the observable process {x t }. Heston (1993) model instantiates a continuous-time stochastic volatility model for univariate processes:

EQUATION

EQUATION

where the correlation between d w

1 t and d w 2 t applies: E[d w 1 t •d w 2 t ] = ρ d t.

We apply Euler's scheme of quantisation (Stoer & Bulirsch, 2013) to obtain the discrete analogue to the continuoustime Heston model (Eqs. ( 3) and ( 4)):

x t = (x t-1 + µ -0.5σ 2 t ) + σ t σ t = (1 + a)σ t-1 + bz t where t z t = N (0, 1 ρ ρ 1

).

(5)

3.3. VOLATILITY MODEL IN GENERAL:

As discussed above, the observable variable x t follows Gaussian distribution of which the mean and variance depend on the history of observable process {x t } and latent {z t }. We presume in addition that the latent process {z t } is an autoregressive model such that z t is (conditionally) Gaussian distributed. Therefore, we formulate the volatility model in general as:

EQUATION

where µ z (x <t , z ≤t ) and Σ z (x <t , z ≤t ) denote the autoregressive time-varying mean and variance of the latent variable z t while µ x (x <t , z ≤t ) and Σ x (x <t , z ≤t ) represent the mean and variance of observable variable x t , which depend on not only history of the observable process {x <t } but that of the latent process {z ≤t }.

These two formulas (Eqs. ( 6) and ( 7)) abstract the generalised formulation of volatility models. Together, they represents a broad family of volatility models with latent variables, where the Heston model for stochastic volatility is merely a special case of the family. Furthermore, it will degenerate to deterministic volatility models such as the well-studied GARCH model if we disable the latent process.

4. NEURAL STOCHASTIC VOLATILITY MODELS:

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

4.1. GENERATING OBSERVABLE SEQUENCE:

Recall that the latent variable z t (Eq. ( 6)) and the observable x t (Eq. ( 7)) are described by autoregressive models (x t has the exogenous input {z ≤t }.) For the distributions of {z t } and {x t }, the following factorisation applies:

EQUATION

EQUATION

where X = {x t } and Z = {z t } are the sequences of observable and latent variables, respectively, while Φ represents the parameter set of the model. The full generative model is defined as the joint distribution:

EQUATION

)

It is observed that the means and variances are conditionally deterministic: given the historical information {z <t }, the current mean µ z t = µ z Φ (z <t ) and variance Σ z t = Σ z Φ (z <t ) of z t is obtained and hence the distribution N (z t ; µ z t , Σ z t ) of z t is specified; after sampling z t from the specified distribution, we incorporate {x <t } and calculate the current mean µ x t = µ x Φ (x <t , z ≤t ) and variance

Σ x t = Σ x Φ (x <t , z ≤t ) of x t and determine its distribution N (x t ; µ x t , Σ x t ) of x t .

It is natural and convenient to present such a procedure in a recurrent fashion because of its autoregressive nature. As is known that RNNs can essentially approximate arbitrary function of recurrent form (Hammer, 2000) , the means and variances, which may be driven by complex non-linear dynamics, can be efficiently computed using RNNs.

It is always a good practice to reparameterise the random variables before we go into RNN architecture. As the covariance matrix Σ is symmetric and positive definite, it can be factorised as Σ = U ΛU , where Λ is a full-rank diagonal matrix with positive diagonal elements. Let A = U Λ 1 2 , we have Σ = AA . Hence we can reparameterise the latent variable z t (Eq. ( 6)) and observable x t (Eq. ( 7)):

EQUATION

where

A z t (A z t ) = Σ z t , A x t (A x t ) = Σ x t and x t ∼ N (0, I x ), z t ∼ N (0, I z )

are auxiliary variables. Note that the randomness within the variables of interest (e.g. z t ) is extracted by the auxiliary variables (e.g. t ) which follow the standard distributions. Hence, the reparameterisation guarantees that gradient-based methods can be applied in learning phase (Kingma & Welling, 2013) .

In this paper, the joint generative model is comprised of two sets of RNN and multilayer perceptron (MLP): RNN z g /MLP z g for the latent variable, while RNN x g /MLP z g for the observables. We stack these two RNN/MLP together according to the causal dependency between those variables. The joint generative model is implemented as the generative network:

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

where h z t and h x t denote the hidden states of the corresponding RNNs. The MLPs map the hidden states of RNNs into the means and deviations of variables of interest. The parameter set Φ is comprised of the weights of RNNs and MLPs.

One should notice that when the latent variable z is obtained, e.g. by inference (details in the next subsection), the conditional distribution p Φ (X|Z) (Eq. ( 9)) will involve in generating the observable x t instead of the joint distribution p Φ (X, Z) (Eq. ( 10)). 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 joint generative model or the conditional one as it can be inferred in context.

4.2. INFERENCING THE LATENT PROCESS:

As the generative model involves latent variable z t , of which the true valus are unaccessible even we have observed x t . Hence, the marginal likelihood p Φ (X) becomes the key that bridges the model and the data. The calculation of marginal likelihood involves the posterior distribution p Φ (Z|X), which is often intractable as complex integrals are involved. We are unable to learn the paramters or to infer the latent variables. 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 flexible to provide good approximations (Bishop, 2006; Kingma & Welling, 2013; Rezende et al., 2014) .

We define the inference model in accordance with the approximate posterior family we have presumed, in a similar fashion as (Chung et al., 2015) , where the factorised distribution is formulated as follows:

EQUATION

where μz Ψ (z <t , x <t ) and Σz Ψ (z <t , x <t ) are functions of the historical information {z <t }, {x <t }, representing the approximated mean and variance of the latent variable z t , respectively. Note that Ψ represents the parameter set of inference model.

The inference model essentially describes an autoregressive model on z t with exogenous input x t . Hence, in a similar fashion as the generative model, we implement the inference model as the inference network using RNN/MLP:

EQUATION

EQUATION

EQUATION

where Ãz t ( Ãz t ) = Σz t = Σz Ψ (z <t , x <t ) while hz t represents the hidden state of RNN and ˜ z t ∼ N (0, I z ) is an auxiliary variable to extract randomness. The inference mean μz t and deviation Ãz t is computed by an MLP from the hidden state hz t . We use the subscript i instead of g to distinguish the architecture used in inference model in contrast to generative model.

4.3. FORECASTING OBSERVATIONS IN FUTURE:

In the realm of time series analysis, we usually pay more attention on forecasting over generating (Box et al., 2015) . It means that we are essentially more interested in the generation procedure conditioning on the historical information rather than generation purely based on a priori belief since the observations in the past of x <t influences our belief of the latent variable z t . Therefore, we apply the approximate posterior distribution of the latent variable z t (Eq. ( 19)) as discussed in previous subsection, in place of the prior distribution (Eq. ( 8)) to build our predictive model.

Inference Generate x t+1 h t+1 x z t+1 ht+1 z x t+1 x t h t x z t ht z x t x t-1 h t-1 x z t-1 μ t-1 z Σ t-1 z μ t z Σ t z μ t+1 z Σ t+1 z μ t-1 x Σ t-1 x μ t x Σ t x μ t+1 x Σ t+1 x ht-1 z x t-1

Given the historical observations x <t , the predictive model infers the current value of latent variable z t using inference network and then generates the prediction of the current observation x t using generative network. The procedure of forecasting is shown in Fig. 1 . NSVM is learned using Stochastic Gradient Variational Bayes following (Kingma & Welling, 2013; Rezende et al., 2014) . For readability, we provide the detailed derivation in Appendix A.

4.4. LINKS TO GARCH(1,1) AND HESTON MODEL:

Although we refer to GARCH and Heston as volatility models, the purposes of them are quite different: GARCH is a predictive model used for volatility forecasting whereas Heston is more of a generative model of the underlying dynamics which facilitate closed-form solutions to SDEs in option pricing. The proposed NSVM has close relations to GARCH(1,1) and Heston model: both of them can be regarded as a special case of the neural network formulation. Recall Eq. (2), GARCH(1,1) is formulated as

σ 2 t = α 0 + α 1 (x t-1 -µ t-1 ) 2 + β 1 σ 2 t-1

, where µ t-1 is the trend estimate of {x t } at time step t calculated by some mean models. A common practice is to assume that µ t follows the ARMA family (Box et al., 2015) , or even simpler, as a constant that µ t ≡ µ. We adopt the constant trend for simplicity as our focus is on volatility estimation.

We define the hidden state as h x t = [µ, σ t ] , and disable the latent variable z t ≡ 0 as the volatility modelled by GARCH(1,1) is conditionally deterministic. Hence, we instantiate the generative network (Eqs. ( 16), ( 17) and ( 18)) as follows:

EQUATION

EQUATION

x t = µ + σ t t where t ∼ N (0, 1).

The set of generative parameters is Φ = {µ, α 0 , α 1 , β 1 }.

Next, we show the link between NSVM and (discrete-time) Heston model (Eq. ( 5)). Let h x t = [x t-1 , µ, σ t ] be the hidden state and z t be i.i.d. standard Gaussian instead of autoregressive vari-able, we represent the Heston model in the framework of NSVM as:

EQUATION

EQUATION

EQUATION

EQUATION

)

The set of generative parameters is Φ = {µ, a, b}.

One should notice that, in practice, the formulation may change in accordance with the specific architecture of neural networks involved in building the model, and hence a closed-form representation may be absent.

5. EXPERIMENTS:

In this section, we present our experiments1 both on the synthetic and real-world datasets to validate the effectiveness of NSVM.

5.1. BASELINES AND EVALUATION METRICS:

To evaluate the performance of volatility modelling, we adopt the standard econometric model GARCH(1,1) Bollerslev (1986) as well as its variants EGARCH(1,1) Nelson (1991), GJR-GARCH(1,1,1) Glosten et al. (1993) , ARCH(5), TARCH(1,1,1), APARCH(1,1,1), AGARCH(1,1,1), NAGARCH(1,1,1), IGARCH(1,1), IAVGARCH(1,1), FIGARCH(1,d,1) as baselines, which incorporate with the corresponding mean model AR(20). We would also compare our NSVM against a MCMC-based model "stochvol" and the recent Gaussian-processes-based model "GPVOL" Wu et al. (2014) , which is a non-parametric model jointly learning the dynamics and hidden states via online inference algorithm. In addition, we setup a naive forecasting model as an alternative baseline referred to as NAIVE, which maintains a sliding window of size 20 on the most recent historical observations and forecasts the current values of mean and volatility by the average mean and variance of the window.

For synthetic data experiments, we take four metrics into consideration for performance evaluation: 1) the negative log-likelihood (NLL) of observing the test sequence with respect to the generative model parameters; 2) the mean-squared error (MSE) between the predicted mean and the ground truth (µ-MSE), 3) MSE of the predicted variance against the true variance (σ-MSE); 4) smoothness of fit, which is the standard deviation of the differences of succesive variance estimates. As for the real-world scenarios, the trend and volatility are implicit such that no ground truth is accessible to compare with, we consider only NLL and smoothness as the metrics for evaluation on real-world data experiment.

5.2. MODEL IMPLEMENTATION:

The implementation of NSVM in experiments is in accordance with the architecture illustrated in Fig. 1 : it consists of two neural networks, namely inference network and generative network. Each network comprises a set of RNN/MLP as we have discussed above: the RNN is instantiated by stacked LSTM layers whereas the MLP is essentially a 1-layer fully-connected feedforward network which splits into two equal-sized sublayers with different activation functions -one sublayer applies exponential function to impose the non-negativity and prevents overshooting of variance estimates while the other uses linear function to calculate mean estimates. During experiment, the model is structured by cascading the inference network and generative network as depicted in Fig. 1 . The input layer is of size 20, which is the same as the embedding dimension D E ; the layer on the interface of inference network and generative network -we call it latent variable layer -represents the latent variable z, where its dimension is 2. The output layer has the same structure as the input one, therefore the latent variable layer acts as a bottleneck of the entire architecture which helps to extract the key factor. The stacked layers between input layer, latent variable layer and output layer are the hidden layers of either inference network or generative network, it consists of 1 or 2 LSTM layers with size 10, which contains recurrent connection for temporal dependencies modelling.

State-of-the-art learning techniques have been applied: we introduce Dropout (Zaremba et al., 2014) into each LSTM recurrent layer and impose L2-norm on the weights of each fully-connected feedforward layer as regularistion; NADAM optimiser (Dozat, 2015) is exploited for fast convergence, which is a variant of ADAM optimiser (Kingma & Ba, 2014) incorporated with Nesterov momentum; stepwise exponential learning rate decay is adopted to anneal the variations of convergence as time goes.

For econometric models, we utilise several widely-used packages for time series analysis: statsmodels (http://statsmodels.sourceforge.net/), arch (https://pypi.python. org/pypi/arch/3.2), Oxford-MFE-toolbox (https://www.kevinsheppard. com/MFE_Toolbox), stochvol (https://cran.r-project.org/web/packages/ stochvol) and fGarch (https://cran.r-project.org/web/packages/fGarch). The implementation of GPVOL is retrived from http://jmhl.org and we adopt the same hyperparameter setting as in Wu et al. (2014) .

5.3. SYNTHETIC DATA EXPERIMENT:

We build up the synthetic dataset by generating 256 heteroskedastic univariate time series, each with 2000 data points i.e. 2000 time steps. At each time step, the observation is drawn from a Gaussian distribution with pre-determined mean and variance, where the tendency of mean and variance is synthesised as linear combinations of sine functions. Specifically, for the trend and variance, we synthesis each using 3 sine functions with randomly chosen amplitudes and frequencies; then the value of the synthesised signal at each timestep is drawn from a Gaussian distribution with the corresponding value of trend and variance at that timestep. A sampled sequence is shown in Fig. 2a . We expect that this limited dataset could well simulate the real-world scenarios: one usually has very limited chances to observe and collect a large amount of data from time-invariant distributions. In addition, it seems that every observable or latent quantity within time series varies from time to time and seldom repeats the old patterns. Hence, we presume that the tendency shows long-term patterns and the period of tendency is longer than observation. In the experiment, we take the former 1500 time steps as the training set whereas the latter 500 as the test set.

For the synthetic data experiment, we simplify the recurrent layers in both inference net and generative net as single LSTM layer of size 10. The actual input { x t } fed to NSVM is D Edimensional time-delay embedding (Kennel et al., 1992) of raw univariate observation {x t } such that x t = [x t+1-D E , . . . , x t ]. 2-dimensional latent variable z t is adopted to capture the latent process, and enforces an orthogonal representation of the process by using diagonal covariance matrix. At each time step, 30 samples of latent variable z t are generated via reparameterisation (Eq. ( 22)).

5.4. REAL-WORLD DATA EXPERIMENT:

We select 162 out of more than 1500 stocks from Chinese stock market and collect the time series of their daily closing prices from 3 institutions in China. We favour those with earlier listing date of trading (from 2006 or earlier) and fewer suspension days (at most 50 suspension days in total during the period of observation) so as to reduce the noise introduced by insufficient observation or missing values, which has significant influences on the performance but is essentially irrelevant to the purpose of volatility forecasting. More specifically, the dataset obtained contains 162 time series, each with 2552 data points (7 years). A sampled sequence is shown in Fig. 2b . We divide the whole dataset into two subsets: the training subset consists of the first 2000 data points while the test subset contains the rest 552 data points.

Similar model configuration is applied to the real-world data experiment: time-delay embedding of dimension D E on the raw univariate time series; 2-dimensional latent variable with diagonal covariance matrix; 30 sampling for the latent variable at each time step. Instead of single LSTM layers, here we adopt stacked LSTM layers composed of 2 × 10 LSTM cells.

5.5. RESULT AND DISCUSSION:

The overall performance of NSVM and baselines is listed in details in Table 1 and case studies on synthetic data and real-world financial data are illustrated in Fig. 2 . The results show that NSVM has higher accuracies for modelling heteroskedastic time series on various metrics: NLL shows the fitness of the model under likelihood measure; the smoothness indicates that NSVM obtains more robust representation of the latent volatility; µ-MSE and σ-MSE in synthetic data experiment imply the ability of recognising the underlying patterns of both trend and volatility, which in fact verifies our claim of NSVM's high flexibility and rich expressive power for volatility (as well as trend) modelling and forecasting compared with the baselines. 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 procedure. The newly proposed NSVM outperforms standard econometric models GARCH(1,1), EGARCH(1,1), GJR-GARCH(1,1,1) and some other variants as well as the MCMC-based model "stochvol" and the recent GP-based model "GPVOL". Apart from the higher accuracy NSVM obtained, it provides us with the ability to simply generalise univariate time series analysis to multivariate cases by extending network dimensions and manipulating the covariance matrices. Furthermore, it allows us to implement and deploy a similar framework on other applications, for example signal processing and denoising. The shortcoming of NSVM comparing to GPVOL is that the training procedure is offline: for short-term prediction, the experiments have shown the accuracy, but for long-term forecasting, the parameters need retraining, which will be rather time consuming. The online algorithm for inference will be one of the work in the future.

Specifically, our NSVM outperforms GARCH(1,1) on 142 out of 162 stocks on the metric of NLL.

In particular, NSVM obtains -2.111, -2.044, -2.609 and -1.939 on the stocks corresponding to 

6. CONCLUSION:

In this paper, a novel volatility model NSVM has been proposed for stochastic volatility estimation and forecast. We integrated statistical models and RNNs, 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 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. Our NSVM has presented higher accuracy compared to GARCH(1,1), EGARCH(1,1) and GJR-GARCH(1,1,1) as well as GPVOL for volatility modelling and forecasting on synthetic data and real-world financial data. Future work on NSVM would be to incorporate well-established models such as ARMA/ARIMA and to investigate the modelling of seasonal time series and correlated sequences.

As we have known, for models that evolve explicitly in terms of the squares of the residuals (e 2 t = (x t -µ t ) 2 ), e.g. GARCH, the multi-step-ahead forecasts have closed-form solutions, which means that those forecasts can be efficiently computed in a recursive fashion due to the linear formulation of the model and the exploitation of relation E t-1 [e 2 t ] = σ 2 t . On the other hand, for models that are not linear or do not explicitly evolve in terms of e 2 , e.g. EGARCH (linear but not evolve in terms of e 2 ), our NSVM (nonlinear and not evolve in terms of e 2 ), the closed-form solutions are absent and thus the analytical forecast is not available. We will instead use simulation-based forecast, which uses random number generator to simulate draws from the predicted distribution and build up a pre-specified number of paths of the variances at 1 step ahead. The draws are then averaged to produce the forecast of the next step. For n-step-ahead forecast, it requires n iterations of 1-step-ahead forecast to get there.

NSVM is designed as an end-to-end model for volatility estimation and forecast. It takes the price of stocks as input and outputs the distribution of the price at next step. It learns the dynamics using RNN, leading to an implicit, highly nonlinear formulation, where only simulation-based forecast is available. In order to obtain reasonably accurate forecasts, the number of draws should be relatively large, which will be very expensive for computation. Moreover, the number of draws will increase exponentially as the forecast horizon grows, so it will be infeasible to forecast several time steps ahead. We have planned to investigate the characteristics of NSVM's long-horizontal forecasts and try to design a model specific sampling method for efficient evaluation in the future.

The corresponding covariance matrix and its determinant is obtained using Woodbury identity and matrix determinant lemma:

Σ = D -1 -D -1 V (I + V D -1 V ) -1 V D -1 (33)

ln det Σ = -ln det (D + V V ) = -ln det D -ln det (I + V D -1 V ) (34)

To calculate the deviation A for the factorisation of covariance matrix Σ = AA , we first consider the rank-1 perturbation where K = 1. It follows that V = v is a column vector, and I + V D -1 V = 1 + v D -1 v is a real number. A particular solution of A is obtain:

EQUATION

)

where γ = v D -1 v, η = (1 + γ) -1 . The computational complexity involved here is merely O(M ).

Observe that V V = K k=1 v k v k , the perturbation of rank K is essentially the superposition of K perturbations of rank 1. Therefore, we can calculate the deviation A iteratively, an algorithm is provided to demonstrate the procedure of calculation. The computational complexity for rank-K perturbation remains to be O(M ) given K M .

Algorithm 1 gives the detailed calculation scheme.

Algorithm 1 Calculation of rank-K perturbation of precision matrices Input: The original diagonal matrix D; The rank-K perturbation V = {v 1 , . . . , v K } Output: A such that the factorisation AA = Σ = (D + V V ) -1 holds 1: A (0) = D -1 2 2: i = 0 3: while i < K do 4:

γ (i) = v (i) A (i) A (i) v (i)

5:

η (i) = (1 + γ (i) ) -1 6:

A (i+1) = A (i) -[γ -1 (i) (1 - √ η (i) )]A (i) A (i) v (i) v (i) A (i) 7: A = A (K)

:

Repeatable experiment code: https://github.com/xxj96/nsvm

Figures and Tables:

Figure 1: Forecasting the future using Neural Stochastic Volatility Model.

Synthetic time series prediction. (up) The data and the predicted µ x and bounds µ x ± σ x . (down) The groundtruth data variance and the corresponding prediction from GARCH(1,1) and NSVM. Real-world stock price prediction. (up) The data and the predicted µ x and bounds µ x ± σ x . (down) The variance prediction from GARCH(1,1) and NSVM. The prediction of NSVM is more smooth and stable than that of GARCH(1,1), also yielding smaller NLL.

Figure 2: A case study of time series prediction.

Fig2(b), Fig 4(a), (b) and (c) respectively, each of which is better than the that of GARCH (0.3433, 0.589, 0.109 and 0.207 lower on NLL).

Results of the experiments.

Table content: <table><tr><td/><td/><td colspan="2">SYNTHETIC DATA</td><td/><td>STOCK DATA</td></tr><tr><td/><td>NLL</td><td>µ-MSE</td><td colspan="2">σ-MSE smoothness</td><td>NLL smoothness</td></tr><tr><td>NSVM</td><td>3.932e-2</td><td colspan="2">2.393e-3 6.178e-4</td><td colspan="2">4.322e-3 -2.184</td><td>3.505e-3</td></tr><tr><td>GARCH(1,1)</td><td colspan="3">6.905e-2 7.594e-3* 8.408e-4</td><td colspan="2">4.616e-3 -1.961</td><td>6.659e-3</td></tr><tr><td colspan="4">GJRGARCH(1,1,1) 6.491e-2 7.594e-3* 7.172e-4</td><td colspan="2">4.426e-3 -2.016</td><td>4.967e-3</td></tr><tr><td>EGARCH(1,1)</td><td colspan="3">5.913e-2 7.594e-3* 8.332e-4</td><td colspan="2">4.546e-3 -2.001</td><td>5.451e-3</td></tr><tr><td>ARCH(5)</td><td colspan="3">7.577e-2 7.594e-3* 1.610e-3</td><td colspan="2">5.880e-3 -1.955</td><td>7.917e-3</td></tr><tr><td>TARCH(1,1,1)</td><td colspan="3">6.365e-2 7.594e-3* 7.284e-4</td><td colspan="2">4.727e-3 -2.012</td><td>3.399e-3</td></tr><tr><td>APARCH(1,1,1)</td><td colspan="3">6.187e-2 7.594e-3* 9.115e-4</td><td colspan="2">4.531e-3 -2.014</td><td>4.214e-3</td></tr><tr><td>AGARCH(1,1)</td><td colspan="3">6.311e-2 7.594e-3* 9.543e-4</td><td colspan="2">4.999e-3 -2.008</td><td>5.847e-3</td></tr><tr><td>NAGARCH(1,1,1)</td><td colspan="3">1.134e-1 7.594e-3* 9.516e-4</td><td colspan="2">4.904e-3 -2.020</td><td>5.224e-3</td></tr><tr><td>IGARCH(1,1)</td><td colspan="3">6.751e-2 7.594e-3* 9.322e-4</td><td colspan="2">4.019e-3 -1.999</td><td>4.284e-3</td></tr><tr><td>IAVGARCH(1,1)</td><td colspan="3">6.901e-2 7.594e-3* 7.174e-4</td><td colspan="2">4.282e-3 -1.984</td><td>4.062e-3</td></tr><tr><td>FIGARCH(1,d,1)</td><td colspan="3">6.666e-2 7.594e-3* 1.055e-3</td><td colspan="2">5.045e-3 -2.002</td><td>5.604e-3</td></tr><tr><td>MCMC-stochvol</td><td colspan="3">0.368 7.594e-3* 3.956e-2</td><td colspan="2">6.421e-4 -0.909</td><td>1.511e-3</td></tr><tr><td>GPVOL</td><td colspan="3">1.273 7.594e-3* 6.457e-1</td><td colspan="2">4.142e-2 -2.052</td><td>5.739e-3</td></tr><tr><td>NAIVE</td><td>2.037e-1</td><td colspan="2">8.423e-3 3.515e-3</td><td colspan="2">2.708e-2 -0.918</td><td>7.459e-3</td></tr><tr><td colspan="3">*the same results obtained from AR(20) mean models</td><td/><td/></tr></table>

Peer review (iclr2017):
The authors propose a recurrent neural network approach for constructing a stochastic volatility model for financial time series. They introduce an inference network based on a recurrent neural network that computes the approximation to the posterior distribution for the latent variables given the past data. This variational approximation is used to maximize the marginal likelihood in order to learn the parameters of the model. The proposed method is validated in experiments with synthetic and real-world time series, showing to outperform parametric GARCH models and a Gaussian process volatility model. Quality: The method proposed seems technically correct, with the exception that in equation (19) the inference model is doing filtering and not smoothing, in the sense that the posterior for z_t* only depends on those other z_t and x_t values with t<t*, but in the true posterior p(Z|X) each z_t depends on all the X. This means the proposed learning method is inefficient. The reference below shows how to perform smoothing too and the results in that paper show indeed that smoothing produces better results for learning the model. Sequential Neural Models with Stochastic Layers Fraccaro, Marco and So nderby, So ren Kaae and Paquet, Ulrich and Winther, Ole In NIPS 2016. It is not clear if the method proposed in the above reference would perform better just because of using smoothing when learning the model parameters. Clarity: The paper is clearly written and easy to read. For the results on real-world data, in Table 1, how is the NLL computed? Is the average NLL across the 162 time series? Originality: The method proposed is not very original. Previous work has already used the variational approach with the reparametrization trick to learn recurrent neural networks with stochastic units (see the reference above). It seems that the main contribution of the authors is to apply this type of techniques to the problem of modeling financial time series. Significance: The results shown are significant, the method proposed by the authors outperforms previous approaches. However, there is a huge amount of techniques available for modeling financial time-series. The number of GARCH variants is probably close to hundreds, each one claiming to be better than the others. This makes difficult to quantify how important the results are.