A Two-Step Framework for Arbitrage-Free Prediction of the Implied Volatility Surface

The dynamics of the IVS has been analyzed in various papers. Skiadopoulos et al. (2000), Fengler et al. (2003) and Cont and Da Fonseca (2002) apply principal component analysis (PCA) to IVS data. The first two papers perform PCA for smiles of different maturities while the third reference performs PCA directly on the surface. It is found in Cont and Da Fonseca (2002) that the first three eigenmodes already explain most of the variations and they are associated with the level, skewness and convexity of the IVS. The authors approximate the IVS on each day using the linear combination of the first three eigensurfaces and the dynamics of each coefficient is modelled by a first-order autoregressive model. Fengler et al. (2007) develop a semiparametric factor model for the IVS, where they assume the IVS is given by a sum of basis functions, each of which uses moneyness and time to maturity as inputs. The factor loadings are modelled by a vector autoregressive process. The basis functions and the loadings are estimated by weighted least squares with the weight given by a kernel function. Bloch (2019) proposes a general modeling framework. He uses several risk factors to represent the surface and model each parameter (corresponding to a risk factor) using a neural network. In particular, he provides three ways of obtaining the risk factors: using the first three eigenmodes from PCA in Cont and Da Fonseca (2002) or the parameters of a polynomial model for the IVS, or the parameters of a stochastic volatility type model (such as the SVI model). Explanatory variables such as the spot price, volume, VIX, etc can be input into the neural networks for these parameters. Cao et al. (2020) model the relationship between the expected daily change in the implied volatility of S&P500 index options by a multilayer feedforward neural network with daily index return, VIX, moneyness and time to maturity as inputs. They find that with the index return and VIX as features, their model can improve an analytic model significantly. One can use all of these models for predicting the IVS on a future date, but their predicative performance is not assessed in these papers. Furthermore, all these papers do not show how to obtain the dynamics of the IVS without arbitrage.

Several other papers directly address the dynamic prediction problem. Dellaportas and Mijatović (2014) predict the implied volatilities of a single maturity by forecasting the parameters in the SABR model from a time series model. As the SABR model is arbitrage free, the predicted implied volatility smile has no arbitrage. However, their approach does not apply to the whole surface because the SABR model cannot fit the surface well. Goncalves and Guidolin (2006) and Bernales and Guidolin (2014) model the IVS as a polynomial of time-adjusted moneyness and time to maturity. To predict the IVS on a future date, they use the forecasted coefficients of the polynomial from a time series model. Audrino and Colangelo (2010) develop a regression tree model through boosting to predict the IVS. Chen and Zhang (2019) apply the long short-term memory (LSTM) model with attention mechanism to predict the IVS. Bloch and Böök (2020) predicts the IVS by predicting the risk factors driving the dynamics of the IVS using temporal difference backpropagation (TDBP) models. Zeng and Klabjan (2019) use tick data on options to construct the implied volatilty surface at high frequency through a support vector regression model. A common issue with all these papers is that the predicted IVS from their models is not necessarily arbitrage free. Recently, Ning et al. (2021) introduce an interesting approach to simulate arbitrage-free IVS over time. They first calibrate an arbitrage-free stochastic model to the IVS data and then generate the model parameters of a future date from a variational autoencoder. The future IVS is then obtained under the stochastic model using the generated model parameters. However, they don’t use their approach to predict future IVS.

The contributions of this paper are twofold. First, we provide a general framework for dynamic prediction and simulation of IVS free of static arbitrage. This is a new feature provided by our approach compared with existing methods for dynamic prediction of IVS. Second, we show how to construct features to represent the IVS and develop some successful models for predicting the IVS in this framework.

Our framework consists of two steps.

• •

  Step 1: We select features to represent the IVS and predict them. The predicted features are mapped to a discrete set of implied volatilities. If the task is simulation, we simulate the features in this step.

• •

  Step 2: We construct the entire IVS from the discrete set of implied volatilities in Step 1 through a deep neural network (DNN) model by taking into account the constraints that rule out static arbitrage.

This framework is completely flexible as users can construct features and predict them in their own ways. Furthemore, any results obtained in Step 1 can be converted to a static-arbitrage free surface through the DNN model in Step 2. However, the predicted surface from our framework may admit dynamic arbitrage opportunities as we don’t enforce constraints that prevent dynamic arbitrage in our model. It would be difficult to do so in our framework as it is data-driven and assumes no stochastic model for the underlying asset.

The accuracy of predicting the IVS obviously hinges on the selected features and the model for prediciting them. In general, one can use features extracted from the IVS data together with exogenous variables to represent the surface. The focus of this paper is on how to extract features and we consider three approaches: using the eigenmodes from the PCA analysis of Cont and Da Fonseca (2002), applying the variational autoencoder (VAE) to extract latent factors for the IVS, and directly sampling the IVS on a discrete grid of moneyness and time to maturity. PCA can be viewed as a parametric approach that approximates the change of the log-IVS using a linear combination of eigensurfaces. The VAE approach is more general than the PCA as it offers a flexible nonlinear factor representation of the surface. The sampling approach is completely nonparametric.

To predict the extracted features, we utilize the long short-term memory (LSTM) model, which is a popular deep learning model for sequential data. Our choice is motivated by the success of deep learning in a range of prediction problems in finance. See e.g., Borovykh et al. (2017) and Sezer et al. (2020) for various financial time series, Sirignano (2019) and Sirignano and Cont (2019) for limit order books, Sirignano et al. (2018) for mortgage risk and Yan et al. (2018) for tail risk in asset returns.

By training our models using data of 9.5 years and putting them to test in a period of 2.5 years, we find that both the sampling and the VAE approach are quite successful in predicting the IVS. The error of the PCA approach is almost three times of the other two, indicating that prediction based on a linear combination of eigensurfaces is not accurate enough.

Another important finding is that the DNN model in Step 2 not only serves the purpose of constructing an arbitrage free surface but it is also crucial for prediction accuracy. Compared with a standard interpolation method for the IVS, using the DNN model can reduce the out-of-sample prediction error substantially for the sampling and VAE approach.

Our paper is related to Bloch (2019) which also provides a general and appealing framework, but they differ in a number of ways. First, Bloch (2019) doesn’t consider how to obtain arbitrage free surfaces. Second, Bloch (2019) constructs features for the IVS using PCA and parametric models (such as the polynomial or SVI model). We offer another two approaches for feature construction in this paper. Our empirical results suggest that the prediction model based on features from PCA is not good enough and we also have some potential issues with predicting features from parametric models (see Remark 1). Third, Bloch (2019) proposes to model each feature by a neural network. We model these features jointly by one model (LSTM in our implementation). Having separate neural network models for the features may fail to capture potential dependence among them unless they are independent. Finally, Bloch (2019) doesn’t provide empirical study to validate his models.

The rest of the paper is organized as follows. Section 2 provides background information on the implied volatility surface, including the definition of implied volatility, conditions that ensure no static arbitrage, an interpolation method and our data. Section 3 presents the two-step framework for prediction and simulation. Section 4 shows empirical results and compares different models. Section 5 concludes with remarks for future research.

Consider a European call option on a dividend paying asset $S_{t}$ with maturity date $T$ and strike price $K$. Set $\tau=T-t$, which is the time to maturity. Denote the risk-free rate by $r$ and the dividend yield by $d$. Let $F_{t}(\tau)$ be the forward price at $t$ for time to maturity $\tau$. It is given by $F_{t}(\tau)=S_{t}\mathrm{e}^{(r-q)\tau}$. We will write $F_{t}$ below to simplify the notation.

Under the Black-Scholes model, the call option price at time $t$ is given by

where

and $N(\cdot)$ is the cumulative distribution function of the standard normal distribution. The variable $m$ defined in (3) is known as the log forward moneyness.

It is well-known that the Black-Scholes model is unrealistic. To use it in practice, one looks for the level of volatility to match an observed option price, i.e., we solve $\sigma$ from the equation

where $C_{mkt}$ is the observed market price for the call option. The solution is called the implied volatility.

The implied volatility surface at a time point is the collection of implied volatilities for options with different $K$ and $\tau$. We prefer to consider the IVS as a function of $m$ and $\tau$, because $m$ is a relative coordinate. Hereafter, the IVS at time $t$ is denoted by $\sigma_{t}(m,\tau)$. Practitioners also like to quote implied volatilities using the Black-Scholes Delta (detnoted by $\delta$) and $\tau$, where $\delta=e^{-q\tau}N(d_{1})$.

Conditions that ensure the implied volatility surface is free of static arbitrage have been well studied in the literature (see e.g., Roper (2010), Gulisashvili (2012)). We summarize them in the proposition below.

Condition 3 implies that $\sigma(m,\tau)$ is free of calendar spread arbitrage and condition 4 guarantees that there is no butterfly arbitrage. In Section 3.4, we will implement these conditions to get an arbitrage free surface.

On a given day, implied volatilities can only be calculated for a discrete set of $(m,\tau)$ pairs, which correspond to options that are listed on that day. Suppose we are given $\{\sigma(m_{i},\tau_{i}):i=1,\cdots,n\}$ on a day. There are various approaches to interplate these given points to obtain the entire IVS. Here, we consider a simple and popular parametric model proposed in Dumas et al. (1998) and hereafter it will simply be called DFW. This model assumes

where a floor of $0.01$ is set to prevent the implied volatility of the model from becoming too small or even negative. The coefficients $a_{0},\cdots,a_{5}$ can be estimated by regression.

Another popular non-parametric approach to estimate the entire implied volatility surface uses the Nadaraya–Watson (NW) estimator (Härdle (1990)), which is given by

where

is a Gaussian kernel. The estimator involves two hyper-parameters $h_{1}$ and $h_{2}$, which are bandwidths and they determine the degree of smoothing. If the parameters are too small, a bumpy surface is generated. If they are too big, important details in the observed data can be lost after smoothing. When implementating this approach, on each day we apply five-fold cross-validation to the implied volatility data of this day to select the best pair of $(h_{1},h_{2})$ from a grid of values for them.

Table 1 presents the average root mean squared error (RMSE) of interpolation of these two methods, where the average is taken over days in our training period. The DFW model is more accurate and it will be our choice for further study. The larger error of the NW estimator is very likely caused by applying the same bandwidth $(h_{1},h_{2})$ to all points, whereas our implied volatility data is non-uniformly distributed in the $(m,\tau)$ space.

The dataset used for this paper contains implied volatility surfaces for the S&P500 index options from January 1, 2009 to December 31, 2020. We obtained the data from OptionMetrics through the Wharton Research Data Services. On each day, we have implied volatilities for a set of $(\delta,\tau)$ pairs with $\delta$ going from $0.1$ to $0.9$ with an increment of $0.05$ and $\tau=10,30,60,91,122,152,182,273,365,547,730$ calendar days. Since we consider the IVS as a function of $m$ and $\tau$, we convert $\delta$ to $m$ using

This results in implied volatilies on different days for different grids of moneyness but the same grid of $\tau$. We denote by $\mathcal{I}_{d,t}$ the set of $(m,\tau)$ pairs for observed implied volatilities at time $t$. In total, we have data on 3021 days and on each day 374 points are observed from the implied volatility surface (i.e., the size of $\mathcal{I}_{d,t}$ is 374). To demonstrate salient features of the implied volatility surface for index options, we calculate the average of $\sigma_{t}(\delta,\tau)$ over $t$ for all observed $(\delta,\tau)$ pairs, and plot the average values as a surface in Figure 1(a) in terms of $\delta$ and $\tau$. In Figure 1(b), we show the implied volatility curves for different maturities as functions of log forward moneyness. A volatility skew is clearly observed for each $\tau$ and it remains quite steep even for large maturities.

As the methods we will apply later cannot be used if the $(m,\tau)$-grid changes from day to day, we have to fix it. We use the following grid for $(m,\tau)$:

where

Since the set of $\tau$ is fixed over time in the data, we simply adopt the set as the grid for $\tau$ but change the time unit to year ($\tau$ was quoted in days initially). For the grid of moneyness, we first get the minimum value and maximum value of $m$ in our dataset and set $[\log(0.6),\log(2)]$ as the range, which is slightly wider than the range from the minimum to the maximum. Then we create a non-uniform grid on this range so that the grid is denser near ATM. As $\mathcal{I}_{0}$ is different from the observed grid for $(m,\tau)$ on a day, we use the DFW model given in (7) to obtain the implied volatilities on $\mathcal{I}_{0}$. Eventually, at every $t$, we have a set of 154 implied volatilities

which can be deemed as a sample of the implied volatility surface.

We consider several methods to extract features from the implied volatility data.

While having a high-dimensional feature vector can better approximate the surface, predicting it may be more difficult. Thus, natually we can consider some dimension reduction techniques to extract features, which leads us to the following two methods.

The loss function for training the VAE has two parts. The first part is the mean squared loss between the synthetic data and the original data given by

where $M$ is the batch size. The second part is the Kullback-Leibler (KL) divergence between the parameterized normal distribution and $N(0,I)$ given by

where $\mu_{k}$ and $\sigma_{k}$ are the mean and standard deviation of the $k$-th latent variable. The loss function is defined as

Adding the KL divergence term encourages the model to encode a distribution that is as close to normal as possible and the hyperparameter $\beta$ measures the extent of regularization.

In our problem, $Y_{t}=\bar{\Sigma}_{t}$ and we set $Z_{t}=\mu(Y_{t})$. With the predicted $\hat{Z}_{T+1}$, $F_{T+1}=h(\hat{Z}_{T+1})=N_{D}(\hat{Z}_{T+1})$, i.e., the $h$ function is given by the decoder FNN.

To predict $Z_{T+1}$ from $\{Z_{0},\cdots,Z_{T}\}$, one can consider all kinds of models. In our experiment, we will use the long short-term memory (LSTM) model (Hochreiter and Schmidhuber (1997)), which is a popular deep learning model for sequential data prediction and its success has been demonstrated in many problems. We use the model in the following way for our problem:

• •

  For any $T$, consider

  $$Z_{T}^{1}=\frac{1}{22}\sum_{t=T-21}^{T}Z_{t},\ Z_{T}^{2}=\frac{1}{5}\sum_{t=T-4}^{T}Z_{t},\ Z_{T}^{3}=Z_{T}.$$

  (24)

  The first two represent monthly and weekly moving averages at $T$, respectively. We predict $Z_{T+1}$ using $Z^{1}_{T},Z^{2}_{T},Z^{3}_{T}$, which can be viewed as long-term, medium-term and short-term features. A similar approach is taken by Corsi (2009) and Chen and Zhang (2019).

• •

  Let $h_{j}$ be a hidden state that represents a summary of information from $\{Z_{T}^{1},\cdots,Z_{T}^{j}\}$. Set $h_{0}=0$. For $j=1,2,3$, calculate

  	$\displaystyle r_{j}$	$\displaystyle=\sigma_{g}\left(W_{r}Z_{T}^{j}+U_{r}h_{j-1}+b_{r}\right),$		(25)
  	$\displaystyle i_{j}$	$\displaystyle=\sigma_{g}\left(W_{i}Z_{T}^{j}+U_{i}h_{j-1}+b_{i}\right),$		(26)
  	$\displaystyle o_{j}$	$\displaystyle=\sigma_{g}\left(W_{o}Z_{T}^{j}+U_{o}h_{j-1}+b_{o}\right),$		(27)
  	$\displaystyle g_{j}$	$\displaystyle=\sigma_{h}\left(W_{g}Z_{T}^{j}+U_{g}h_{j-1}+b_{g}\right),$		(28)
  	$\displaystyle c_{j}$	$\displaystyle=r_{j}\odot c_{j-1}+i_{j}\odot g_{j},$		(29)
  	$\displaystyle h_{j}$	$\displaystyle=o_{j}\odot\sigma_{h}\left(c_{j}\right),$		(30)
  	$\displaystyle y_{j}$	$\displaystyle=\sigma_{h}\left(W_{y}h_{j}+b_{y}\right).$		(31)

  All $W$, $U$, $b$ are parameters and $\sigma_{g}$, $\sigma_{h}$ are the sigmoid function and the tanh function, respectively, for activation. At $j$, $i_{j}$, $r_{j}$ and $o_{j}$ represent input gate, forget gate, and output gate.

• •

  Finally, we predict $Z_{T+1}$ as

  $$\hat{Z}_{T+1}=\sigma_{out}(W_{out}y_{3}+b_{out}).$$

  (32)

  The range of $Z_{T+1}$ varies in our framwork depending on the feature extraction method. For SAM, we use relu for $\sigma_{out}$ as $Z_{T+1}$ is positive. For VAE and PCA, since $Z_{T+1}$ can take any real value, we don’t use any nonlinear activitation function and simply set $\sigma_{out}$ as the identity function.

In the following, we will write the feature prediction model as

where $\theta_{\mathcal{P}}$ is the vector of parameters involved.

With $F_{T+1}=h(\hat{Z}_{T+1})$, we construct the entire implied volatility surface from $F_{T+1}$ using a DNN illustrated in Figure 4. The neural network is a feedforward one with inputs $F_{T+1},m,\tau$. The output is an implied volatility for the input $(m,\tau)$ pair. We use the Softplus function $\ln(1+\exp(x))$ as the activation function of the output layer, because it makes the output nonnegative and twice differentiable, so that the first two no-arbitrage conditions in Proposition 1 are fulfilled.

Suppose the time horizon in our data is given by $\mathcal{T}=\left\{1,2,...,T\right\}$ and let $q_{t}$ be the number of observed implied volatilities on the surface at $t$ (in our data $q_{t}=374$ for all $t$ but in general it could change over time). The loss function of the featuer prediction part is given by

We minimize $\mathcal{L}_{\mathcal{P}}(\theta_{\mathcal{P}})$ to train the LSTM model. For the construction of the implied volatility surface, one can set the loss function as

However, minizing $\mathcal{L}_{\mathcal{S}}(\theta_{\mathcal{S}})$ to train the surface construction model cannot guarantee the output surface is arbitrage free. By design the output of the DNN model satisfies the first two conditions in Proposition 1, but does not necessarily fulfil the other three. Inspired by Zheng et al. (2019) and Ackerer et al. (2020), we incorporate Conditions 3,4,5 for no arbitrage into our training by formulating them as penalties in the loss function.

First, we create the following synthetic grids to faciliate the calculation of the penalty functions:

where $m_{\min}=\log(0.6)$, $m_{\max}=\log(2)$, $[a,b]_{40}$ means a uniform grid over the interval $[a,b]$ with it divided into $40$ equal parts,

and $\tau_{\max}=730/365$. The grid $\mathcal{I}_{\mathrm{C}34}$ is used for the penalty calculation associated with Condition 3 and 4 while $\mathcal{I}_{\mathrm{C}5}$ is used for Condition 5. These grids are different from the $(m,\tau)$-grid in (11) used for sampling. In particular, $\mathcal{I}_{\mathrm{C}34}$ has 1600 points which is a lot more than the 154 points on the sampling grid and it also covers a much wider range for both $m$ and $\tau$. We use such a dense grid on a wide region to reduce the chance of missing points on the surface at which there is significant violation of the no-arbitrage conditions. As Condition 5 considers the large moneyness behavior, we analyze moneyness levels which are extremely negative or positive.

We denote by $\mathcal{L}_{\mathrm{C}j}(\theta_{\mathcal{S}})$ the penalty function for the $j$-th condition ($j=3,4,5$). For Conditions 3 and 4, they are given by

where $\ell_{\mathrm{cal}}\left(m_{i},\tau_{i},F_{t};\theta_{\mathcal{S}}\right)$ and $\ell_{\mathrm{but}}\left(m_{i},\tau_{i},F_{t};\theta_{\mathcal{S}}\right)$ are defined as in (5) and (6) with $\sigma$ replaced $f$. For Condition 5, it is equivalent to that the second-order derivative of $\sigma^{2}(m,\tau)$ goes to zero as $|m|\to\infty$, where $\partial^{2}_{mm}\sigma^{2}(m,\tau)=\sigma(m,\tau)\partial^{2}_{mm}\sigma(m,\tau)+(\partial_{m}\sigma(m,\tau))^{2}$. Hence, the penalty is

Finally, we obtain our loss function for training the DNN model as

for some $\lambda>0$. One could use a separate penalization parameter for each penalty, but for simplicy we assume they are the same. We minimize $\mathcal{L}_{\mathrm{C}}(\theta_{\mathcal{S}})$ to train the DNN model. In our implementation, we choose $\lambda=1$, which is used in Ackerer et al. (2020) in their penalized loss function. We also tried other values for $\lambda$ and found that using $\lambda=1$ results in the smallest error for the IVS on the training data and the penalties converge to zero quickly.

Our framework can also be used to simulate the IVS over time. We can write the feature transition equation as

or

where $\varepsilon_{T+1}$ is the error vector at $T+1$. In (43) we assume additive error and in (44) we assume multiplicative error. The multiplicative formulation is more convenient to use than the additive one when the positivity of $Z_{T+1}$ is required.

We assume the error process $\varepsilon_{1},\varepsilon_{2},\cdots$ is an i.i.d. white noise with mean zero and covariance matrix $\Sigma_{\varepsilon}$. After obtaining the estimate of $\theta_{\mathcal{P}}$ by minimizing the loss function $\mathcal{L}_{\mathcal{P}}(\theta_{\mathcal{P}})$, one can calculate the error vector on each day and hence obtain a sample for the errors. We can assume the error vector follows a multivariate parametric distribution $F_{\eta}$ with parameter vector $\eta$ (e.g., Gaussian) and estimate $\eta$ from the error sample.

The simulation of $Z_{T+1}$ given the available information at $T$ consists of the following steps.

• •

  Step 1: Calculate $p(Z^{1}_{T},Z^{2}_{T},Z^{3}_{T};\theta_{\mathcal{P}})$.

• •

  Step 2: Simulate $\varepsilon_{T+1}$ from $F_{\eta}$ or by bootstrapping from the error sample.

• •

  Step 3: Calculate $Z_{T+1}$ by (43) or (44).

• •

  Step 4: Calculate $\sigma_{T+1}(m,\tau)=f(m,\tau,h(Z_{T+1}))$.

The DNN model $f$ ensures that the output IVS is free of static arbitrage.

The details of the three feature extraction methods can be found in Section 3.1. For each day in the dateset, we extract features using these three methods and some details are as follows:

• •

  For SAM, we use $\bar{\Sigma}_{t}$ as the feature vector (see (14)), which is a set of implied volatilities on a $(m,\tau)$-grid with 154 points, to represent the entire surface at $t$.

• •

  For PCA, we follow Cont and Da Fonseca (2002) to use the first three eigenmodes (i.e., $K=3$ in (20)), which already explain over $96\%$ of the variations in our data.

• •

  For VAE, the FNNs for the encoder and the decoder both have three hidden layers with 128 nodes per layer. We try five values for the latent dimension $d$: $2,5,10,15,20$. Their performance on the test data is shown in Table 4 and the difference is small, indicating the performance of the VAE model is quite robust to the choice of $d$. The VAE model with $d=10$ achieves the smallest out-of-sample prediction error.

We use the LSTM model to predict the extracted features as discussed in Section 3.2. For the DNN model for surface construction, we use three hidden layers with fifty neurons on each layer. In the training of both models, we do the following:

• •

  We initialize the parameters using Xavier initialization (Glorot and Bengio (2010)), which can prevent initial weights in a deep network from being either too large or too small. This method sets the weight of the $j$-th layer to follow a uniform distribution given by

  $$W^{j}\sim U\left[-\frac{1}{\sqrt{n_{j}}},\frac{1}{\sqrt{n_{j}}}\right],$$

  (45)

  where $n_{j}$ is the number of neurons on the $j$-th layer.

• •

  We use the Adam optimizer with minibatches (Kingma and Ba (2014)) to minimize the loss function. Calculating the gradient of the loss function using all the samples can be computationally expensive, so in each iteration we only use a minibatch (i.e., subset) of samples for the gradient evaluation. The Adam optimizer is a popular gradient-descent algorithm, which utilizes the exponentially weighted average of the gradients to accelerate convergence to the minimum.

• •

  We apply batch normalization to the inputs of the neural network (Ioffe and Szegedy (2015)). For all the samples in a minibatch, we first estimate the mean and standard deviation of each input in this minibatch and then normalize it by subtracting its estimated mean and dividing by its estimated standard deviation.

Values of the hyperparameters associated with training and the hidden size of the models (i.e., number of neurons on a hidden layer) are displayed in Table 2. We train LSTM and DNN for 200 and 20 epochs, respectively. An epoch consists of all the iterations required to work through all the samples in the training set, so it is given by the size of the training data divided by the size of a minibatch.

Figures 5 and 6 show the results of loss on the training data and test data as the number of epochs increases. In Figure 6, we only plot the DNN results for the model with features extracted from the sampling approach and results for the other two feature extraction approaches are similar. From Figure 5, we can see that there is no overfitting for the LSTM model as the test loss is close to the training loss. Similarly, there is no overfitting for the DNN model as shown by Figure 6(a). The value of the penalty for three no-arbitrage conditions also become zero in the test data eventually, so there is no violation of these conditions on the synthetic grids. It should be noted that although there are some spikes for the calender spread penalty in the training data, the largest value is still very small, so the violation is insignificant.

Let $\hat{\theta}_{\mathcal{P}}$ and $\hat{\theta}_{\mathcal{S}}$ be the estimated parameters from the training data for the LSTM model and the DNN model, respectively. Suppose the time index of the last day in the training period is $T_{\text{train}}$ and of the last day in the whole dataset is $T_{\text{total}}$. Set $T_{\text{test}}=T_{\text{total}}-T_{\text{train}}$, which is the number of days in the test period. We do out-of-sample test as follows: for every $t>T_{\text{train}}$,

• •

  obtain $\hat{Z}_{t}=p(Z^{1}_{t-1},Z^{2}_{t-1},Z^{3}_{t-1};\hat{\theta}_{\mathcal{P}})$ and $F_{t}=h(\hat{Z}_{t})$;

• •

  calculate $\hat{\sigma}_{t}(m,\tau)=f(m,\tau,F_{t})$ for $(m,\tau)\in\mathcal{I}_{d,t}$.

Here, $\mathcal{I}_{d,t}$ is the set of $(m,\tau)$-pairs in the observed implied volatility data at $t$, which contains 374 points. It is important to note that it is different from $\mathcal{I}_{0}$, the set of $(m,\tau)$-pairs used for sampling the surface, which has only 154 points. The error for a pair of $(m,\tau)$ is given by

where $\sigma_{t}(m,\tau)$ is the observed implied volatility at $t$ for this pair (the ground truth). The error not only reflects the prediction error of the LSTM model for the features, but also the interpolation error of the DNN model.

To evaluate the overall out-of-sample prediction performance, we consider two commonly used error measures: root mean squared error (RMSE) and the mean absolute percentage error (MAPE). They are defined as

In our data, $\mathcal{I}_{d,t}$ contains different $(m,\tau)$ pairs for a different $t$, but $|\mathcal{I}_{d,t}|=374$ for all $t$. ¹¹1One should be cautious in comparing the errors reported in different papers. Some papers only evaluate the error on a limited set of $(m,\tau)$ pairs. For example, Chen and Zhang (2019) only consider the errors at 45 points in the $(m,\tau)$ space.

We examine various models. In the first step, there are three feature extraction approaches: SAM, PCA and VAE. In the second step, we consider two methods: the DNN model and the DFW model in (7) applied to $F_{t}$ to predict $\sigma_{t}(m,\tau)$ at time $t$. The DNN model yields an arbitrage free surface whereas the DFW interpolation model cannot. This leads to six models for comparison: SAM-DNN, SAM-DFW, PCA-DNN, PCA-DFW, VAE-DNN, VAE-DFW. We also consider a classical benchmark given by the DFW model. At time $T$, we simply forecast the IVS at $T+1$ from the DFW model with its coefficients given by their estimates at $T$.

The performance of these models on the test dataset is shown in Table 3. For any pair of models, we also perform the Diebold-Mariano (DM) test (Diebold and Mariano (2002)) to assess the statistical significance of the difference in the forecast performance as measured by RMSE and the p-value is shown in Table 5. Consider model 1 and model 2. In the DM test, the null hypothesis is that the forecast error is equal while the alternative hypothesis is that the forecast error of model 1 is less than model 2. Table 5 should be read in the following way. For any entry of the table, the model on its row is model 1 and the model on its column is model 2. We consider 1% as the significance level.

Several observations can be made from Table 3 and 5. (1) The best performers in out-of-sample prediction are SAM-DNN and VAE-DNN. The DM test shows their difference is statistically insignificant, so they can be considered as equally good. Both of them outperfom the other models with overwhelmingly strong statistical evidence. In particular, these two models constructed in the proposed two-step framework beat the classical DFW model by a large margin. (2) The out-of-sample error of SAM-DNN and VAE-DNN is only about one third of the error of PCA-DNN. This result highlights the importance of feature selection in step 1 for predicting the IVS. While the PCA approach is good for understanding the main factors that drive the IVS movements, the approximation based on a linear combination of eigensurfaces is not sufficiently accurate for predicting the IVS. In contrast, the VAE approach is more flexible as it combines the latent factors in a nonlinear way. Thus, its improvement over PCA can be expected and this is confirmed by the results. The sampling approach can be deemed as a nonparametric way to represent the surface, which can lead to a high-dimensional feature vector (in our data its dimension is 154). Thanks to the power of LSTM, we are able to predict it quite accurately. Using a powerful model like LSTM for feature prediction is key to the success of the sampling approach. (3) The DNN model for IVS construction in step 2 also makes significant contributions to improving the prediction accuracy. For each feature extraction method, using the DNN model outperforms using the DFW model for surface construction in step 2 with statistical significance. The improvement in prediction accuracy is already considerable for SAM and even more substantial for VAE. (4) It’s worth noting that the SAM-DFW model also performs quite well in prediction. It’s simpler than the SAM-DNN model as it uses the simple DFW model instead of the complex DNN model for surface construction in step 2.

We further plot the RMSE and MAPE of each day in the test period for four models in Figure 7. Both SAM-DNN and VAE-DNN are better than PCA-DNN throughout the preiod and they also outperform DFW on most days. However, the error of all four models spikes up in March, 2020, during which the US stock market suffered a meltdown due to the pandemic. The relatively big error signals a shift in the market regime in that period. The features we use in all the models are extracted from the implied volatility data, which do not provide a direct representation of the market regime. To improve their performance, one can further augment the feature vector with exogenous variables like index return and VIX which are proxies of the market regime.

For the predicted IVS on the test days, we check for violation of the constraints on calendar spread arbitrage and butterfly arbitrage. Consider

The no-arbitrage constraints require $\ell_{\mathrm{cal}}(m,\tau)$ and $\ell_{\mathrm{but}}(m,\tau)$ to be nonnegative for any $(m,\tau)$. In the above quantities, we check $\ell_{\mathrm{cal}}(m,\tau)$ and $\ell_{\mathrm{but}}(m,\tau)$ on the $(m,\tau)$ grid for the observed implied volatilities for each test day. A negative value at a pair of $(m,\tau)$ indicates violation at this point and we calculate the average of the negative values over all test days. The results are reported for four models in Table 6. While no violation is detected for the three models that use DNN in step 2, prediction under the DFW model yields implied volatility surfaces with arbitrage opportunities and calendar spread arbitrage in particular.