# DEEP WEIGHTED MONTE CARLO: A HYBRID OPTION PRICING FRAMEWORK USING NEURAL
NETWORKS

Sándor Kunsági-Máté∗, Gábor Fáth, István Csabai, Gábor Molnár-Sáska

ELTE Eötvös Loránd University, Budapest, Hungary
*Corresponding author
E-mail<EMAIL_ADDRESS>

ABSTRACT. Recent studies have demonstrated the efficiency of Variational
Autoencoders (VAE) to compress high-dimensional implied volatility surfaces
into a low dimensional representation. Although this method can be effectively
used for pricing vanilla options, it does not provide any explicit information
about the dynamics of the underlying asset. In our work we present an
effective way to overcome this problem. We use a Weighted Monte Carlo approach
to first generate paths from a simple a priori Brownian dynamics, and then
calculate path weights to price options correctly. We develop and successfully
train a neural network that is able to assign these weights directly from the
latent space. Combining the encoder network of the VAE and this new "weight
assigner" module, we are able to build a dynamic pricing framework which
cleanses the volatility surface from irrelevant noise fluctuations, and then
can price not just vanillas, but also exotic options on this idealized vol
surface. This pricing method can provide relative value signals for option
traders.

Key findings of the paper:

  * •

We confirmed the success of Variational Autoencoders in compressing volatility
surfaces into low dimensional latent spaces in the swaption market.

  * •

We successfully extracted the physically meaningful information from the
latent representation of the market by translating the latent coordinates into
weights of Monte Carlo paths.

  * •

We developed a hybrid approach by combining Deep Learning and Weighted Monte
Carlo that can price any type of option consistently with the latent model
represented by the fully trained VAE.

Keywords: Volatility surface - Deep Learning - Option pricing - Variational
Autoencoders - Monte Carlo

##  1 Introduction

In the early ‘70s, the Black-Scholes model (Black and Scholes, (1973)) put the
pricing of options on a rigorous conceptual and mathematical basis. Assuming
that the underlying asset volatility is constant was a reasonable
simplification to make, and indeed, it was a good approximation for equity
markets before the 1987 Black Monday crash. Since then the empirical
observation is that options at different strike prices (K) imply different
volatilities. This behaviour of the market is called the volatility skew or
volatility smile. During the ‘90s several local volatility models appeared
which were generalisations of Black-Scholes and tried to capture the observed
volatility smile (see e.g. Dupire et al., (1994), Derman et al., (1996), Dumas
et al., (1998)). While these models could well reproduce the static pattern of
the smiles they were unable to predict their dynamics. An effective solution
to this problem was provided by stochastic volatility models (Heston, (1993),
Kim et al., (1998), Aguilar and West, (2000)). One frequently used approach
proposed by Hagan and coworkers is the SABR stochastic volatility model (Hagan
et al., (2002)). SABR uses two coupled stochastic differential equations to
describe the time development of the underlying asset price and its
volatility. However, the dynamics of the volatility smile is still not
described properly by the SABR model either. Traders have to recalibrate
relatively often the SABR parameters (especially $\rho$ and $\nu$) to fit the
model to observed data.

As in the case of many fields of science, models used in the financial
industry became more and more complicated, having larger and larger degrees of
freedom to describe a specific problem with increasing accuracy. The usual
model development process consists of two main steps: 1) Find the most
relevant parameters of the problem and create a mathematical model that
describes the relation between them; 2) Compare the model behaviour with the
observations, and go back to step 1 if there is a non-negligible discrepancy.
This clearly shows that during this iterative process it is really hard to
find the number of minimally required parameters and the correct relations
among them.

Neural networks made a big impact on several areas of sciences and everyday
life in recent years, and brought a completely different approach to solve
problems. Deep Learning provides a data-driven approach where the model is
adapted directly to the observed data. We do not have to care about all
possible events anymore because the neural network can itself learn the
complex relation between the real parameters, hence providing a better
representation of the problem. Neural networks are flexible non-linear
functions that can be trained to find an adequate low dimensional
representation behind the data sometimes called the active subspace (Tripathy
and Bilionis, (2018)).

There are two main types of applications of neural networks in option pricing:
model replication and model calibration. Several stochastic models do not have
an explicit formula for option prices, and therefore we need to run
computationally expensive Monte Carlo simulations. Neural networks are able to
find the connection between the model parameters and the outcome, and hence
they can replicate the model providing an approximation function that can be
evaluated immediately (Liu et al., (2019), Pagnottoni, (2019), Gaspar et al.,
(2020)). Another common problem is the calibration of the stochastic models to
the observed market. In practice, finding the best fitting pricing model means
to run again the time-consuming Monte Carlo simulation for several different
sets of model parameters. Neural networks are also able to predict the model
parameters directly from the observations (Horvath et al., (2019), Bayer et
al., (2019)), hence they reduce the duration of the calibration process by
several orders of magnitude. However, this kind of model replication is only a
speed-up, it is not able to alleviate systematic problems of the model if it
is not particularly appropriate to the actual market it aims to describe. One
may wonder if it might be more reasonable to apply the neural network model
directly to the observed data.

Recently, Bergeron et al., (2021) demonstrated that Variational Autoencoders
(VAEs) (Kingma and Welling, (2013)) are able to find a good low dimensional
representation of a large variety of volatility surfaces. They analyzed the FX
market where they combined data from several currency pairs during the
training process. They showed that a sufficiently trained VAE is not only able
to interpolate between points on the volatility surface but in the case of
sparse market data the decoder part of the VAE can efficiently reconstruct
those missing points and provide a reliable extended surface. Although the
results are very promising the trained neural network still remains a black
box model which is relatively hard to interpret. Unlike the traditional
approach, when we fit an analytic (e.g. SABR) model on the observed market and
get the dynamics of the underlying asset, in case of VAEs we do not have any
information about the stochastic processes under the hood. However, this
information is crucial when we need to price more complex, path-dependent
options (such as barriers). In our work we provide a hybrid solution for this
problem where we combine Deep Learning with a method called Weighted Monte
Carlo (Avellaneda et al., (2001)).

This paper is organised as follows: In Section 2 we provide the details of our
financial data and the conventional pricing methodology typically used by
market participants. In Section 3 we explain our methodology as well as the
implemented model architecture. We also show the results of our approach made
on the training and test sets. In Section 4 we discuss the key messages of our
study, and finally in Section 5 we summarize our work.

##  2 Financial data and traditional pricing methodology

A swaption (also known as a swap option) is an option to enter into an
interest rate swap. Swaptions are used on the market to manage interest rate
risk, either to hedge interest rate exposures or to trade on a directional
view. A swaption has two key properties: the expiry and the tenor. Expiry
defines the time interval in which one can enter into a swap, and tenor is the
duration of the underlying swap. Expiries can range from a few months to 30+
years, while tenors are in the range of 1-30 years. This shows that there is a
large variety of swaptions. To keep our analysis easy to interpret we chose a
specific type of swaption with a 2-year tenor. We used historical data of GBP
swaptions from ICAP (www.icap.com) with implied normal volatility values at 7
strike points around the at-the-money (ATM) point: ATM + [-50, -25, -12.5,
0.0, 12.5, 25, 50] basis points shifts, and at 7 expiries: 9M, 1Y, 18M, 2Y,
3Y, 4Y and 5Y. The observed time interval is between 2015 September 1 and 2020
October 21, hence we have volatility surface data for 1338 business days.

The traditional approach (Standard Model) to valuate swaptions uses a variant
of Black’s pricing formula originally developed to price options on futures
Black, (1976). For swaptions Black’s approach is justified if we work in the
annuity measure (Andersen and Piterbarg, (2010)). This uses the annuity factor

$A_{k,m}(t)=\sum_{n=k}^{k+m-1}P(t,T_{n+1})(T_{n+1}-T_{n})$ (1)

as the numeraire. Here $t$ is the time of observation, $T_{k}$ is the starting
date and $T_{k+m}$ is the maturity date of the (forward starting) swap.
$P(t,T)$ are zero coupon bond prices (discount factors) defined for the
underlying swap fixing dates.

The requirement of no arbitrage connects the forward swap rate $S_{k,m}(t)$,
i.e., the par rate of the forward starting swap, to the term structure
$P(t,T)$ observable at the time of pricing (Andersen and Piterbarg, (2010)):

$S_{k,m}(t)=\frac{P(t,T_{k})-P(t,T_{k+m})}{A_{k,m}(t)}.$ (2)

At time $t=0$, the forward swap rate $S_{0}=S_{k,m}(0)$ is known, but future
interest rates are uncertain, so for any $t>0$ the swap rate $S_{k,m}(t)$ is a
random process. Lack of arbitrage only dictates that in the annuity measure
the forward swap rate should be a martingale (see Andersen and Piterbarg,
(2010)):

$S_{k,m}(t_{1})=E_{t_{1}}^{A}[S_{k,m}(t_{2})],\quad\forall t1<t2.$ (3)

Any model that makes the swap rate a martingale is arbitrage free and can be
used to price other derivatives. In particular, in the annuity measure the
European swaption price can be expressed as

$V_{k,m}(0)=NA_{k,m}(0)E_{0}^{A}[S_{k,m}(T_{k})-K]$ (4)

where $N$ is the notional of the swap, $T_{k}$ is the maturity of the swaption
(the same as the starting time of the swap), $K$ is the strike rate, and
expectation should be calculated under the annuity measure. If we assume that
$S_{k,m}(T_{k})$ is lognormally distributed we arrive at the annuity scaled
Black formula typically used by practitioners. Note that lognormality can
arise from the assumption that the swap rate follows a geometric Brownian
motion, but to price a European option there is no need to know the dynamics
in detail, only the marginal distribution counts. Similarly, if we assume a
normal marginal (e.g., stemming from an arithmetic Brownian dynamics) we
arrive at a Bachelier type pricing formula (Choia et al., (2022)):

$\displaystyle c_{payer}$
$\displaystyle=A(0)(S_{0}-K)\Phi(d)+\sigma\sqrt{T-t}\phi(d)$ (5)
$\displaystyle c_{receiver}$
$\displaystyle=A(0)(K-S_{0})\Phi(-d)+\sigma\sqrt{T-t}\phi(d),$
$\displaystyle\text{where}\quad d=\frac{S_{0}-K}{\sigma\sqrt{T-t}}.$ (6)

containing a normal implied volatility $\sigma$, and using $\Phi()$ and
$\phi()$ as the probability distribution and density function of the standard
normal distribution, respectively. Note that payer/receiver are the swaption
nomenclature for call/put options. Since some interest rates went negative on
the market recently, the quoting convention for swaptions typically works with
normal implied volatilies. All vol surface data we obtained are quoted as
normal implied volatilities, and should be interpreted in the Bachelier
framework.

In the following we will use GBP swaptions with a 2-year tenor, but this is
only an illustration; our method is expected to work with other markets as
well. During our research we also investigated swaption products related to
other expiries and tenors as well as EUR swaptions, where we got similar
results to those presented here. Emphasis in this paper is on demonstrating
the applicability of the proposed pricing method, and the GBP swaption market
is used only as an example.

##  3 Methods and results

###  3.1 Building a Variational Autoencoder

After some investigation of the data we could clearly see that all volatility
surfaces have a relatively similar shape which indicates that the 7x7=49
points are correlated and can be expressed with much less parameters.
Similarly to the work of Bergeron et al., (2021) we applied a Variational
Autoencoder (VAE) to find the most compact latent representation behind the
vol surface manifold.

VAEs have been first introduced in Kingma and Welling, (2013) and they are
essentially the sophisticated version of Autoencoders. The input and output of
these models are the same, hence they belong to the unsupervised model family
meaning that they don’t need labeled data for training. A simple Autoencoder
neural network contains an encoder and a decoder model, where the encoder
projects the input parameters (here the volatility surface) into a lower
dimensional subspace, while the decoder tries to reproduce the original input
parameters as accurately as possible given the latent representation. Although
the reproduction efficiency of an Autoencoder can be high, there is no
guarantee that the latent space is smooth enough meaning that two similar
volatility surfaces are also close in the latent coordinates. This issue has
been solved by the Variational Autoencoders in which the encoder projects the
input on a distribution instead of a single point in the latent space and we
force the distribution to be as close as possible to a multivariate normal
distribution with a diagonal correlation matrix. In practice this discrepancy
is measured by calculating the Kullback-Leibler divergence (also known as
relative entropy) between the two distributions. This extra loss term in the
training process will help us to find the best latent representation of the
data in which the latent coordinates are independent from each other (they do
not correlate) and nearby latent data points represent similar volatility
surfaces.

The loss function we minimize is

$\mathcal{L}=\mathcal{L}_{rec}+\beta\mathcal{L}_{KL},$ (7)

where $\mathcal{L}_{rec}$ is the reconstruction loss between the real
($\mathbf{\sigma}$) and reconstructed ($\mathbf{\hat{\sigma}}$) volatility
vectors and $\mathcal{L}_{KL}$ is the Kullback-Leibler divergence. We used the
grid search approach to finetune the model hyperparameters over the
$N_{l}-N_{u}$ plane, where $N_{l}$ and $N_{u}$ denote the number of hidden
layers and the number of units in each layer, respectively, and searched for
an optimal model in the $N_{l}\times N_{u}\in[2,3,4]\times[5,10,15,20]$
domain. With regard to the model complexity and accuracy trade-off, we finally
selected the model whose key properties and training settings are summarized
in Table 1.

Table 1: Model architecture of the VAE as well as the training settings. Number of hidden layers (encoder) |  2
---|---
Number of hidden layers (decoder) |  2
Number of units in each layer |  10
Activation function of hidden layers |  elu
Batch size |  32
Initial learning rate |  0.001
Optimizer |  Adam
$\beta$ |  5e-5

As mentioned in Bayer et al., (2019) there are two main approaches in the
volatility surface reconstruction: the grid-based and the pointwise
approaches. While in the grid-based approach the decoder predicts the 49
points of the volatility surface directly, in the pointwise approach the
decoder is similar to a traditional pricer formula: it predicts a single
implied volatility value from the model parameters (latent coordinates) and
the $\Delta K$ strike offset and $T$ expiry. We chose the pointwise approach
since it has two main advantages: 1) Once the VAE has been trained
sufficiently, we can easily interpolate between the grid-points and can create
arbitrary fine resolution of the volatility surface; 2) We can easily
calculate the greeks using backpropagation, which is the most important
parameter in hedging applications. The general scheme of the pointwise
approach is shown in Figure 1.

Figure 1: Illustration of the data-flow through the Variational Autoencoder.
From left to right: the encoder reads in and transforms the 49 points of the
volatility surface into a 3 dimensional latent space; the decoder predicts a
given point of the volatility surface using the latent representation of the
volatility surface and the corresponding strike offset and expiry values.

###  3.2 Training the Variational Autoencoder

We split our data set into training and test sets using a split ratio of 7:3,
where the two sets were separated chronologically at the date of 2019 April 5
to avoid leakage of future information in the training process. Hence, the
blind test set contains normal as well as extreme (the pandemic crisis) market
environments. According to our data analysis we defined the 2019-04-9 -
2020-01-15 and 2020-01-16 - 2020-10-21 time intervals as the normal and
extreme periods of the test set, respectively. We also separated a validation
set – by randomly choosing 100 business days – from the training set to
continuously trigger the generalization performance during the training
process and save only those model weights of the VAE that produce the lowest
Mean Squared Error (MSE) in the validation set. This is a common practice to
avoid model overfitting. Similarly to the work of Bergeron et al., (2021) we
tried out the latent dimensions of 2, 3 and 4. In Table 2 we summarized the
model reconstruction performance (expressed in the Mean Squared Error) on the
training set with respect to each dimensionality of the latent space.

Table 2: Mean Squared Error between the reconstructed and real volatility surfaces of the training set using different number of latent dimensions. Number of latent dimensions |  2 |  3 |  4
---|---|---|---
Mean Squared Error |  1.068e-4 |  6.152e-5 |  4.339e-5

As we can observe increasing the dimensions from two to three reduced the MSE
by half. However, using 4 dimensions did not help significantly, therefore
considering the trade-off between model complexity and performance we used
three latent dimensions for our study. In Figure 2 we plotted the
reconstruction efficiency of the VAE regarding to the training set in the best
and worst cases considering the Mean Relative Absolute Error (MRAE) between
the real ($\boldsymbol{\sigma}_{real}$) and reconstructed
($\boldsymbol{\sigma}_{rec}$) implied volatility surfaces (see Eq. 8).

$MRAE=\frac{1}{N}\sum\limits_{i=1}^{N}\frac{\left|\sigma_{real,i}-\sigma_{rec,i}\right|}{\sigma_{real,i}}$
(8)

Figure 2: Illustration of the implied volatility surface reproduction
efficiency applying the fine tuned Variational Autoencoder network on the
training set. The top and bottom rows correspond to the best and worst
performing cases, respectively. From left to right we plotted the observed and
reconstructed implied volatility surfaces as well as the relative error
between them, respectively.

We can observe that in the best case the relative error is below 1% at almost
every point of the strike offset - expiry grid, and even in the worst case the
errors increase mostly to a few percents. This demonstrates that we found a
model that can well describe the large variety of volatility surfaces in the
training set with only three parameters. Furthermore, this model is able to
reliably interpolate between the observed strike offset - expiry grid points.
In Figure 3 we demonstrate this capability of the decoder network by
continuously increasing the resolution.

Figure 3: Illustration of the interpolation capability of the decoder network.
From left to right we increased the resolution from the original unequally
spaced 7x7 grid to equally spaced 10x10, 25x25 and 50x50 grids.

We can notice that the interpolated volatility surfaces remained smooth which
means that we did not overfit the decoder network on the data. Now, let us
investigate the effect of each latent dimension on the volatility surfaces,
where we used the finest, 50x50 resolution (see Figure 4).

Figure 4: Effect of each latent dimension on the volatility surface. Each row
represents a movement in the latent space from -2 to +2 along the
corresponding latent coordinate. The related position in the latent space can
be seen in the title of each plot. Color coding represents the difference
between the actual volatility surface and the volatility surface at the origin
[0,0,0] of the latent space.

We can see that each latent dimension has a clear meaning. Moving along the
first latent dimension increases the height of the volatility surface.
Modifying the second latent coordinate adjusts the skewness between the bottom
left and top right corner of the volatility surfaces. Finally, the last latent
coordinate has only a slight effect: it changes the slope of the volatility
smiles related to the shortest expiries.

Now, let us investigate the reproduction efficiency of the VAE on the blind
test set. Before diving into the details, we would like to highlight again
that the main purpose of training a Variational Autoencoder on historical data
is to catch the correlations between the volatility surface points and find a
low dimensional latent space from which the decoder neural network can
reproduce the observed volatility surfaces the best. After the VAE is trained
successfully, we can say that the network is good at interpolating between the
volatility surfaces of the training set. When we analyse the performance of
the VAE on unseen data we are essentially interested in its extrapolation
capability that depends mainly on how reliable the latent space is.

There are two ways to get the latent representation of a given volatility
surface: 1) We can use the encoder network to map directly to the latent
parameters (we will call this "direct encoding") or 2) we can fit the decoder
network using an optimization algorithm searching in the latent space (to be
called "calibration"). Direct encoding is certainly faster, but we have seen
cases when the encoder network failed to predict the correct latent
coordinates. In many cases we could achieve a smaller reproduction error using
the calibration approach. The observed effect may originate from the fact that
we intentionally introduce some noise in the training process and therefore
omitting the encoder will reduce the reconstruction error. However, deeper
analysis of this phenomenon is out of the scope of the current work.

To improve direct encoding we decided to boost training with synthetic data.
Given that our latent space is designed to be normally distributed, we
randomly sampled 10,000 points of the latent space from a Gaussian
distribution with zero mean and a variance of 4 to cover a broader space than
that of the training set (see Figure 5). We decoded these points into
synthetic volatility surfaces,111We found 44 unrealistic volatility surfaces
of the 10,000 that contained negative values. These data were removed from the
training set. and used this new synthetic data set to retrain the encoder part
of the VAE. During this fine tuning of the encoder, the decoder layers were
kept fixed.

Figure 5: Distribution of the synthetic, training and test set in the latent
space.

We determined the mean relative absolute error (MRAE) of the whole historical
data set for both the original and fine tuned versions of direct encoding and
also for the calibration approach. In this latter case we used the Trust
Region Reflective algorithm of Coleman and Li, (1996). In Figure 6 we plotted
the time evolution of MRAE related to the three approaches, and we determined
for reference the standard deviation of volatility surfaces from their mean
value within a 5-days long sliding window (gray thin line). This reference
quantity illustrates how large is the fluctuation in the market from day to
day. We can observe that all approaches produce almost the same accuracy,
0.66-0.79% in the training set. However, in the blind test set the MRAE curves
start to deviate and the lowest values are produced by the optimization method
(blue dashed line). Fine tuned direct mapping (green line) performs better at
several time regions (e.g.: 2020-01-15 - 2020-06-30) than its simply trained
version (red line), which means that the retraining process successfully
increased the accuracy. In the normal period of the data set the MRAE is about
0.86–1.23% while in the extreme period it increases to 2.14–3.30%, which is
still acceptable. We found that the trained model worked well for almost the
whole first year of the unseen test set producing the same MRAE as in the
training set, and only started to deviate at the beginning of the coronavirus
crisis in the middle of March 2020. This good performance is an indication
that the model is not overfitting on the training set, and therefore it is
able to extrapolate for a relatively long time horizon with only three latent
parameters. The increase of the error after one year is due to a drastic
change in the market. The low dimensional submanifold which represents the
physical market starts to move away from the active subspace of the model
learnt from the training set.

Figure 6: Time evolution of the mean relative absolute error on the training
and test sets. For reference we also show the standard deviation of volatility
surfaces within a 5-day long sliding window (gray thin line).

###  3.3 Weighted Monte Carlo

After successfully training the Variational Autoencoder on our swaption data
we get a generative model with three model parameters that can describe with
good accuracy the shape of every volatility surface that may typically arise
on the market. Now imagine that we currently observe the volatility surface of
swaptions today. Let’s investigate how we can use the trained VAE on this new
observation:

  * •

If we have all the points of the volatility surface, we can get the model
parameters by simply transmitting the data through the encoder.

  * •

If we have only a few points of the volatility surface we can still get the
related latent representation of the market by calibrating the decoder on the
data applying a given optimization method.

  * •

Now using the decoder we can price any vanilla option related to a given
($\Delta K$, $T$) pair.

  * •

We can also generate arbitrary number of synthetic, but realistic volatility
surfaces.

The shortcoming of our description is that the decoder network predicts the
implied volatility of the market directly, and does not give us any
information about the dynamics of the underlying asset, which we would need to
price any other non-vanilla derivatives on this market. Any information about
the dynamics minimally required to describe the volatility surface is there
implicitly in the latent coordinates. Therefore, to go beyond vanillas, we
need to build another "decoder" neural network that is able to translate the
latent representation of the volatility surface into the possible dynamics of
the underlying asset that the market participants expect. For this purpose we
used the so-called Weighted Monte Carlo method that was first introduced in
the field of option pricing by Avellaneda et al., (2001). In this framework we
start with an initial set of paths. We generated $\nu=20,000$ paths having a
5-year time horizon with 25,000 timesteps. In case of regular Monte Carlo we
take all paths with the same weight, $q_{i}=1/\nu,i\in[1,...,\nu]$, into
account. Given a security (in our case a swaption) with a path dependent
payoff $g_{i}$, we can calculate its discounted price $\Pi_{g}$ as

$\Pi_{g}=\frac{1}{\nu}\sum\limits_{i=1}^{\nu}g_{i}.$ (9)

In the Weighted Monte Carlo approach we assign a unique weight to each path
and modify them in a way that we can reproduce the observed market,

$\Pi_{g}=\sum\limits_{i=1}^{\nu}g_{i}p_{i}.$ (10)

This means however that we have a vast number of free parameters and therefore
we should introduce some regularization. Similar to Avellaneda et al., (2001)
we used the relative entropy measure, $D(\mathbf{p}/\mathbf{q})$, between the
$\mathbf{q}$ prior (here: the uniform distribution) and $\mathbf{p}$ posterior
distributions

$D(\mathbf{p}/\mathbf{q})=\sum\limits_{i=1}^{\nu}p_{i}\text{ln}\left(\frac{p_{i}}{q_{i}}\right).$
(11)

If we have observation of prices $c_{j}$ of $N$ securities, we can define a
$\mathbf{G}$ payoff matrix with $\nu$ rows and $N$ columns

$\mathbf{G}=\begin{pmatrix}g_{1,1}&g_{1,2}&\cdots&g_{1,N}\\\
g_{2,1}&g_{2,2}&\cdots&g_{2,N}\\\ \vdots&\vdots&\ddots&\vdots\\\
g_{\nu,1}&g_{\nu,2}&\cdots&g_{\nu,N}.\end{pmatrix}$

The calibration algorithm then goes for finding those $p_{i}$ values that
satisfy

$\sum\limits_{i=1}^{\nu}p_{i}g_{i,j}=c_{j},\quad j\in(1,2,...,N).$ (12)

We can also rewrite Eq. 12 in matrix form,

$\mathbf{p}\mathbf{G}=\mathbf{c}$ (13)

In the traditional Weighted Monte Carlo approach we solve this problem using
Lagrange multipliers $\lambda_{j}$, $j\in(1,2,...,N)$), and minimize the loss

$\mathcal{L}=D(\mathbf{p}/\mathbf{q})+\sum\limits_{j=1}^{N}\lambda_{j}\left(\sum\limits_{i=1}^{\nu}p_{i}g_{i,j}-c_{j}\right).$
(14)

According to Avellaneda et al., (2001) the weights can be expressed with the
Lagrange multipliers as

$p_{i}=\frac{1}{Z(\mathbf{\lambda})}\exp\left(\sum\limits_{j=1}^{N}g_{i,j}\lambda_{j}\right),\quad\text{where}\quad
Z(\mathbf{\lambda})=\frac{1}{\nu}\sum\limits_{i=1}^{\nu}\exp\left(\sum\limits_{j=1}^{N}g_{i,j}\lambda_{j}\right).$
(15)

In this way the weight vector is a function $\mathbf{p}=f(\mathbf{\lambda})$
of as many parameters as observations we have on the market. In our case, we
have 7x7=49 data points on the vol surface, which means that we can calculate
the price of at least 49 securities. However, using the Variational
Autoencoder we found that this 49-dimensional space is degenerate and the
market can be described with good accuracy by only three parameters. Therefore
we can try to build a new Weight Decoder neural network, $\Phi_{WD}$, which is
able to translate the latent representation of the market into the weights of
possible paths,

$\mathbf{p}=\Phi_{WD}(\mathbf{z}).$ (16)

In Figure 7 we illustrate the data processing from the latent space to the
vanilla price surface using the Weight Decoder. Note that we work with the
option price surface, and not the vol surface here, as the Monte Carlo method
gives prices, not IVs directly.

Figure 7: Illustration of the data-flow through the Weight Decoder neural
network. From left to right: the weight decoder reads in the latent
coordinates and transforms them into the weights of different paths; we
calculate then the prices of securities using Weighted Monte Carlo. During the
training phase we try to reproduce the prices of securities (that we
previously determined from the observed implied volatility values) as
accurately as possible.

First we apply the Weight Decoder on the latent space to get the weights
related to the paths. Then we calculate the prices of securities under
consideration using Weighted Monte Carlo. Since we do not have an explicit
formula from prices of securities to implied volatilities, we generated call
($c_{call}$) and put ($c_{put}$) options related to the observed volatility
surface using the Bachelier formulae, Eqs. 5 and 6. We ran the training
process to reproduce these prices as accurately as possible.

It is worth mentioning that we use the call $\hat{c}_{call}$ and put prices
$\hat{c}_{put}$ related to the reconstructed implied volatilities of the VAE
decoder network. This is a convenient choice since this way we train the
weight decoder on smooth price surfaces that do not contain noise anymore. The
key properties of the model architecture is summarized in Table 3.

Table 3: Model architecture of the weight decoder as well as the training settings. Number of hidden layers |  2
---|---
Number of units in each layer |  20/49
Activation function of hidden layers |  elu
Activation function of the readout layer |  softmax
Batch size |  32
Initial learning rate |  0.01
Optimizer |  Adam
$\gamma$ |  1e-8

###  3.4 Put-call parity

To avoid arbitrage opportunities it is necessary to ensure that the put-call
parity holds for all predicted prices:

$K-S_{0}+c_{call}-c_{put}=\Delta K+c_{call}-c_{put}=0$ (17)

Instead of introducing an extra loss term for the put-call parity we embedded
this relation as a constraint into the put price estimation from call prices:

$\hat{\hat{c}}_{put}=\Delta K+\hat{\hat{c}}_{call}$ (18)

The loss of call and put price reconstruction is then considered during the
training, where we denoted the reconstructed call and put prices by
$\hat{\hat{c}}_{call}$ and $\hat{\hat{c}}_{put}$, respectively. The final form
of the loss function is

$\mathcal{L}=\mathcal{L}_{rec,call}+\mathcal{L}_{rec,put}+\gamma
D(\mathbf{p}/\mathbf{q}),$ (19)

where $\mathcal{L}_{rec,call}$ and $\mathcal{L}_{rec,put}$ are the
reconstruction losses of the prices of put and call options,
$D(\mathbf{p}/\mathbf{q})$ is the relative entropy between the predicted
weights and the uniform distribution, and $\gamma$ is a weighting constant.

###  3.5 Martingale condition

The main purpose of the calibration method described above is to reproduce the
observed option prices. However, to price derivatives without arbitrage, the
paths have to satisfy the martingale condition Eq. 3. Although the original
paths have been generated from a martingale process, the re-weighting process
can break the martingale condition as was found in Elices and Giménez, (2011).
They argue that in case of calibrating on several expiries with steep smiles,
the resulting underlying process may fail to satisfy the martingale condition.

Note that including the put-call parity breaking penalty per Section 3.4 is
already a way to curb deviations from the martingale condition. However to
force this requirement more rigorously in a discrete Monte Carlo environment
we would need more control. Consider the subset of paths that go through a
narrow window at time $t$ around the point $S=s$, denoted
$W_{t}(s)=[s-\Delta,s+\Delta]$, where $\Delta=0.001$ by our choice. In Figure
8 we illustrate these paths going through two different windows.

Figure 8: Illustration of the martingale condition in the Monte Carlo
environment. The blue and red paths correspond to two different windows,
$W_{0.75}(-0.05)$ and $W_{0.75}(0.05)$ at time $t_{1}=0.75$. The grey paths
represent the full path set and the green vertical lines denote the range
within we defined the windows.

To be a martingale process, for every window $W=W_{t_{1}}(s)$ and later time
$t_{2}$ the relationship

$E_{t_{1}}^{A}[S_{t_{1}}|S_{t_{1}}\in W]=E_{t_{1}}^{A}[S_{t_{2}}|S_{t_{1}}\in
W]$ (20)

should hold. This can be rewritten in a path discrete form

$\sum\limits_{j|S_{j,t_{1}}\in
W}\left(S_{j,t_{2}}-S_{j,t_{1}}\right)\frac{p_{j}}{\sum\limits_{i|S_{i,t_{1}}\in
W}p_{i}}=0.$ (21)

One way to introduce the martingale condition into the training process is to
create additional fictitious derivatives with zero price ($c_{f}=0$) that have
the payoff function

$g_{\text{mart},j}=\left(S_{j,t_{2}}-S_{j,t_{1}}\right),\quad\text{where}\quad
S_{j,t_{1}}\in W$ (22)

for a set of windows $W_{t_{i}}(s)$ scattered around time and rate space.
Considering only the consecutive expiries (as in Elices and Giménez, (2011))
using windows at 20 different positions between -0.01 and 0.01222Note that the
range of strikes is between -0.005 and 0.005, hence we cover the important
region. we have 6x20=120 additional derivatives. During the training process
we force the model to reproduce the zero prices of these derivatives by adding
the same loss term as for the call and put options.

During our analysis we found that this extra loss term is computationally very
expensive, increasing the time of a single training step by an order and the
training loss was still fluctuating, the model did not converge. Therefore, we
decided to train the weight decoder without these loss terms and investigate
only the deviation of our emerging measure from the martingale property. In
Figure 9 we plotted the time evolution of the mean martingale loss relative to
the support of window positions (0.01-(-0.01) = 0.02). In this way we give an
intuition about the severity of the observed deviation from the martingale
condition.

Figure 9: Time evolution of the observed martingale loss (red line) relative
to the range of window positions ([-0.01,0.01]). For reference we plotted the
level of martingale loss (green dashed line) related to the noise originated
from the path generation process.

We can see that the discrepancy is relatively small and stays in the 2-6%
range for the whole data set. Note, however, that part of this noise comes
from the Monte Carlo method itself. To estimate this contribution we
regenerated 50 times the used set of 20,000 paths and calculated the standard
deviation of the martingale loss for the unweighted paths. The relative value
of this amplitude was found to be about 1% (see green dashed line in Figure
9). This means that a non-negligible part of the martingale loss comes from
the sampling noise of the Monte Carlo method. We can conclude that using the
put-call parity relation during the training process is sufficient to keep the
martingale loss at a sufficient level. The addition of a more complex set of
constraints makes the approach difficult to manage and may not have tangible
benefits. One should however keep the option of these extra controls in mind
if dealing with securities particularly sensitive to the martingale condition.

###  3.6 Training the weight decoder

At this stage we have a latent representation of the swaption market that
contains all necessary information. Since the observable market essentially
shows us the expectations of the market participants, the derived latent space
must be in a direct relation with the possible paths of the underlying asset.
Therefore we built a new neural network, the weight decoder that is able to
translate the latent coordinates into the weights of the paths to be used
later in the Weighted Monte Carlo method. We trained this model on the same
training set as was used in the fine tuning process of the encoder network. To
continuously trigger the performance of the network we chose the original
training set (the period of 2015-09-01 - 2019-04-09) as the validation set. In
Figure 10 we plot the reconstruction efficiency of the weight decoder, where
we made a comparison with the best performing calibration method (see Fig. 6).
The reconstructed implied volatilities have been iteratively calculated from
the calculated call price surface using the Bachelier formula.

Figure 10: Time evolution of the Mean Relative Absolute Error of the weight
decoder related to the training and test sets (orange line). We compared our
results to the best performing calibration method (blue dashed line). For
reference we calculated the relative standard deviation of volatility surfaces
within a 5-days long sliding window (gray thin line).

We can observe that the reproduction efficiency achieved by the weight decoder
is very close to the calibration method (1.06% and 1.63% for training and test
sets, respectively). This means that the weight decoder predicts the correct
weights from the latent space, that can be assigned to the paths to reproduce
the observed market accurately.

In Figure 11 we plotted the effect of each latent dimension on the risk-
neutral densities using the trained weight decoder where we analyzed the same
positions in the latent space as in Figure 4. We can observe that the first
latent dimension has the most significant effect, that modifies the volatility
of the paths. The other two coordinates have much less influence on the risk
neutral densities as we expected from Figure 4.

Figure 11: Effect of each latent dimension on the risk neutral densities. Each
row represents a movement in the latent space from -2 to +2 along the
corresponding latent coordinate. The related position in the latent space can
be seen in the title of each plot. All of the risk-neutral densities are
scaled the same for comparison purposes.

###  3.7 Pricing a barrier call option

In the following we will price a path-dependent option, namely an up barrier
call option with an expiry of five years. The payoff value of this product is
zero in case of paths that exceed the predefined barrier before the expiry and
is equal to the price of a regular call option otherwise. We chose the
swaption market of 2016-03-30, where we got an accurate fit using the Weight
Decoder. For comparison we fitted the analytic $\beta=0$ SABR formula as well
on the volatility smile related to the 5Y expiry333The calibration of
stochastic models on the whole volatility surface is often impossible and
therefore the usual way for calibrating e.g. the SABR model is to fit the
approximation formula to each volatility smile, individually. to get the
$\alpha$, $\rho$ and $\nu$ model parameters driving the two stochastic
differential equations

$\displaystyle dS_{t}$ $\displaystyle=\sigma_{t}S_{t}^{\beta}dW_{t}^{1}$ (23)
$\displaystyle d\sigma_{t}$ $\displaystyle=\nu\sigma_{t}dW_{t}^{2},$
$\displaystyle\text{where}\quad dW_{t}^{1}dW_{t}^{2}=\rho
dt\quad\text{and}\quad\alpha=\sigma_{0}$ (24)

In Figure 12 we plotted the real and reconstructed volatility surfaces using
the Weight Decoder neural network as well as the SABR model.

Figure 12: Illustration of the implied volatility surface reproduction
efficiency applying the weight decoder as well as the SABR model. From left to
right we plotted the observed and reconstructed implied volatility surfaces as
well as the relative error between them, respectively. Note, that the SABR
model was calibrated on the last expiry.

We can see that the Weight Decoder is able to accurately reproduce the
observed market with a mean relative absolute error of 0.41%, while the SABR
model fails to explain the volatilities of expiries shorter than 5 years, the
discrepancy between the fitted and observed values can be more than 30%. To
illustrate the difference between the predicted dynamics of the underlying, we
plot the related paths as well as the risk neutral distributions at each
expiry in Figure 13. We can notice that the weighted paths on the left hand
side follow a slightly asymmetric probability density function, similar to a
log-normal distribution. The paths related to the SABR model on the right hand
side contain however more extreme realizations that can be seen in the wider
probability densities as well.

Figure 13: Top row: Stock price dynamics predicted by the Bachelier and the
SABR models. In case of the Bachelier paths we marked the estimated weights of
paths using the alpha value of the lines (more possible paths are less
transparent). Bottom row: Risk neutral densities at each expiry.

Now, we will set the strike price of the barrier option to ATM and
continuously modify the barrier from 0.01 to 0.1. The calculated option prices
based on the two path sets can be seen in Figure 14. We can see that as long
as the barrier level is under 0.02 the prices are in agreement with each
other. However, there is a significant difference in the prices in the range
of [0.02, 0.08]. This implies that for such types of barrier options it is
crucial to find the correct dynamics of the underlying asset. In our case the
weighted Bachelier paths are able to calibrate to the term structure of the
market, therefore the calculated barrier option prices are more reliable than
using the single term SABR model.

Figure 14: Price of ATM barrier call option as a function of the barrier
level. Prices have been calculated using i) weighted Bachelier and ii)
unweighted SABR paths.

##  4 Discussion

During our analysis we made two main findings: 1) We have confirmed that
Variational Autoencoders are extremely effective in finding the correlations
of the swaption market and hence they can massively reduce the number of
degrees of freedom; 2) We demonstrated how to extract physically meaningful
information from the derived latent space by translating the latent
coordinates into weights of paths. We believe that our methodology can be used
more generally. At many scientific fields we often face the problem of
describing a specific phenomenon of the world by finding a mathematical model.
Usually, as we observe more and more data the created analytic model is less
and less powerful in explaining the new observations and therefore becomes
more and more complex to keep its original performance. Contrarily, we claim
that finding the minimally needed number of free parameters for the first step
of the analysis is very important and Variational Autoencoders are great tools
in doing that. Creating an interpretable mathematical model with some model
parameters makes sense only after the previous step. We can use then neural
networks to find the non-trivial relations between the latent space and our
model parameters. This final process essentially means to extract the
physically meaningful information from the latent space.

##  5 Conclusions

Our aim in this work was to find low dimensional representations of the
swaption market. We built a Variational Autoencoder that was able to represent
the market for more than five years using only three model parameters. We
proposed a method to translate this static representation into a consistent
dynamic model under the Weighted Monte Carlo framework initially staring from
simple Brownian paths. Our Weight Decoder neural network was able to map the
latent parameter space into the required path weights directly with good
accuracy. The resulting (weighted) Monte Carlo measure can be applied to price
non-vanilla derivatives such as barrier options. We compared prices to those
resulting from a SABR model approach and demonstrated that the lack of
matching the term structure of volatilities in the latter can lead to huge
pricing differences. Clearly, finding a consistent dynamical model is crucial
for adequate pricing of path dependent products. Our codes are publicly
available at https://github.com/skunsagimate/DeepWMC.git.

##  6 Acknowledgements

Prepared with the support of the Doctoral Student Scholarship Program of the
Cooperative Doctoral Program and the MILAB Artificial Intelligence National
Laboratory Program of the Ministry of Innovation and Technology from the
National Research, Development and Innovation Fund.

##  7 Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible
for the content and writing of the paper.

## References

  * Aguilar and West, (2000) Aguilar, O. and West, M. (2000).  Bayesian dynamic factor models and portfolio allocation.  Journal of Business & Economic Statistics, 18(3):338–57.
  * Andersen and Piterbarg, (2010) Andersen, L. B. G. and Piterbarg, V. V. (2010).  Interest Rate Modeling, Volume 1.  Atlantic Financial Press.
  * Avellaneda et al., (2001) Avellaneda, M., Friedman, C., Buff, R., and Branchamp, N. (2001).  Weighted Monte-Carlo: A new technique for calibrating asset-pricing models.  International Journal of Theoretical and Applied Finance, pages 91–119.
  * Bayer et al., (2019) Bayer, C., Horvath, B., Muguruza, A., Stemper, B., and Tomas, M. (2019).  On deep calibration of (rough) stochastic volatility models.  arXiv e-prints, page arXiv:1908.08806.
  * Bergeron et al., (2021) Bergeron, M., Fung, N., Hull, J., and Poulos, Z. (2021).  Variational Autoencoders: A Hands-Off Approach to Volatility.  arXiv e-prints, page arXiv:2102.03945.
  * Black, (1976) Black, F. (1976).  The pricing of commodity contracts.  Journal of Financial Economics, 3(1-2):167–179.
  * Black and Scholes, (1973) Black, F. and Scholes, M. (1973).  The pricing of options and corporate liabilities.  Journal of Political Economy, 81(3):637–654.
  * Choia et al., (2022) Choia, J., Kwakb, M., Teec, C. W., and Wang, Y. (2022).  A black–scholes user’s guide to the bachelier model.  arXiv e-prints, page arXiv:2104.08686.
  * Coleman and Li, (1996) Coleman, T. F. and Li, Y. (1996).  An interior trust region approach for nonlinear minimization subject to bounds.  SIAM Journal on Optimization, 6(2):418–445.
  * Derman et al., (1996) Derman, E., Kani, I., and Zou, J. Z. (1996).  The local volatility surface: Unlocking the information in index option prices.  Financial Analysts Journal, 52(4):25–36.
  * Dumas et al., (1998) Dumas, B., Fleming, J., and Whaley, R. E. (1998).  Implied volatility functions: Empirical tests.  The Journal of Finance, 53(6):2059–2106.
  * Dupire et al., (1994) Dupire, B., (see Black, T. B. M., and Options, G. (1994).  Pricing with a smile.  Risk Magazine, pages 18–20.
  * Elices and Giménez, (2011) Elices, A. and Giménez, E. (2011).  Weighted Monte Carlo: Calibrating the Smile and Preserving Martingale Condition.  arXiv e-prints, page arXiv:1102.3541.
  * Gaspar et al., (2020) Gaspar, R. M., Lopes, S. D., and Sequeira, B. (2020).  Neural network pricing of american put options.  Risks, 8(3).
  * Hagan et al., (2002) Hagan, P., Kumar, D., Lesniewski, A., and Woodward, D. (2002).  Managing smile risk.  Wilmott Magazine, 1:84–108.
  * Heston, (1993) Heston, S. L. (1993).  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  The Review of Financial Studies, 6(2):327–343.
  * Horvath et al., (2019) Horvath, B., Muguruza, A., and Tomas, M. (2019).  Deep Learning Volatility.  arXiv e-prints, page arXiv:1901.09647.
  * Kim et al., (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 Welling, (2013) Kingma, D. P. and Welling, M. (2013).  Auto-Encoding Variational Bayes.  arXiv e-prints, page arXiv:1312.6114.
  * Liu et al., (2019) Liu, S., Oosterlee, C. W., and Bohte, S. M. (2019).  Pricing options and computing implied volatilities using neural networks.  Risks, 7(1).
  * Pagnottoni, (2019) Pagnottoni, P. (2019).  Neural network models for bitcoin option pricing.  Frontiers in Artificial Intelligence, 2.
  * Tripathy and Bilionis, (2018) Tripathy, R. K. and Bilionis, I. (2018).  Deep uq: Learning deep neural network surrogate models for high dimensional uncertainty quantification.  Journal of Computational Physics, 375:565–588.