Deep Stochastic Volatility Model

The GARCH (Bollerslev 1986) assumes a linear dependence of the future variance $\sigma_{t}^{2}$ on the previous $p$ squared returns and previous $q$ variances, computed as as follows:

Volatility clustering is captured by the parameters $\alpha_{i},i=1,\cdots,p$. However, the GARCH model assumes symmetric effects of negative and positive returns. (Glosten, Jagannathan, and Runkle 1993) propose the GJR-GARCH model which extends the GARCH model to account for the asymmetry effects of positive and negative returns as follows:

where the asymmetry effect is captured by the indicator function $I\left\{r_{t-i}<0\right\}$ and controlled by the leverage parameters $\gamma_{i}$. Similarly, the Threshold GARCH (TGARCH) proposed by (Zakoian 1994) also use leverage parameter $\gamma_{i}$ to reflect the asymmetry effects. The difference is that the TGARCH directly works with the volatility $\sigma_{t}$ instead of the variance $\sigma_{t}^{2}$:

Another popular extension of the GARCH model is the Exponential GARCH (EGARCH, (Nelson 1991)), which models the evolution of the log variance $\log(\sigma_{t}^{2})$:

where the leverage parameter $\gamma_{i}$ is used to account for the asymmetry effect.

For analysis convenience, stochastic volatility models are usually formulated with stochastic differential equations in continuous-time (Hull and White 1987; Heston 1993). The Heston model for a univariate process is

where $s(t)$ is the logrithm of stock price at $t$, $W^{\sigma}(t)$ and $W^{s}(t)$ are two correlated Wiener processes, with correlation represented as $E[\mathrm{d}W^{\sigma}(t)\mathrm{d}W^{s}(t)]=\rho\mathrm{d}t$. For practical use, stochastic volatility can be formulated in a discrete-time fashion for regularly spaced data. The discrete stochastic volatility without leverage (S. Kim, N. Sheppard, and Chibb 1997) is defined as

In order to accommodate the asymmetric effect of positive and negative returns, (Omori et al. 2007) propose the discrete stochastic volatility with leverage as follows:

where the leverage $\rho<0$ is trying to capture the increase in volatility that follows a drop in the return.

(Wu, Lobato, and Ghahramani 2014) proposed a Gaussian process volatility model (GP-Vol):

where $f$ can be modeled as a sample path of a Gaussian Process and $z_{t}$ is the stochastic noise. The model can be trained with Sequential Monte Carlo method.

The neural stochastic volatility model (NSVM) proposed by (Luo et al. 2018) can be formulated as follows:

where $f$ is a deep neural network model. As will be detailed in the following, one of the main differences of our DSVM model and the NSVM is that we use a different generative probabilistic model to account for the effect of past volatility and model the volatility as $\sigma_{t}=f_{\sigma}(\sigma_{1:t-1},r_{1:t-1},z_{1:t})$. Another difference is that, (Luo et al. 2018) model the volatility utilizing the information from multiple stock return time series which are correlated to each other and thus utilize more information. However, how to properly group multiple stocks is not a easy task in the first place. In this paper, we aim to propose a volatility model using less information, i.e. only based on the knowledge of the univariate stock return series.

Besides, there exist other machine learning volatility models using more exogenous variables, such as sentiment data (Xing, Cambria, and Zhang 2019) and using other advanced deep learning architectures, see (Zhang et al. 2018).

Denote $q_{\phi}\left(z_{1:T}|r_{1:T}\right)$ as any approximated posterior for the true posterior $p_{\theta}\left(z_{1:T}|r_{1:T}\right)$, it can be derived that

When $q_{\phi}\left(z_{1:T}|r_{1:T}\right)$ is equal to the true posterior $p_{\theta}\left(z_{1:T}|r_{1:T}\right)$, the lower bound is tight and we have $\mathcal{L}(\theta)=ELBO(\theta,\phi)$. However, the true posterior is intractable. In order to make the lower bound $ELBO$ tight, we design the inference network which represents $q_{\phi}\left(z_{1:T}|r_{1:T}\right)$ to best approximate the true posterior $p_{\theta}\left(z_{1:T}|r_{1:T}\right)$. In fact, the true posterior can be factorized as

The factorization can be obtained from the Figure 1 using d-seperation (Geiger, Verma, and Pearl 1990). The factorization shows that, the posterior of ${z_{t}}$ depends on the past information encoded in $z_{t-1}$ as well as the future information in ${r}_{t:T}$. We can design our inference network according to this factorization form of the true posterior. Specifically,

In addition, we have the following two factorization

The past information $r_{1:t-1},{z}_{1:t}$ is reflected in $\sigma_{t}$. Therefore, the conditional density of return at time $t$ on the past returns and the stochastic noises can be replaced by the conditional density of return on its simultaneous volatility.

The ELBO then can be derived as

where $q_{\phi}^{*}({z}_{t})$ is the marginal probability density of ${z}_{t}$:

We approximate the ELBO in Equation (3) using the Monte Carlo estimation. Specifically, we samples from $q_{\phi}^{*}\left({z}_{t}\right)$ using ancestral sampling according to the Equation (2). Given a sample $z_{t-1}^{(s)}$ from $q_{\phi}^{*}\left(z_{t-1}\right)$, a sample ${z}_{t}^{(s)}$ from $q_{\phi}\left(z_{t}|z_{t-1}^{(s)},{r}_{t:T}\right)$ follows the distribution of $q_{\phi}^{*}\left({z}_{t}\right)$. After we get the samples $z_{t}^{(s)}\text{ for }t=1\cdots,T$ from the $q_{\phi}\left(z_{1:T}|r_{1:T}\right)$, the ELBO in (3) can be approximated with the generated samples as:

Recall that, instead of factorizing $q_{\phi}$ as mean-field approximation across different time step, i.e. $z_{t},t=1,\cdots,T$ are assumed to independent and does not use the information of the dependence structure in the generative network, we choose the variational posterior to have the same factorization as the true posterior in Equation (2), where the posterior of $z_{t}$ depend on $z_{t-1}$ and $r_{t:T}$.

We additionally assume that the information in $r_{t:T}$ can be encoded with a backward RNN, denoted as $g_{A}$ with hidden state, denoted as $A_{t}$, as follows:

where $g_{A}$ is a RNN model. The graphical model of the inference network is shown in Figure 2. $z_{t-1}$ represents the information coming from the past, whereas $A_{t}$ encodes the information from present and the future, thus the information from all time steps are used to approximate the posterior at each $t$. Further, $q_{\phi_{z}}\left(z_{t}|z_{t-1},A_{t}\right)$ is parameterized to be a Gaussian distribution density:

$z_{t}$ is a Gaussian variable with parameter $m^{(q)}_{t},v^{(q)}_{t}$ given by MLP models $g_{1}$ and $g_{2}$ that depend on $z_{t-1},A_{t}$. The parameters in the inference network include parameters in the following functions: $\phi=\{g_{1},g_{2},g_{A}\}$.

The derivative of the $ELBO(\theta,\phi)$ with respect to $\theta$ can be calculated as:

The derivative of the $ELBO(\theta,\phi)$ with respect to $\phi$ is more tricky as $\phi$ appears in the expectation in $ELBO(\theta,\phi)$. Score function gradient estimator (Glynn 1987; Williams 1992) can be used to approximate the gradient, but the obtained gradients suffer from high variance. Reparameterization trick (Rezende, Mohamed, and Wierstra 2014; Kingma and Welling 2013) is often used to obtain a lower variance gradient estimator instead. Specifically, in order to obtain a sample $z_{t}$ from $N(m_{t}^{(q)},(v_{t}^{(q)})^{2})$, we first generate a sample $\eta_{t}\sim N(0,1)$, and then use $m_{t}^{(q)}+\eta_{t}v_{t}^{(q)}$ as a sample of $z_{t}$, denote as $z_{t}=k_{t}(\eta_{t})$. The gradients then can be backpropogated through the gaussian random variables and $\nabla_{\phi}ELBO(\theta,\phi)$ can be calculated as

The above inference and learning algorithm for training DSVM are summarized in Algorithm 1.

After training the model, we can use the estimated model to do one-step-ahead forecasting. Specifically, the probability of $r_{T+1}$ conditional on the information of $r_{1:T}$ can be derived as:

We approximate the $p(z_{1:T}|r_{1:T})$ with $q(z_{1:T}|r_{1:T})$ and use ancestral sampling to generate samples from the approximated $p(r_{T+1}|r_{1:t})$. Specifically, we generate $\{z_{1:T}^{(s)}\}_{s=1}^{S}$ from $q(z_{1:T}|r_{1:T})$ and then generate samples $\{z_{T+1}^{(s)}\}_{s=1}^{S}$ from $p(z_{T+1}|z_{T})$ based on $\{z_{T}^{(s)}\}_{s=1}^{S}$. The approximated predictive distribution for $p(r_{T+1}|r_{1:t})$ is then a mixture of $S$ Gaussian models. We can then generate $\{r_{T+1}^{(s)}\}_{s=1}^{S}$ from the predictive distribution (5) and use the sample standard deviation of $\{r_{T+1}^{(s)}\}_{s=1}^{S}$ as the predicted volatility for $r_{T+1}$.

We choose the daily adjusted return of the 39 biggest stocks in terms of market values among SP500 stocks from 2001-01-02 to 2018-12-31. For each stock, the training and validation data are from 2001-01-02 to 2015-10-18 (in total of 3721 days) and the testing data are from 2015-10-19 to 2018-12-31 (in total of 806 days). We extract sequences of 10 timestep data (corresponding to half of a trading month) from the original univariate time series and the training, validation, testing data size is 105924, 31434 and 31473 respectively (a proportion of around 6:2:2).

We train one model for all the 39 stocks. The dimension of the latent variable $z$ is chosen to be 1. $f_{1},f_{2},f_{3},g_{1},g_{2}$ are all parameterized by two-hidden-layer MLP models with different activation functions at last layer (Linear for mean function $f_{1}$, $g_{1}$, Softplus for standard deviation function $f_{2},f_{3},g_{2}$). $f_{h},g_{{A}}$ are both parameterized by one-layer GRU models with 10 hidden states. Monte Carlo samples size used for forecasting is 1000. The model is trained on a Tesla V100 GPU for 1 hours (300 epochs) before it converges. The final model is chosen by minimizing the validation loss over the 300 epochs. The DSVM is only trained once on the training and validation data, and then is used to do recursive one-day-ahead forecasting without retraining.

We compare the performance of the propsed model with seven popular alternative volatility models described in Section 2.

1. 1.

   Deterministic volatility models:
   We consider four models: GARCH, GJR-GARCH, TGARCH and EGARCH. The parameter estimation is done with a widely used R package “rugarch” ¹¹1https://cran.r-project.org/web/packages/rugarch.

2. 2.

   Stochastic volatility models:
   We choose two models: the discrete stochastic volatility model without leverage and discrete stochastic volatility model with leverage, and use the R package “stochvo” ²²2https://cran.r-project.org/web/packages/stochvol to estimate the parameters.

3. 3.

   Gaussian process volatility model (GP-Vol)
   The implementation of the model is obtained from the authors’ website³³3https://bitbucket.org/jmh233/code-gpvol-nips2014/src/master and the hyperparameters are chosen as suggested in (Wu, Lobato, and Ghahramani 2014).

For each type of the alternatives, we train one model for each stock as the alternatives are designed to be so. We will perform recursive one-day-ahead forecasting and all the alternative models are retrained at every forecasting time step with a rolling window of 1000.

We first evaluate each model in terms of negative loglikelihood (NLL) for the testing samples for each stock. The smaller the NLL value, the better the model can explain the distribution of the observed testing samples. The results of NLL evaluations are listed in Table 1. We first apply the Friedman rank sum test and the test indicates that there exist significant differences of the NLL values among different models (the p-value is less than $e^{-35}$). Furthermore, a pairwise pos thoc Nemenyi-Friedman test (Demšar 2006) is performed to compare different models, and it confirms that the DSVM performs best among all the methods. Our DSVM is the only stochastic volatility model that is able to beat deterministic volatility models and is the best for 36 out of the 38 stocks.

We further plot the predicted volatility of DSVM for training data and testing data and compare with the volatility produced by the EGARCH model (the best performing alternative model in terms of mean NLL). We use the two largest stocks (in terms of market value), the MSFT and the AAPL, for illustration. Other stocks perform similarly and are omitted in the manuscript due to page limit. The predicted volatility for training and testing data are displayed in Figure 3 (with a randomly selected zoom-in period) and Figure 4 respectively. In general, the changes of the DSVM is more adaptive to the market than the EGARCH model, which is reflected by several spikes in the predicted volatility of DSVM. In addition, the DSVM produces higher predicted volatility when the market is more volatile (see the circle areas) and gives lower volatility when the market is more stable (see the rectangle areas) for both training and testing data. In contrast, the predicted volatility of EGARCH is more smooth and the changes of the market status could not be quickly captured. Take the AAPL stock for example (see Figure 4), the volatility of the EGARCH decays more slowly around 2016 June when the market was calm, while it increases more slowly on 2018 September when the market became more volatile. In summary, the DSVM produces predicted volatility more responsive to the market status, which is more desirable as a risk measure for the financial market.