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 xfollows from the Gaussian distribution with mean 0 and variance ; the (conditional) variance is fully de- termined by a linear function (Eq. (1)) of previous observations and variances . 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 Iis the indicator function: Iif x, and 0 otherwise, which allows for asymmetric reactions of volatility in terms of the sign of previous .

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 :

where (Eq. (5)) accommodates the asymmetric relation between observations and volatility changes. If setting q in 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 sat time t, and represent two correlated Wiener processes and the correlation between and is expressed as t.

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 depends on not only the historical log-variances but a latent stochastic process . The latent process is, according to Eq. (9), white noise process with i.i.d. Gaussian variables.

Notably, the volatility is no longer conditionally deter- ministic (i.e. deterministic given the complete history ) 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 and the latent stochastic process as . As seen in previous sections, the dynamics of the volatility process can be abstracted in the form of

The latter equality holds when we recurrently substitute with f for all t. For models within the GARCH family, we discard zin (Eq. (10)); on the other hand, for the primitive SV model, xis ignored instead. We can loosen the constraint that xis zero-mean to a time-varying mean for more flexibility.

Recall that the latent stochastic process (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 and variance 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 xand the latent variable zare normally distributed. Note that the autoregressive process degenerates to i.i.d. white noise process when and . It should be emphasised that the purpose of reinforcing an autoregressive structure (11) of the latent variable zis 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 zwith an autoregressive structure provides a possibility of decoupling the internal influential factors from the external ones, which we believe is the essence of introducing z.

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

Generating Observable Sequence Recall that the observable variable x(Eq. (12)) and the latent variable z(Eq. (11)) are described by autoregressive models (as xalso involves an exogenous input z). Let pand pdenote the probability distributions of xand zat time t. The factorisation on the joint distributions of sequences and applies as follow:

where X and Z 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 , the current mean and variance of zis obtained and hence the distribu- tion of zis specified; after sampling zfrom the specified distribution, we incorporate and calculate the current mean and variance of xand determine its distribution of x. 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 for the latent variable and for 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 hand hdenote 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 zis N(0, 1) in a way that zis 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 zis obtained, e.g. by inference (see details in the next subsection), the conditional distribution p(Eq. (14)) will be involved in generating the observable xinstead of the joint distribution p(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 z, of which the true values are inaccessible even we have observed x, the marginal distribution pbecomes the key that bridges the model and the data. However, the calcula- tion of pitself or its complement, the posterior distribution p, is often intractable as complex integrals are involved. We are unable to learn the parameters by differentiating the marginal log-likelihood or to infer the latent variables through the true posterior. Therefore, we consider instead a restricted family of tractable distributions q, referred to as the approximate posterior family, as approximations to the true posterior psuch 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 , for any T, zis 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 and are functions of the given observation sequence , representing the approximated mean and variance of the latent variable zdenotes 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 and represent the hidden states of the forward and backward directions of the bidirectional RNN. The autoregressive RNN with hidden state takes the joint state of the bidirectional RNN and the previous value of zas input. The inference mean and variance is computed by an MLP from the hidden state of 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 up to time stepT, 1-step-ahead prediction of either or xis 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 pinvolved in Eq. (28) is intractable, the exact evaluation of conditional predictive distribution pis difficult.

A straightforward solution is that we substitute the true posterior pwith the approximation q(see Eq. (22)) and leverage qto inference S sample paths of the latent variables according to the historical observations . 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(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 sinto log-returns xand 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 . Specifically, the actual dataset for training and evaluation comprises a collection of 2000 series of d-dimensional normalised log-return vectors of length 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 and the latent variables . 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 ’s output layer is of size 4 + 4 for whereas the size of ’s output layer is 6 + 6 for . 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 zinferred 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 , 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 , we target at maximising the marginal log-likelihood pw.r.t. , where the actual posterior pis involved. Because of the intractability of p, the exact inference is not applicable; we have to seek for approximate solutions instead.

We factorise the marginal log-likelihood as

Note that we have introduced a tractable, -parameterised distribution qfrom a flexible family of distributions to approximate the actual posterior p. The evidence lower bound (ELBO) 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 z. 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 denotes the collection of S sample paths.

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

where is the standard Gaussian variable.

Figure 2: Generalised stochastic volatility model.

The reparameterisation extracts the randomness out of the latent variable zvia , leaving and , i.e. the mean and standard deviation of z, being deterministic functions. It guarantees that the gradient-based optimisation techniques are applicable by isolating the model parameters ( and ) from the sampling procedure (involving ).

Figure 3: Two components of the generative network.

By introducing hidden state hand has 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 pinto the marginal distribution pin Eq. (13) and the conditional distribution pin Eq. (14).

The marginal pis implemented by

which represents an autoregressive network for the latent zas illustrated in Fig. 3a. The conditional pis built as

which corresponds to a conditional generative network for the observable xas 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 is fed into the inference network, which outputs the inferred sequence of the causing latent variable. The latent sequence is then put into the conditional generative network (in Fig. 3b) to generate the predictions . The likelihood of the predictions regarding the actual observations 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 into the inference network (in Fig. 5) to infer the latent variables that might have caused that observations, 2. evolving the latent dynamics using the autoregressive network (in Fig. 3b) with the inferred sequence to produce the latent variable zfor the next time step, and 3. invoking the conditional generative network (in Fig. 3b) with and to predict the next possible x. 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 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 VV, where V is the rank-K perturbation with each being 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 returns a real number. A solution of A reads

where and . The complexity of computation is just of order O(M). Observe that VV, 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.