# Financial option valuation by unsupervised learning with artificial neural
networks

Beatriz Salvador1, Cornelis W. Oosterlee1,2, Remco van der Meer1,2 1 CWI –
Centrum Wiskunde & Informatica, Amsterdam, the Netherlands
2 DIAM, Delft University of Technology, Delft, the Netherlands

###### Abstract

Artificial neural networks (ANNs) have recently also been applied to solve
partial differential equations (PDEs). In this work, the classical problem of
pricing European and American financial options, based on the corresponding
PDE formulations, is studied. Instead of using numerical techniques based on
finite element or difference methods, we address the problem using ANNs in the
context of unsupervised learning. As a result, the ANN learns the option
values for all possible underlying stock values at future time points, based
on the minimization of a suitable loss function. For the European option, we
solve the linear Black-Scholes equation, whereas for the American option, we
solve the linear complementarity problem formulation. Two-asset exotic option
values are also computed, since ANNs enable the accurate valuation of high-
dimensional options. The resulting errors of the ANN approach are assessed by
comparing to the analytic option values or to numerical reference solutions
(for American options, computed by finite elements).

###### keywords:

(non)linear PDEs, Black-Scholes model, artificial neural network, loss
function, multi-asset options

##  1 Introduction

The interest in machine learning techniques, due to the remarkable successes
in different application areas, is growing exponentially. Impressive results
have been achieved in image recognition or natural language processing
problems, among others. The availability of large data sets and powerful
compute units has brought the broad field of data science to a next level.
ANNs are learning systems based on a collection of artificial neurons that
constitute a connected network [26]. Such systems “learn” to perform tasks,
generally without being programmed with task-specific rules. The neurons are
organized in multiple layers; The input layer receives external data, the
output layer produces the final result. The layers in between input and output
are the so-called hidden layers [25]. Many different financial problems have
also been addressed with machine learning, like stock price prediction, where
ANNs are trained to detect patterns in historical data sets to predict future
trends [12, 20], or bond rating predictions, see [6, 17, 21].

Motivated by the universal approximation theorems [3, 4], nowadays ANNs are
also being used to approximate solutions to ordinary differential equations
(ODEs) or partial differential equations (PDEs) [1,5-8]. We will contribute to
this field by solving some PDEs that appear in computational finance
applications with ANNs, following the unsupervised learning methodology
introduced by [19] and refined in [16]. The resulting ANN-based methods do not
require a discretization of the differential equation, and mesh generation is
therefore not required.

The financial application on which we focus is the valuation of financial
derivatives with PDEs. Generally, we can distinguish between supervised and
unsupervised machine learning techniques. Research so far has mainly focused
on supervised machine learning, i.e. given input variables $x$ and labeled
output variables $y$, the ANN is employed to learn the mapping function from
the input to the output. The goal is then to approximate the mapping function
accurately, so that for new input data $x^{\prime}$, the corresponding output
$y^{\prime}$ is well approximated. Such ANN methodology usually consists of
two phases. During the training phase, the ANN should learn the PDE solver
with input parameters and output. This (off-line) phase usually takes
substantial computing time. In the testing phase, the trained model is used to
very rapidly approximate solutions for other parameter sets. In [8], the
authors showed that ANNs efficiently approximate the solution to the Black-
Scholes equation. In [14], option values as well as the corresponding implied
volatilities were directly computed with one neural network in a supervised
learning approach. The authors in [1] examined whether an ANN could derive
option pricing formulas based on market prices. ANN studies for American
options are also found, like in [11], and in [22], where the option was
formulated as a free boundary problem. In [13] the American option implied
volatility and implied dividend were assessed with the help of ANNs.

The goal of the current work is to solve the financial PDEs by applying
unsupervised machine learning techniques. In such a case, only the inputs of
the network are known, and based on a suitable loss function that needs to be
minimized, the ANN should “converge” to the solution of the PDE problem. The
ANN should learn solutions that satisfy constraints that are imposed by the
PDE and the boundary conditions, without using any information about the true
solution. These constraints are typically formulated as soft constraints, that
are satisfied by minimizing some loss function. The potential advantage of
applying ANNs to address PDE problems, instead of using classical numerical
methods, is found in the problem’s dimensionality. An ANN-based methodology
does not suffer much from the curse of dimensionality. The authors of [1, 6,
7] provide evidence that for the well-known Poisson and Burgers equations,
these unsupervised learning methods yield accurate results. The authors in
[16] extended the class of PDE solutions that may be approximated by these
unsupervised learning methods, by translating the PDEs to a suitably weighted
minimization problem for the ANNs to solve. Moreover, in [4, 5] American
options were formulated as optimal stopping problems, where optimal stopping
decisions were learned and so-called ANN regression was used to estimate the
continuation values. This is an example of the unsupervised learning approach
to solve a specific formulation of options with early-exercise features.

We will price European and American options modeled by the Black–Scholes PDE
and look for solutions for all future time points and stock values. So, linear
and nonlinear partial differential equations need to be solved. We will solve
European and American option problems based on one and two underlying assets,
as the methodology is easily extended to solving multi-asset options. For the
European problems, the accuracy of the network can be measured as we have the
analytic Black-Scholes solution as a reference. American options will be
formulated as linear complementarity problems. Since an analytic solution is
not known in this case, the reference solutions are obtained by finite element
computations on fine meshes.

This paper is organized as follows. In Section 2, the methodology to train the
neural network is introduced. In Section 3, the financial PDE problems are
formulated, for the linear and the nonlinear case. Numerical results, ANN
convergence and solution accuracy, are presented in Section 4. Finally,
Section 5 concludes.

##  2 Artificial Neural Networks Solving PDEs

In this section, we introduce the methodology following [16] to solve linear
and nonlinear time-dependent PDEs by ANNs. With this aim, we write a general
PDE problem as follows:

$\displaystyle\mathcal{N}_{I}(v(t,x))$ $\displaystyle=0,\quad
x\in\widetilde{\Omega},\,t\in[0,T],$ $\displaystyle\mathcal{N}_{B}(v(t,x))$
$\displaystyle=0\quad\text{on}\,\,\partial\widetilde{\Omega},$ (1)
$\displaystyle\mathcal{N}_{0}(v(t^{\ast},x))$ $\displaystyle=0\quad
x\in\widetilde{\Omega}\,\,\text{and}\,\,t^{\ast}=0\,\,\text{or}\,\,t^{\ast}=T,$

where $v(t,x)$ denotes the solution of the PDE, $\mathcal{N}_{I}(\cdot)$ is a
linear or nonlinear time-dependent differential operator,
$\mathcal{N}_{B}(\cdot)$ is a boundary operator, $\mathcal{N}_{0}(\cdot)$ is
an initial or final time operator, $\widetilde{\Omega}$ is a subset of
$\mathbb{R}^{D}$ and $\partial\widetilde{\Omega}$ denotes the boundary on the
domain $\widetilde{\Omega}$.

As mentioned in the introduction, we will compute European and American option
values for one and two underlying assets by unsupervised learning. The goal is
to obtain $\hat{v}(t,x)$ by minimizing a suitable loss function $L(v)$ over
the space of $k$-times differentiable functions, where $k$ depends on the
order of the derivatives in the PDE, i.e

$arg\min_{v\in\mathcal{C}^{k}}L(v)=\hat{v}\,,$

where we denote by $\hat{v}(t,x)$ the true solution of the PDE.

Results are available that establish a relation between the value of the loss
function and the accuracy of the approximated solution. A general expression
for the loss function, defined in terms of the $L^{p}$ norm, including a
weighting, is defined as follows [19, 16]:

$\displaystyle L(v)$
$\displaystyle=\lambda\int_{\Omega}\mid\mathcal{N}_{I}(v(t,x))\mid^{p}d\Omega$
(2)
$\displaystyle+(1-\lambda)\int_{\partial\Omega}\left(\mid\mathcal{N}_{B}(v(t,x))\mid^{p}+\mid\mathcal{N}_{0}(v(t,x))\mid^{p}\right)d\gamma,$
(3)

where $\Omega=\widetilde{\Omega}\times[0,T]$, $\partial\Omega$ the boundary of
$\Omega$ and

$\displaystyle\mathcal{N}_{I}(v(t,x))\equiv
N(v(t,x))-F(t,x)\quad\text{in}\,\,\Omega\,,$
$\displaystyle\mathcal{N}_{B}(v(t,x))\equiv
B(v(t,x))-G(t,x)\quad\text{on}\,\,\partial\widetilde{\Omega}\,,$
$\displaystyle\mathcal{N}_{0}(v(t^{\ast},x))\equiv
H(x)-v(t^{\ast},x)\quad\text{in}\,\,\widetilde{\Omega}\times
t^{\ast},\,\,\text{with}\,\,t^{\ast}=0\,\,\text{or}\,\,t^{\ast}=T.$

The integrals of the loss function are labeled as:

$L_{I}(v)\equiv\int_{\Omega}\mid\mathcal{N}_{I}(v(t,x))\mid^{p}d\Omega,$

and

$L_{B}(v)\equiv\int_{\partial\Omega}\left(\mid\mathcal{N}_{B}(v(t,x))\mid^{p}+\mid\mathcal{N}_{0}(v(t,x))\mid^{p}\right)d\gamma,$

which are denoted as the interior and the boundary loss functions,
respectively.

Financial options with early-exercise features give rise to free boundary PDE
problems. Free boundary problems are well-known and often appearing in a
variety of engineering problems. We recall some classical formulations of the
free boundary problems that we encounter here:

  * 1.

An optimal stopping time problem,

  * 2.

A linear complementarity problem (LCP),

  * 3.

A parabolic variational inequality,

  * 4.

A penalty problem.

We will focus on the reformulation of the free boundary problem as an LCP, and
aim to solve this formulation by ANNs and unsupervised learning. The generic
LCP formulation reads,

$\displaystyle\mathcal{N}_{I}(v(t,x))\cdot\mathcal{N}_{0}(v(t,x))$
$\displaystyle=0,\quad x\in\widetilde{\Omega},t\in[0,T].$ (4)

or, equivalently,

$\displaystyle\max(\mathcal{N}_{0}(v(t,x)),\mathcal{N}_{I}(v(t,x)))$
$\displaystyle=0,\quad x\in\widetilde{\Omega},\,t\in[0,T],$
$\displaystyle\mathcal{N}_{B}(v(t,x))$
$\displaystyle=0,\quad\text{on}\,\,\partial\widetilde{\Omega},$
$\displaystyle\mathcal{N}_{0}(v(t^{\ast},x))$ $\displaystyle=0,\quad
x\in\widetilde{\Omega}\,\,\text{and}\,\,t^{\ast}=0\,\,\text{or}\,\,t^{\ast}=T.$

Our expression for the loss function, to solve the linear complementarity
problem, is as follows:

$\displaystyle L(v)$
$\displaystyle=\lambda\int_{\Omega}\mid\max(\mathcal{N}_{0}(t,x,v),\mathcal{N}_{I}(t,x,v))\mid^{p}d\Omega$
$\displaystyle+(1-\lambda)\int_{\partial\Omega}\left(\mid\mathcal{N}_{B}(t,x,v)\mid^{p}+\mid\mathcal{N}_{0}(t,x,v)\mid^{p}\right)d\gamma\,.$
(5)

As an alternative loss function for the LCP, a variance normalization loss
function has also been considered [16], which is defined as:

$\displaystyle L(v)$ $\displaystyle=$
$\displaystyle\frac{\int_{\Omega}\mid\max(\mathcal{N}_{0}(t,x,v),\mathcal{N}_{I}(t,x,v))\mid^{p}dx}{\int_{\Omega}(\max(\mid\mathcal{N}_{0}(t,x,v)\mid,\hat{\mathcal{N}}_{I}(t,x,v)))^{p}dx}$
(6) $\displaystyle+$
$\displaystyle\frac{\int_{\partial\Omega}\left(\mid\mathcal{N}_{B}(t,x,v)\mid^{p}+\mid\mathcal{N}_{0}(t,x,v)\mid^{p}\right)d\gamma}{\int_{\partial\Omega}\mid
v(t,x)-\bar{v}\mid^{p}d\gamma},$

where $\hat{\mathcal{N}}_{I}$ is defined as $\mathcal{N}_{I}$ but considering
each term in absolute value and $\bar{v}$ is the mean of $v$ over the
corresponding domain.

The parameter $\lambda\in(0,1)$ in the loss functions represents the relative
importance of the interior and boundary functions in the minimization process.
The choice of such value can be addressed in different ways, see [19, 16]. In
this work, the loss weight is, in most of the tests, set equal to
$\lambda=0.5$. It was found in [16] that this choice works very well for PDE
problems with smooth, non-oscillatory solutions (as we also encounter them in
the option valuation problems under consideration). For some linear
complementarity problems, we will compare the basic choice with the variance
normalization loss function. In addition, for some other cases, we will
compute the loss function considering a so-called optimal loss weight (as in
[16]).

Based on the loss function, the ANN has been trained with the Broyden-
Fletcher-Goldfarb-Shanno optimization (BFGS). This is a quasi-Newton method
which employs an approximate Hessian matrix. Particularly, we use the L-BFGS
algorithm to optimize the vector $\theta$, which contains all parameters
defining the neural network. The activation function used in the ANN is the
hyperbolic tangent function $\tanh(x)$, however, other choices of the
activation function can also be used, like the sigmoid function (resulting in
very similar results in this work). We will work with relatively small neural
networks formed by four hidden layers with $20$ neurons each for the European
and American options. Increasing the number of layers did not improve the
accuracy of the solution significantly for these particular problems. Finally,
the integral terms in the loss function are approximated by Monte Carlo
techniques.

##  3 Financial derivative pricing partial differential equations

In this section the option pricing partial differential equation problems are
presented. We briefly introduce the models.

###  3.1 European options, one underlying asset

The reference option pricing PDE for the valuation of a plain vanilla
European, put or call, option is the Black-Scholes equation. The underlying
asset $S_{t}$ is assumed to pay a constant dividend yield $\delta$, and
follows the geometric Brownian motion:

$dS_{t}=(\mu-\delta)S_{t}dt+\sigma S_{t}dW^{P}_{t}\,,$ (7)

where $W^{P}_{t}$ is a Brownian motion. The drift term $\mu$, the risk-free
interest rate, $r$, and the asset volatility, $\sigma$, are known functions.
Assuming there are no arbitrage opportunities, the European option value
follows from the Black–Scholes equation,

$\displaystyle\begin{cases}\mathcal{L}(v)=\partial_{t}v+\mathcal{A}v-rv=0\,,\quad\quad
S\in\widetilde{\Omega}\,\,,t\in[0,T)\,,\\\ v(T,S)=H(S)\,,\end{cases}$ (8)

where operator $\mathcal{A}$ is defined as,

$\mathcal{A}v\equiv\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2}v}{\partial
S^{2}}+(r-\delta)S\frac{\partial v}{\partial S}$ (9)

and function $H$ denotes the option’s payoff, which is given by:

$\displaystyle\begin{cases}(K-S)^{+}\quad\quad\text{for a put option}\\\
(S-K)^{+}\quad\quad\text{for a call option}\,,\end{cases}$ (10)

with $K$ the strike price in the option contract.

In order to apply numerical methods to solve the PDE, a bounded domain should
be considered and a proper set of boundary conditions should be imposed. We
assume a domain large enough being $[0,S_{\infty}]$, with $S_{\infty}$ four
times the strike $K$. Depending on the kind of option, call $v_{c}$ or put
$v_{p}$, the problem (8) is subject to the conditions:

$\displaystyle\begin{cases}v_{c}(t,0)=0\\\
v_{c}(t,S_{\max})=S_{\max}-Ke^{-r(T-t)}\,,\end{cases}\quad\quad\begin{cases}v_{p}(t,0)=Ke^{-r(T-t)}\\\
v_{p}(t,S_{\max})=0\,.\end{cases}$ (11)

The analytic solution for (8) is known:

$\displaystyle v_{c}(t,S)$ $\displaystyle=$ $\displaystyle
S\exp(-\delta(T-t))N_{0,1}(d_{1})-K\exp(-r(T-t))N_{0,1}(d_{2}),$
$\displaystyle v_{p}(t,S)$ $\displaystyle=$ $\displaystyle
K\exp(-r(T-t))N_{0,1}(-d_{2})-S\exp(-\delta(T-t))N_{0,1}(-d_{1}),$

with,

$d_{1}=\frac{\log(S/K)+(r-\delta+\sigma^{2}/2)(T-t)}{\sigma\sqrt{T-t}}\,\,,\,\,d_{2}=\frac{\log(S/K)+(r-\delta-\sigma^{2}/2)(T-t)}{\sigma\sqrt{T-t}}$

and $N_{0,1}(x)$ the distribution function of a standard $\mathcal{N}(0,1)$
random variable. Regarding the numerical solution with ANNs, we will use the
methodology introduced in the previous section. In particular, the loss
function is defined as:

$\displaystyle L(v)$
$\displaystyle=\lambda\int_{\Omega}\mid\mathcal{L}(v(t,x))\mid^{p}d\Omega$
$\displaystyle+(1-\lambda)\int_{\partial\Omega}\left(\mid
v(t,x)-G(t,x)\mid^{p}+\mid v(t,x)-H(x)\mid^{p}\right)d\gamma,$ (12)

where functions $G$ and $H$ denote the values of the spatial boundary
conditions and final condition, respectively. The integral terms in the loss
function are approximated by Monte Carlo techniques, as a result, we obtain
the following interior and boundary loss function for the parameter vector
$\theta$:

$\displaystyle\widehat{L}(\theta)=$
$\displaystyle\lambda\frac{1}{n_{I}}\sum_{i=1}^{n_{I}}\mid\mathcal{L}(v(\mathbf{y}^{I}_{i},\theta)))\mid^{p}+$
$\displaystyle(1-\lambda)\left(\frac{1}{n_{B}}\sum_{i=1}^{n_{B}}\mid
v(\mathbf{y}^{B}_{i},\theta)-G(\mathbf{y}^{B}_{i})\mid^{p}+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mid
v(\mathbf{y}^{0}_{i},\theta)-H(\mathbf{x}^{0}_{i})\mid^{p}\right).$ (13)

The collocation points $\\{\mathbf{y}_{i}^{I}\\}_{i=1}^{n_{I}}$ and
$\\{\mathbf{y}_{i}^{B}\\}_{i=1}^{n_{B}}$ are uniformly distributed over the
domain $\Omega$ and the boundary $\partial\widetilde{\Omega}$ and
$\\{\mathbf{y}_{i}^{0}\\}_{i=1}^{n_{0}}$ are uniformly distributed over the
domain $T\times\widetilde{\Omega}$ , respectively and $\mathbf{y}=(t,x)$.

###  3.2 Two underlying assets

We extend the model for one underlying asset to valuing basket options with
two underlying assets. The two-asset prices follow the following dynamics,

$\displaystyle dS_{1_{t}}$
$\displaystyle=(\mu_{1}-\delta_{1})S_{1_{t}}dt+\sigma_{1}S_{1_{t}}dW_{t}^{1},$
$\displaystyle dS_{2_{t}}$
$\displaystyle=(\mu_{2}-\delta_{2})S_{2_{t}}dt+\sigma_{2}S_{2_{t}}dW_{t}^{2},$

where $\mu_{1},\mu_{2}$ are drift terms, $\delta_{1},\delta_{2}$ dividend
yields, the Brownian increments, $dW^{i}$ for $i=1,2$, satisfy
$\mathbb{E}(dW^{i})=0$, and the underlying assets are correlated:

$\text{corr}(W^{1},W^{2})=\rho
t\quad\quad\text{or}\quad\quad\mathbb{E}(dW^{1},dW^{2})=\rho dt\,.$

In the Black-Scholes framework, the two-asset European option price,
$v(t,S_{1},S_{2})$, satisfies the following PDE:

$\displaystyle\begin{cases}\mathcal{L}_{2}(v)=\partial_{t}v+\mathcal{B}v-rv=0\quad(S_{1},S_{2})\in\widetilde{\Omega}\,,\quad
t\in[0,T),\\\ v(T,S_{1},S_{2})=H_{2}(S_{1},S_{2})\,,\end{cases}$ (14)

where the operator $\mathcal{B}$ is defined as follows:

$\displaystyle\mathcal{B}v$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\sigma_{1}^{2}S_{1}^{2}\frac{\partial^{2}v}{\partial
S_{1}^{2}}+\frac{1}{2}\sigma_{2}^{2}S_{2}^{2}\frac{\partial^{2}v}{\partial
S_{2}^{2}}+\rho\sigma_{1}\sigma_{2}S_{1}S_{2}\frac{\partial^{2}v}{\partial
S_{1}\partial S_{2}}$ (15) $\displaystyle+$
$\displaystyle(r-\delta_{1})S_{1}\frac{\partial v}{\partial
S_{1}}+(r-\delta_{2})S_{2}\frac{\partial v}{\partial S_{2}},$

and function $H_{2}(S_{1},S_{2})$ denotes the payoff function. By prescribing
different payoff functions, different options can be defined, like an exchange
option, rainbow option or an average put option. We will deal with the
exchange option, for which an analytic solution is given by the Margrabe’s
formula [15] and the max-on-call rainbow option, for which a closed-form
expression was introduced in [10] and [23]. These particular options are
defined by their payoff functions:

$\displaystyle H_{2}(S_{1},S_{2})$ $\displaystyle=$
$\displaystyle(S_{1}-S_{2})^{+}\,\mbox{ exchange option},$ $\displaystyle
H_{2}(S_{1},S_{2})$ $\displaystyle=$
$\displaystyle(\max(S_{1},S_{2})-K)^{+}\,\mbox{ max-on-call rainbow option}.$

According to the Margrabe’s formula, the fair value of a European exchange
option at time $t$ is given by:

$v(t,S_{1},S_{2})=e^{-\delta_{1}(T-t)}S_{1}(t)N_{0,1}(d_{1})-e^{-\delta_{2}(T-t)}S_{2}(t)N_{01}(d_{2})$
(16)

where $N_{0,1}$ again denotes the cumulative distribution function for the
standard normal,
$\sigma=\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}-2\sigma_{1}\sigma_{2}\rho}$ and

$d_{1}=(\log(S_{1}(t)/S_{2}(t))+(\delta_{2}-\delta_{1}+\sigma^{2}/2)T/\sigma\sqrt{T-t}\quad,\quad
d_{2}=d_{1}-\sigma\sqrt{T-t}\,.$

With the following parameters:

$d_{i}=\frac{\log(S_{i}/k)+(r-\delta_{i}+\frac{\sigma_{i}^{2}}{2})(T-t)}{\sigma_{i}\sqrt{T-t}}\,,$
$\rho_{1}=\frac{\sigma_{1}-\rho\sigma_{2}}{\sigma}\quad\text{and}\quad\rho_{2}=\frac{\sigma_{2}-\rho\sigma_{1}}{\sigma}\,,\quad
i=1,2,$

the closed-form formula for a call on the maximum is given by:

$\displaystyle v_{c}^{\max}(t,S_{1},S_{2})$
$\displaystyle=S_{1}e^{-\delta_{1}(T-t)}M(d_{1},d;\rho_{1})+S_{2}e^{-\delta_{2}(T-t)}M(d_{2},-d+\sigma\sqrt{T-t};\rho_{2})$
$\displaystyle-
Ke^{-r(T-t)}(1-M(-d_{1}+\sigma_{1}\sqrt{T-t},-d_{2}+\sigma_{2}\sqrt{T-t};\rho)),$
(17)

where $M$ is the cumulative bivariate normal distribution

$M(a,b;\rho)=\frac{1}{2\pi\sqrt{1-\rho^{2}}}\int_{-\infty}^{a}\int_{-\infty}^{b}e^{-\frac{x^{2}-2\rho
xy+y^{2}}{2(1-\rho)}}dxdy\,.$

To obtain a numerical solution of the PDE (14), we bound the domain and impose
appropriate boundary conditions. The computational domain should be
sufficiently large, $[0,S_{1\infty}]\times[0,S_{2\infty}]$, where
$S_{1\infty}=S_{2\infty}=4K$ ($K$ the option strike). In the particular case
of the exchange and rainbow max-on-call options, where the analytic solutions
are known, we impose as boundary conditions the analytic option value on each
boundary.

Similar to the one-dimensional problem, we address the European exchange
option problem building the loss function as a sum of the interior and
boundary loss functions, using $\lambda=0.5$.

###  3.3 American options, one underlying asset

As we have introduced in Section 2, we also address the problem for an
American option depending on one underlying asset price. With this aim, we
focus on the linear complementarity formulation.

####  3.3.1 Linear complementarity formulation

We will here consider the linear complementarity problem (LCP) American option
valuation formulation, see, for example, [24, 9], as follows,

$\displaystyle\begin{cases}\mathcal{L}(v)=\partial_{t}v+\mathcal{A}v-rv\leq
0\,,\quad\quad S\in\widetilde{\Omega}\,,t\in[0,T)\,,\\\ v(t,S)\geq H(S),\\\
\mathcal{L}(v)(v-H)=0,\\\ v(T,S)=H(S)\,.\end{cases}$ (18)

This LCP can be rewritten as a nonlinear PDE as follows

$\displaystyle\begin{cases}\max\\{H(S)-v(t,S),\mathcal{L}(v)\\}=0\,,\quad\quad
S\in\widetilde{\Omega}\,,t\in[0,T)\,,\\\ v(T,S)=H(S)\,.\end{cases}$ (19)

Essentially, using the same methodology for solving the European option PDEs,
we address the linear complementarity formulation and its equivalent
formulation as a nonlinear PDE given by (19).

As we introduced in Section 2, the loss function can be formulated using
variance normalization. Moreover, in case of the American option we will also
compute $\lambda$ as the optimal loss weight.

The loss function based on variance normalization depends on the variance of
the network output. For the Black-Scholes American option problem, the loss
function following variance normalization is given by

$\displaystyle L(v)$ $\displaystyle=$
$\displaystyle\frac{\int_{\Omega}\mid\max(H(x)-v(t,x),\mathcal{L}(v(t,x)))\mid^{p}dx}{\int_{\Omega}(\max(\mid
H(x)-v(t,x)\mid,\tilde{\mathcal{L}}(v(t,x))))^{p}dx}$ (20) $\displaystyle+$
$\displaystyle\frac{\int_{\partial\Omega}\left(\mid v(t,x)-G(t,x)\mid^{p}+\mid
v(t,x)-H(x)\mid^{p}\right)d\gamma}{\int_{\partial\Omega}\mid
v(t,x)-\bar{v}\mid^{p}d\gamma},$

where $\tilde{\mathcal{L}}(v(t,x)))$ is defined as follows

$\mathcal{\tilde{L}}(\hat{v})=\mid\partial_{t}\hat{v}\mid+\mid\frac{1}{2}\sigma^{2}S^{2}\partial^{2}_{SS}\hat{v}\mid+\mid(r-\delta)S\partial_{S}\hat{v}\mid+\mid
r\hat{v}\mid\,,$ (21)

function $G$ refers to the boundary conditions imposed in a bounded domain
which are defined as in (11) and function $H$ denotes the final condition.
Moreover, $\bar{v}$ is the mean of $v$ over the corresponding domain, which is
given as

$\bar{v}=\frac{1}{\|\partial\Omega\|}\int_{\partial\Omega}v(t,x)d\Omega\,.$
(22)

Then, approximating each integral term by Monte Carlo techniques the resulting
function is defined as follows

$\displaystyle\widehat{L}(\theta)=$
$\displaystyle\frac{\sum_{i=1}^{n_{I}}\mid\max(H(x^{I}_{i})-v(\mathbf{y}^{I}_{i},\theta),\mathcal{L}(v(\mathbf{y}^{I}_{i},\theta)))\mid^{p}}{\sum_{i=1}^{n_{I}}\max(\mid
H(x^{I}_{i})-v(\mathbf{y}^{I}_{i},\theta)\mid,\mathcal{\tilde{L}}(v(\mathbf{y}^{I}_{i},\theta)))^{p}}+$
$\displaystyle\frac{\frac{1}{n_{B}}\sum_{i=1}^{n_{B}}\mid
v(\mathbf{y}^{B}_{i},\theta)-G(\mathbf{y}^{B}_{i})\mid^{p}+\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\mid
v(\mathbf{y}^{0}_{i},\theta)-H(\mathbf{x}^{0}_{i})\mid^{p}}{\frac{1}{n_{\ast}}\sum_{i=1}^{n_{\ast}}\mid
v(\mathbf{y}^{\ast}_{i},\theta)-\frac{1}{n_{\ast}}\sum_{j=1}^{n_{\ast}}v(\mathbf{y}^{\ast}_{j})\mid^{p}},$
(23)

with $\theta$ containing all parameters of the neural network, vector
$\mathbf{y}=(t,x)$, and the collocation points
$\\{\mathbf{y}_{i}^{\ast}\\}_{i=1}^{n_{\ast}}$ are uniformly distributed over
the boundary $\partial\Omega$.

An alternative is to build the loss function based on an optimal loss weight.
However, optimizing $\lambda$ can be nontrivial.

In order to find the optimal loss weight, we may look for a so-called
$\epsilon$-close solution to the true solution $\hat{v}$, see [16],

$\biggl{|}\frac{\partial^{n}v}{\partial
y_{i}^{n}}-\frac{\partial^{n}\hat{v}}{\partial
y_{i}^{n}}\biggr{|}\leq\epsilon\frac{\partial^{n}\hat{v}}{\partial
y_{i}^{n}},$

for all $n\geq 0$ and $i\in{1,\ldots,d}$, where $d$ is the dimension of the
problem. Satisfying such condition, the value of the optimal loss weight
$\lambda^{*}$ should be:

$\lambda^{*}=\frac{\int_{\partial\Omega}\mid\hat{v}(t,x)\mid^{p}d\gamma}{\int_{\Omega}(\hat{\mathcal{N}}_{I}(t,x,\hat{v}))^{p}d\Omega+\int_{\partial\Omega}\mid\hat{v}(t,x)\mid^{p}d\gamma},$
(24)

where function $\hat{\mathcal{N}}_{I}(x,\hat{v})$ is defined as the function
$\mathcal{N}_{I}(x,\hat{v})$ with each term in absolute value. This expression
of $\lambda^{*}$ is constant when the analytical solution is known. However,
for the American options where the analytical solution is not known, the
optimal loss weight can be computed by approximating the value of $\hat{v}$ in
(24) by the trained solution. Note that in this case, the loss weight is a
function instead of a constant value and is optimized by the neural network.

As a result, the loss function is built in the following way:

$\displaystyle L(v)=$
$\displaystyle\lambda^{*}L_{I}(v)+(1-\lambda^{*})L_{B}(v)$ $\displaystyle=$
$\displaystyle\lambda^{*}\int_{\Omega}\mid\max(H(x)-v(t,x),\mathcal{L}(v(t,x)))\mid^{p}d\Omega+$
$\displaystyle(1-\lambda^{*})\int_{\partial\Omega}\left(\mid
v(t,x)-G(t,x)\mid^{p}+\mid v(t,x)-H(x)\mid^{p}\right)d\gamma,$ (25)

and the optimal loss weight is given in terms of the trained solution $v$, as
follows:

$\lambda^{*}=\frac{\int_{\partial\Omega}\mid
v(t,x)\mid^{p}d\gamma}{\int_{\Omega}(\mathcal{\tilde{L}}_{1}(v))^{p}d\Omega+\int_{\partial\Omega}\mid
v(t,x)\mid^{p}d\gamma},$

with

$\mathcal{\tilde{L}}_{1}(v)=\max(\mid
H(x)-v(t,x)\mid,\mathcal{\tilde{L}}(v)),$

where $\mathcal{\tilde{L}}(v)$ is defined as in (21).

###  3.4 Two-Asset American option

The one underlying asset American option pricing problem is extended to also
price multi-asset American options. We focus on two underlying assets and
formulate the problem as a linear complementarity problem. Based on two asset
prices following correlated geometric Brownian motion, the American option
value can be modeled by the following linear complementarity problem:

$\displaystyle\begin{cases}\mathcal{L}_{2}(v)=\partial_{t}v+\mathcal{B}v-rv\leq
0\,,\quad\quad(S_{1},S_{2})\in\widetilde{\Omega}\,,t\in[0,T),\\\
v(t,S_{1},S_{2})\geq H_{2}(S_{1},S_{2}),\\\ \mathcal{L}_{2}(v)(v-H_{2})=0,\\\
v(T,S_{1},S_{2})=H_{2}(S_{1},S_{2})\,.\end{cases}$ (26)

Operator $\mathcal{B}$ is defined as in (15) and function $H_{2}(S_{1},S_{2})$
denotes the payoff function. In order to compare with the European option
problem, an American call on the maximum is also priced, moreover, we address
a two-asset spread option and a put arithmetic average option. Then, the
payoff functions are defined as:

$\displaystyle H_{2}(S_{1},S_{2})$ $\displaystyle=$
$\displaystyle(\max(S_{1},S_{2})-K)^{+},\;\,\mbox{max-on-call rainbow
option},$ $\displaystyle H_{2}(S_{1},S_{2})$ $\displaystyle=$
$\displaystyle(S_{1}-S_{2}-K)^{+},\;\,\mbox{ asset spread option},$
$\displaystyle H_{2}(S_{1},S_{2})$ $\displaystyle=$
$\displaystyle(K-(S_{1}+S_{2})/2)^{+},\;\,\mbox{arithmetic average put}.$

In order to solve the linear complementarity formulation using numerical
methods, a bounded domain should be considered and appropriate boundary
conditions should be imposed. In particular, we consider a domain large enough
to avoid that the solution is affected by the conditions, in the interested
regions of the asset prices. Whereas for the European option problem the
analytical solution is known and imposed as a boundary condition, for the
American options problem, where the analytical solution is not known, we
should define the appropriate boundary conditions. Then, we start studying at
which boundaries a condition should be imposed. Following [18], that includes
the theory of Fichera [7], we introduce the notation $x_{0}=\tau$,
$x_{1}=S_{1}$, $x_{2}=S_{2}$, and the domain
$\Omega^{\ast}=(0,x_{0}^{\infty})\times(0,x_{1}^{\infty})\times(0,x_{2}^{\infty})$,
where $x_{0}^{\infty}=T$, $x_{1}^{\infty}=S_{1\infty}$ and
$x_{2}^{\infty}=S_{2\infty}$. The boundary of $\Omega^{\ast}$ is,

$\partial\Omega^{\ast}=\bigcup_{i=0}^{2}(\Gamma_{i}^{\ast,-}\cup\Gamma_{i}^{\ast,+}),$

where we use the notation:

$\Gamma_{i}^{\ast,-}=\\{(x_{0},x_{1},x_{2})\in\partial\Omega^{\ast}\,,\/x_{i}=0\\},$
$\Gamma_{i}^{\ast,+}=\\{(x_{0},x_{1},x_{2})\in\partial\Omega^{\ast}\,,\/x_{i}=x_{i}^{\infty}\\}\,.$

Then, the PDE in (26) can be written in the form:

$\sum_{i,j=0}^{2}b_{i,j}\frac{\partial^{2}v}{\partial x_{i}\partial
x_{j}}+\sum_{j=0}^{2}p_{j}\frac{\partial v}{\partial x_{j}}+c_{0}v\leq g_{0},$

where the involved data are defined as follows:

$B(x_{0},x_{1},x_{2})=(b_{ij})=\left(\begin{matrix}0&0&0\\\
0&\frac{1}{2}\sigma_{1}^{2}x_{1}^{2}&\frac{\rho\sigma_{1}\sigma_{2}x_{1}x_{2}}{2}\\\
0&\frac{\rho\sigma_{1}\sigma_{2}x_{1}x_{2}}{2}&\frac{1}{2}\sigma_{2}^{2}x_{2}^{2}\end{matrix}\right)\,,\quad\quad
c_{0}=r\,,$ $p(x_{0},x_{1},x_{2})=(p_{j})=\left(\begin{array}[]{c}-1\\\
(r-\delta_{1})x_{1}\\\ (r-\delta_{2})x_{2}\end{array}\right)\,,\quad\quad
g(x_{0},x_{1},x_{2})=0\,.$

Next, we introduce the following subset of $\Gamma^{\ast}$ in terms of the
normal vector to the boundary pointing inwards $\Omega^{\ast}$,
$\overrightarrow{m}=(m_{0},m_{1},m_{2})$

$\Sigma^{0}=\left\\{x\in\partial\Omega^{\ast}/\sum_{i,j=0}^{2}b_{ij}m_{i}m_{j}=0\right\\}\quad,\quad\Sigma^{1}=\partial\Omega^{0}-\Sigma^{0}\,,$
$\Sigma^{2}=\left\\{x\in\Sigma^{0}/\sum_{i=0}^{2}\left(p_{i}-\sum_{j=0}^{2}\frac{\partial
b_{ij}}{\partial x_{j}}\right)m_{i}\leq 0\right\\}\,.$

In this particular problem, we have:
$\Sigma^{0}=\Gamma_{0}^{\ast,-}\cup\Gamma_{0}^{\ast,+}\cup\Gamma_{1}^{\ast,-}\cup\Gamma_{2}^{\ast,-}$,
$\Sigma^{1}=\Gamma_{1}^{\ast,+}\cup\Gamma_{2}^{\ast,+}$ and
$\Sigma^{2}=\Gamma_{0}^{\ast,+}$. Thus, following [18], the boundary
conditions must be imposed over the subset $\Sigma^{1}\cup\Sigma^{2}$ which
matches with the set
$\Gamma_{0}^{\ast,-}\cup\Gamma_{1}^{\ast,+}\cup\Gamma_{2}^{\ast,+}$. Then, is
not necessary to impose boundary conditions above the boundary where the asset
prices $S_{1}$ and $S_{2}$ are equal zero. Moreover, for simplicity, we assume
that the option value is equal to the payoff when the asset prices $S_{1}$ and
$S_{2}$ take the maximum values.

Next, taking into account the methodologies proposed to solve the one-
dimensional American problem and the two-dimensional European problem, we
propose the loss function to solve the multi-asset American option by
artificial neural network. First of all, we rewrite the linear complementarity
problem (26) as the equivalent nonlinear PDE:

$\displaystyle\begin{cases}\max\\{H_{2}(S_{1},S_{2})-v(t,S_{1},S_{2}),\mathcal{L}_{2}(v)\\}=0,\\\
v(T,S_{1},S_{2})=H_{2}(S_{1},S_{2})\,.\end{cases}$ (27)

Similar to the previous problems, we build the loss function as the sum of the
interior and boundary loss functions as follows:

$\displaystyle L(v)=$ $\displaystyle\lambda L_{I}(v)+(1-\lambda)L_{B}(v)$
$\displaystyle=$
$\displaystyle\lambda\int_{\Omega}\mid\max(H_{2}(\mathbf{x})-v,\mathcal{L}(v))\mid^{p}d\Omega+$
$\displaystyle(1-\lambda)\int_{\partial\Omega}\left(\mid
v(t,\mathbf{x})-G(t,\mathbf{x})\mid^{p}+\mid
v(t,\mathbf{x})-H(\mathbf{x})\mid^{p}\right)d\gamma,$ (28)

where function $G$ refers to the boundary conditions and function $H$ denotes
the final condition imposed for the problems. Note that the loss function is a
generalization of the loss function introduced for the one asset problem and
the integral terms are also approximated by Monte Carlo techniques.

##  4 ANN Option Pricing Results

In this section the European and American options values are computed with the
ANNs based on the loss functions introduced. We apply the unsupervised
learning methodology from the previous section to compute the solutions and
show some results. For the following tests, we have considered the parameter
$p=2$ in the loss functions.

###  4.1 European options

First of all, we discuss the European option single asset results obtained
solving the PDE problem (8) by ANNs.

The results are presented with the loss function introduced in (3.1) and here
we use the basic choice $\lambda=0.5$. Recall that optimal loss weight-based
loss function or the variance normalization technique are especially useful in
the case of nontrivial solutions.

We start with a European put option, with the following parameters values:
$\sigma=0.25$, $r=0.04$, $T=1$, $K=15$, $S_{\infty}=4K$, $\delta=0.0$. In
Figure 1, the ANN-based, trained and the analytical solution are plotted for
two time instances.

Figure 1: European put option for different times instances, $t=0,t=0.5$, with
$\lambda=0.5$.

We measure the accuracy of the solution generated by the ANN by comparing the
relative error of the trained solution $v_{ANN}$ with the analytic solution
$v_{BS}$, as follows:

$error=\frac{\|v_{BS}-v_{ANN}\|_{L^{2}}}{\|v_{BS}\|_{L^{2}}}\,.$ (29)

Figure 2: Error surface for the ANN solution.

In Figure 2 the error throughout the domain is plotted. Clearly, the biggest
error in the ANN solution is found close to the strike price at maturity time
$t=T$, where the payoff is non-smooth. The relative error according to (29)
with $\lambda=0.5$ is equal to $2.23\times 10^{-4}$.

Next, we show some results for a European option depending on two underlying
assets.

The corresponding loss function has been optimized by means of the L-BFGS
algorithm and choosing the $\tanh$ as the activation function. In the last
layer a linear activation function is considered. We have plotted in Figure 3
the ANN solution for the European exchange option. The error, comparing the
approximated ANN solution with the analytical solution given by (16) is also
plotted. Note that the maximum error is obtained for the minimum value of both
asset prices, which is related to $S_{1}/S_{2}$ in the expression of $d_{1}$
in (16).

Figure 3: European exchange option, with parameters:
$\sigma_{1}=\sigma_{2}=0.25$, $\rho=0.1$, $r=0.05$, $\delta_{i}=0.1$,
$S_{1\infty}=S_{2\infty}=60$ and loss weight $\lambda=0.5$.

Due to the relatively big differences in the asset prices $S_{1}$, $S_{2}$ and
the time $t$-values, we have scaled the inputs of the artificial neural
network, i.e. the original domain
$\widetilde{\Omega}=[0,S_{1\infty}]\times[0,S_{2\infty}]$, is scaled to a
dimensionless computational domain, i.e.,
$\widetilde{\Omega}^{\ast}=[0,1]\times[0,1]$. By pricing the option with the
parameters, $S_{1\infty}=S_{2\infty}=4K$, $\sigma_{1}=\sigma_{2}=0.25$,
$\rho=0.1$, $r_{R}=0.04$, $r=0.3$ and $T=1$, for several values of $K$,
modifying the original domain, we found that scaling the input parameters is
not sufficient to obtain highly accurate results for large domain sizes. In
Table 1, the error for a European max-call option is presented, based on
different unscaled domain sizes. It can be observed that as the domain
increases the accuracy of the neural network solution decreases.

($S_{1,\infty},S_{2,\infty}$) | Relative error
---|---
(10 , 10) | $2.58\times 10^{-4}$
(60 , 60) | $3.17\times 10^{-4}$
(120 , 120) | $8.08\times 10^{-4}$
(180 , 180) | $1.71\times 10^{-2}$
(240 , 240) | $2.75\times 10^{-1}$
(300 , 300) | $3.96\times 10^{-1}$
(360 , 360) | $4.30\times 10^{-1}$
Table 1: Relative error for different domains.

In order to understand, the reasons for the degraded accuracy with an
increasing domain size, we have computed the gradients of the interior and
boundary loss functions. In Table 2, we present these values for the European
max-call. The gradient of the interior loss remains constant, note that the
domain is always $[0,1]\times[0,1]$, however, the gradient of the boundary
loss increases with the size domain. Clearly, the interior and boundary loss
functions do not have the same dependency on the domain size.

($S_{1,\infty},S_{2,\infty}$) | $\|\partial L_{I}/\partial\omega\|_{L^{2}}$ | $\|\partial L_{B}/\partial\omega\|_{L^{2}}$
---|---|---
(10 , 10) | $0.4325$ | $8.8515$
(60 , 60) | $0.4325$ | $52.6274$
(120 , 120) | $0.4325$ | $105.1598$
(180 , 180) | $0.4325$ | $157.6923$
(240 , 240) | $0.4325$ | $210.2249$
(300 , 300) | $0.4325$ | $262.7574$
(360 , 360) | $0.4325$ | $315.2899$
Table 2: Gradient values for different domain sizes with standard weights.

We wish to compute accurate approximations of the solution independent the
domain size, and therefore the ANN needs to be modified. The initialization of
the weights is adapted by using a variation of the Xavier initialization. In
particular, the initial values of the weight values in the last layer of the
ANN will be scaled, by multiplying them by the maximum option value. As a
result, we obtain a solution which is accurate independent of the size of the
domain, see Table 3. This adaptation, i.e. the weights having similar
magnitude as the expected largest option value in the output, forms a robust
weight initialization. Moreover, such initialization helps for the interior
and boundary loss functions to have similar sensitivity to the domain size. In
Table 4, we can observe such behaviour, where the rate between both gradients
remains close to $1/3$ when the size of the domain increases. Our results show
that the BFGS optimization doesn’t seem to pick up the gradient if the initial
weights are not sufficiently large. Moreover, similar results can be observed
when the inputs are not scaled.

($S_{1,\infty},S_{2,\infty}$) | Relative error
---|---
(10 , 10) | $3.60\times 10^{-4}$
(60 , 60) | $3.19\times 10^{-4}$
(120 , 120) | $3.56\times 10^{-4}$
(180 , 180) | $4.00\times 10^{-4}$
(240 , 240) | $2.65\times 10^{-4}$
(300 , 300) | $3.14\times 10^{-4}$
(360 , 360) | $4.29\times 10^{-4}$
Table 3: Relative error with scaled weights. ($S_{1,\infty},S_{2,\infty}$) | $\|\partial L_{I}/\partial\omega\|_{L^{2}}$ | $\|\partial L_{B}/\partial\omega\|_{L^{2}}$
---|---|---
(10 , 10) | $26.372$ | $74.790$
(60 , 60) | $949.424$ | $2692.448$
(120 , 120) | $3797.696$ | $10769.792$
(180 , 180) | $8544.816$ | $24232.031$
(240 , 240) | $15190.74$ | $43079.168$
(300 , 300) | $23735.602$ | $67311.2$
(360 , 360) | $34179.266$ | $96928.125$
Table 4: Gradient value for different domains with scaled weights.

Based on the adapted weights initialization, in Table 5, the results for the
European max-call option are presented, and we compare the solution computed
by the ANN with the analytical solution given by (3.2) for some specific asset
prices and different strike values, based on the corresponding loss function
and $\lambda=0.5$. The parameters considered are, $\sigma_{1}=\sigma_{2}=0.2$,
$\rho=0.1$, $\delta_{i}=0.1$, $r=0.05$ and $T=1$. Moreover, the maximum value
of the asset prices is $S_{1\infty}=S_{2\infty}=4K$. Note that the accuracy of
the trained solution is not affected by the size of the domain, in addition,
similar to the one-dimensional case, the maximum error is obtained when the
underlying value is close to the strike price.

strike | ($S_{1}$, $S_{2}$) | ANN | Analytical
---|---|---|---
| ($15,15$) | $-3.16\times 10^{-2}$ | $8.90\times 10^{-5}$
15 | ($10,20$) | $4.58\times 10^{-2}$ | $1.16\times 10^{-2}$
| ($25,5$) | $2.02\times 10^{-1}$ | $2.11\times 10^{-1}$
| ($30,30$) | $-4.23\times 10^{-2}$ | $1.78\times 10^{-4}$
30 | ($20,40$) | $4.75\times 10^{-2}$ | $2.32\times 10^{-2}$
| ($50,10$) | $4.10\times 10^{-1}$ | $4.21\times 10^{-1}$
| ($60,60$) | $1.035\times 10^{-1}$ | $3.56\times 10^{-4}$
60 | ($40,80$) | $1.79\times 10^{-1}$ | $4.63\times 10^{-2}$
| ($100,20$) | $7.79\times 10^{-1}$ | $8.42\times 10^{-1}$
Table 5: European max-call option value.

In Table 6, we present the error for the two-asset European options. The
values are computed based on the expression in (29). We observe very fine
accuracy for the problems.

Option | $\lambda$ | Error
---|---|---
Asset exchange | $0.5$ | $4.16\times 10^{-4}$
Max-call $K=15$ | $0.5$ | $4.55\times 10^{-4}$
Max-call $K=30$ | $0.5$ | $3.51\times 10^{-4}$
Max-call $K=60$ | $0.5$ | $3.83\times 10^{-4}$
Table 6: Error according to the loss weight values.

###  4.2 American options

The goal of this section is to address the American option problem by using
unsupervised learning with the ANN. As for the European option, we compute the
value for one and two underlying assets. However, whereas for the European
option, an analytical solution is known, for the American case, we will use
the option values computed by finite elements (FEM) using the numerical
methods in [2] and [3] to solve the linear complementarity problem for the
American options as the reference.

Similar to the European options problem, the loss function has been optimized
using the L-BFGS algorithm, moreover, the ANN is based on the activation
function $\tanh(x)$. With the aim of comparing both methodologies, we price an
American option with the same parameter data that in the previous example,
considering now, the optimal loss weight, which equals $\lambda\approx 0.90$.
We determine first American options with the following parameter data:
$\sigma=0.25$, $\delta=0.26$, $r=0.3$, $T=1$, $K=15$, $S_{\infty}=4K$. Figure
4 shows the trained solution and the error related to a reference FEM
solution, for all time points. As for the European options, the maximum error
is reached when the asset price is equal to the strike price, close to the
maturity time, where the payoff function is not smooth. In Figure 5, a
comparison of the American option value computed by ANNs or FEM is presented
for different time points. Moreover, the payoff is plotted to demonstrate that
the obstacle condition is satisfied.

Figure 4: American option price with dividends (left). Error surface comparing
with the solution obtained by finite element method (right). Solution obtained
by variance normalization method

Figure 5: American price and the payoff function. Finite element method (green
line) and Neural networks (red line) for different time points. Solution
obtained by variance normalization method

The accuracy of two loss functions for the LCP, one based on optimal loss
weight and another based on variance normalization, is compared by means of
the relative error of the solution, computed in terms of the $L^{2}$-norm,
similar to (29), i.e.

$error=\frac{\|v_{FEM}-v_{ANN}\|_{L^{2}}}{\|v_{FEM}\|_{L^{2}}}\,.$

Very similar accuracy is obtained with both loss functions, $5.38\times
10^{-4}$ (for optimal loss weight) versus $5.82\times 10^{-4}$ (variance
normalization). However, comparing the convergence of both methodologies,
which is presented in Figure 6, we clearly observe that defining the loss
function with a variance normalization (left) the neural network converges
faster than using the optimal loss weight (right). In this figure, the
relative error is plotted for different numbers of iterations of the L-BFGS
algorithm.

Figure 6: Error value (represented in log scale) obtained for different
iterations with the variance normalization method (left) and using the optimal
loss weight (right). The reference solution has been obtained solving the PDE
by the finite element method.

Next, we value the American options depending on two underlying assets.
Optimizing the loss function with the L-BFGS algorithm, with the $\tanh(x)$ as
the activation function, and equal weighting of boundary and interior losses,
$\lambda=0.5$, the following results have been obtained for the three types of
options.

In Table 7, a comparison between the ANN and FEM solutions is shown. We focus
on an American max-call option with several strike values and the following
parameter data, $\rho=0.1$, $\sigma_{1}=\sigma_{2}=0.25$, $r=0.04$,
$\delta=0.01$ and $T=0.5$. Moreover, for the FEM, $75$ time steps have been
considered for the time discretization and the spatial discretization is based
on a $101\times 101$ mesh.

strike | ($S_{1}$, $S_{2}$) | ANN | FEM
---|---|---|---
| ($15,15$) | $2.021$ | $2.066$
15 | ($10,20$) | $5.703$ | $5.643$
| ($25,5$) | $10.996$ | $10.969$
| ($30,30$) | $4.102$ | $4.133$
30 | ($20,40$) | $11.405$ | $11.29$
| ($50,10$) | $21.998$ | $21.938$
| ($60,60$) | $7.916$ | $8.266$
60 | ($40,80$) | $22.753$ | $22.573$
| ($100,20$) | $43.994$ | $43.877$
Table 7: Comparison of American max-call option values.

Figure 7 shows the trained solution for the spread American option with the
following parameter values, $K=15$, $S_{1\infty}=S_{2\infty}=4K$,
$\sigma_{1}=\sigma_{2}=0.25$, $\rho=0.0$, $r_{R}=0.04$, $r=0.3$ and $T=1$, and
the error surface using the FEM solution following [3] as a reference. In
Figure 8, the option value and the difference with the payoff function are
shown.

Figure 7: Two-asset spread American option value in the whole domain (left).
Error surface between the FEM and the ANN solution (right).

Figure 8: Two-Asset spread American option value in a reduced domain (left).
Difference between the ANN solution and the payoff function (right).

Finally, in order to show the accuracy of the method applied to train the ANN
to price American options depending on two asset prices, the relative error is
presented in Table 8.

| Error
---|---
Max-call | $1.73\times 10^{-3}$
Spread | $2.45\times 10^{-3}$
Arithmetic average put | $6.42\times 10^{-3}$
Table 8: Error for different multi-asset American options

Note that the accuracy of the neural network for the American options
depending on two stochastic factors is lower than for the European options.
However, this may be because here a numerical solution is our reference and
not a closed-form expression.

##  5 Conclusions

In this work, classical problems in financial option pricing have been
addressed with artificial neural networks. In particular, following the
classical Black-Scholes model, European and American options depending on one
and two underlying assets have been valued. A new unsupervised learning
methodology is introduced to solve the option value problems based on the PDE
formulation. With this aim, we proposed appropriate loss functions. The
classical Black-Scholes American option pricing problem has been formulated as
a linear complementarity problem.

For the European option problem, the accuracy of the methods was compared to
the analytical solution, whereas, for American options, solutions computed by
the finite element method were used as reference values. For all problems
considered, the final error in the ANN solution was highly satisfactory.
Needless to mention that ANNs can be easily extended to solving higher-
dimensional problems, as they are not drastically affected by the curse of
dimensionality. Finally, the PDE problem formulation can be easily generalized
by introducing counterparty risk which gives rise to nonlinear option
valuation PDEs.

## References

## References

  * Amilon [2003] Amilon, H., 2003.  A neural network versus Black–Scholes:a comparison of pricing and hedging performances.  Journal of Forecasting 22, 317–335.
  * Arregui et al. [2017] Arregui, I., Salvador, B., Vázquez, C., 2017.  PDE models and numerical methods for total value adjustment in European and American options with counterparty risk.  Applied Mathematics and Computation 308, 31–53.
  * Arregui et al. [2020] Arregui, I., Salvador, B., Ševčovič, D., Vázquez, C., 2020\.  Pde models for American options with counterparty risk and two stochastic factors: mathematical analysis and numerical solution.  Computers & Mathematics with Applications DOI: https://doi.org/10.1016/j.camwa.2019.09.014.
  * Becker et al. [2019a] Becker, S., Cheridito, P., Jentzen, A., 2019a.  Pricing and hedging American-style options with deep learning.  ArXiv:1912.11060.
  * Becker et al. [2019b] Becker, S., Cheridito, P., Jentzen, A., Welti, T., 2019b.  Solving high-dimensional optimal stopping problems using deep learning.  ArXiv:1908.01602.
  * Dutta and Shekar [1988] Dutta, S., Shekar, S., 1988\.  Bond rating: a non-conservative application of neural networks.  Proceedings of the IEEE International Conference on Neural Networks 2, 443–450.
  * Fichera [1963] Fichera, G., 1963.  On a unified theory of boundary value problems for elliptic–parabolic equations of second order.  Mathematika 7, 99–122.  Also in Boundary problems in differential equations, University of Wisconsin Press, Madison, 97–120, 1960.
  * Grohs et al. [2019] Grohs, P., Hornung, F., Jentzen, A., Wurstemberger, P.V., 2019\.  A proof that artificial neural networks overcome the curse of dimensionality in the numerical approximation of black-scholes partial differential equations.  Memoirs of the American Mathematical Society ArXiv:1809.02362v1 To appear.
  * Ikonen and Toivanen [2004] Ikonen, S., Toivanen, J., 2004\.  Operator spliting methods for American option pricing.  Applied Matematics Letters 17, 809–814.
  * Johnson [1987] Johnson, H., 1987.  Options on the maximum or the minimum of several assets.  The Journal of Financial and Quantitative Analysis 22, 277–283.
  * Kohler et al. [2010] Kohler, M., Krzyzak, A., Todorovic, N., 2010.  Pricing of highdimensional American options by neural networks.  Mathematical Finance 20, 383–410.
  * Kryzanowski et al. [1993] Kryzanowski, L., Galler, M., Wright, D.W., 1993.  Pricing of high-dimensional American options by neural networks.  Fynancial Analysts Journal 49, 21–27.
  * Liu et al. [2020] Liu, S., Leitao, A., Borovykh, A., Oosterlee, C., 2020\.  On calibration neural networks for extracting implied information from american options ArXiv:2001.11786v1.
  * Liu et al. [2019] Liu, S., Oosterlee, C.W., Bothe, S.M., 2019.  Pricing options and computing implied volatilities using neural networks.
  * Margrabe [1978] Margrabe, W., 1978.  The value of an option to exchange one asset for another.  The Journal of Finance 33, 177–186.
  * van der Meer et al. [2020] van der Meer, R., Oosterlee, C., Borovykh, A., 2020.  Optimally weighted loss functions for solving PDEs with Neural Netwoks .
  * Moody and Utans [1994] Moody, J., Utans, J., 1994.  Architecture selection strategies for neural networks.  Neural Networks in the Capital Markets In Refenes, A. P. (ed.), 141–183.
  * Oleinik and Radkevic [1973] Oleinik, O., Radkevic, E., 1973\.  Second Order Equations with Nonnegative Characteristics Form.  AMS, Plenum Press, New York-London.
  * Raissi et al. [2017] Raissi, M., Karniadakis, G.E., Perdikaris, P., 2017.  Physics informed deep learning (part I): Data-driven solutions of nonlinear partial differential equations.  ArXiv:1711.10561v1.
  * Refenes et al. [1995] Refenes, A.P., Zapranis, A., Francis, G., 1995.  Modelling stock returns in the framework of APT.  Neural Networks in the Capital Markets , 101–125.
  * Singleton and Surkan [1995] Singleton, J.C., Surkan, A.J., 1995\.  Bond rating with neural networks.  Neural Networks in the Capital Markets , 301–307.
  * Sirignano and Spiliopoulos [2018] Sirignano, J., Spiliopoulos, K., 2018\.  DGM:a deep learning algorithm for solving partial differential equations.
  * Stulz [1982] Stulz, R., 1982.  Options on the minimum or the maximum of two risky assets: Analysis and applications.  Journal of Financial Economics 10, 161–185.
  * Wilmott et al. [1996] Wilmott, P., Howison, S., Dewynne, J., 1996.  The Mathematics of Financial Derivatives. A Students Introduction.  Cambridge University Press, Cambridge.
  * Yadav et al. [2015] Yadav, N., Yadav, A., Kumar, M., 2015.  Neural Network Methods for Solving Differential Equations.  Springer, Dordrecht.
  * Zurada [1992] Zurada, J.M., 1992.  Introduction to Artificial Neural Systems.  Wadsworth Publishing Co Inc,.


# EVulHunter: Detecting Fake Transfer Vulnerabilities for EOSIO’s Smart
Contracts at Webassembly-level

Lijin Quan1, Lei Wu2, Haoyu Wang1 1 Beijing University of Posts and
Telecommunications, 2 PeckShield, Inc.

###### Abstract

As one of the representative Delegated Proof-of-Stake (DPoS) blockchain
platforms, EOSIO’s ecosystem grows rapidly in recent years. A number of
vulnerabilities and corresponding attacks of EOSIO’s smart contracts have been
discovered and observed in the wild, which caused a large amount of financial
damages. However, the majority of EOSIO’s smart contracts are not open-
sourced. As a result, the WebAssembly code may become the only available
object to be analyzed in most cases. Unfortunately, current tools are web-
application oriented and cannot be applied to EOSIO WebAssembly code directly,
which makes it more difficult to detect vulnerabilities from those smart
contracts. In this paper, we propose EVulHunter, a static analysis tool that
can be used to detect vulnerabilities from EOSIO WASM code automatically. We
focus on one particular type of vulnerabilities named fake-transfer, and the
exploitation of such vulnerabilities has led to millions of dollars in
damages. To the best of our knowledge, it is the first attempt to build an
automatic tool to detect vulnerabilities of EOSIO’s smart contracts. The
experimental results demonstrate that our tool is able to detect fake transfer
vulnerabilities quickly and precisely. EVulHunter is available on
GitHub111Tool and benchmarks: https://github.com/EVulHunter/EVulHunter and
YouTube222Demo video: https://youtu.be/5SJ0ZJKVZvw.

##  I Introduction

With the growing prosperity of cryptocurrencies like Bitcoin, blockchain
techniques has become more attractive and been adopted by different aeras,
such as financial and logistical systems. Due to the limited throughput (e.g.,
Transaction Per Second, aka TPS) derived from the essence of the Proof-of-Work
(PoW) consensus, traditional blockchain platforms (e.g., Bitcoin and
Ethereum), cannot be used to support high-performing applications. Researchers
proposed different consensuses, such as Proof-of-Share (PoS) and Delegated
Proof-of-Stake (DPoS), to resolve the performance issue.

As one of the most representative DPoS platforms and the first decentralized
operating system, EOSIO [2] has become one of the most active communities all
over the world. EOSIO adopts a multi-threaded mechanism based on its DPoS
consensus protocol [6]. EOSIO claims that it is capable of achieving millions
of TPS. The main native token of EOSIO is called EOS. The amount of EOS tokens
of the holders are used to allocate system resources (e.g., bandwidth and
storage), vote and participate in the on-chain governance, according to the
corresponding proportion of the total stake. Briefly speaking, 21 Block
Producers (BPs) are voted during its launch, and a validate block is generated
every 500 ms block time by those BPs in a round-robin schedule.

Smart contract is a computer protocol which allows users to digitally
negotiate an agreement in a convenient and secure way [3]. In contrast to the
traditional contract law, the transaction costs of the smart contract are
dramatically reduced, and the correctness of its execution is ensured by the
consensus protocol. EOSIO’s smart contracts can be written in C++ , which will
be compiled down to WebAssembly (aka WASM) and executed in EOSIO WASM virtual
machine. WASM is a web standards specifying the binary instruction format for
a stack-based virtual machine, and it can be run in modern web browsers and
other environments [5].

However, it is not easy to guarantee the security of the implementation of
smart contracts, new platforms in particular. A number of vulnerabilities have
been discovered from EOSIO’s smart contracts, including fake EOS transfer,
fake transfer notice and flawed random numbers generators. As such, severe
attacks have been observed in the wild, which caused a large amount of
financial damages.

Unfortunately, most Decentralized Applications (DApps) projects on EOSIO are
not open-sourced, while the web-application oriented analyzed tools cannot be
applied to EOSIO WASM code directly, which makes it more difficult to detect
vulnerabilities from those smart contracts. Specifically, there does exist two
major challenges must be overcome. First of all, since the EOSIO WASM code
does not contain context information (e.g., symbols, function names and type
information), we have to infer the indexes of the functions we are interested
in, including relevant system functions and functions implemented by the
developer. Secondly, the Contract Development Toolkit (CDT) of EOSIO has been
evolving ever since the launch of the mainnet, and we have to cover variants
because of the changes caused by CDT(s).

In this paper, we propose a static analysis tool named EVulHunter which can be
used to detect vulnerabilities from EOSIO WASM code automatically.
Specifically, it first traverses the WASM code and constructs the
corresponding Control Flow Graphs (CFGs), and then detects the existence of
vulnerabilities based on predefined patterns. To facilitate the deep analyses,
EVulHunter also contains a WASM Simulator module which is capable of
performing the simulated code execution. Besides, we focus on one particular
type of vulnerabilities named fake-transfer with two variants, fake EOS
transfer and fake transfer notice. The exploitation of fake-transfer
vulnerability has led to millions of dollars in damages from gambling DApps.

To the best of our knowledge, our work is the first attempt to build an
automatic tool to detect vulnerabilities of EOSIO’s smart contracts. The
experimental results demonstrate that our system is able to detect target
vulnerabilities quickly and precisely with low false positives. Our efforts
can positively shed some light on the exploration of vulnerability detection
of EOSIO’s smart contracts.

##  II Background

###  II-A EOSIO’s Smart Contract

In EOSIO, an action (a base32 encoded 64-bit integer) is used to represent a
single operation, which can be sent individually, or in a combined form, to
serve the purpose of communication between a smart contact and an account [4].
A transaction is composed of one or more actions.

In order to handle requested actions, an apply function is necessary to
dispatch corresponding action handlers for any validate smart contract.
Specifically, apply function listens to all incoming actions and execute the
concrete action handlers accordingly. To this end, the input parameters of
apply function, i.e., the receiver (the account that is currently processing
the action), code (the account that authorized the contract), and action (the
ID of the currently running action), are used as filters to map to the desired
functions that implement particular actions [4].

As one of the 5 fundamental smart contracts of EOSIO333Namely EOSIO.contracts,
i.e., eosio.bios, eosio.token, exchange, eosio.msig, eosio.system, eosio.token
contract is a token standard (like Ethereum’s ERC20 standard) contract which
enables the creation of many different tokens all running on the same contract
but potentially managed by different users [1]. It can be used to create,
issue and transfer tokens [4], and the transfer functionality is related to
the target fake-transfer vulnerabilities.

###  II-B Fake-Transfer Vulnerabilities

Figure 1 shows the lifecycle of a transfer action. Firstly a user sends the
transfer action to the eosio.token contract. After that, the eosio.token
contract will modify the balances accordingly, and then launch the following
two invocations: require_recipient (from) and require_recipient (to), from and
to are source account and destination account respectively. Finally, if any
account has been deployed with a contract, and there does exist a pre-defined
transfer action handler in this contract, the handler would be invoked.

Figure 1: The lifecycle of a transfer action.

Fake EOS Transfer. As mentioned earlier, the code parameter of the apply
function represents the account that authorized the contract. Obviously, the
code parameter of the apply function of the victim (i.e., the recipient to)
should be eosio.token in a normal transfer action. However, if an
implementation of the apply function did not verify the code parameter
properly, it might be deceived into receiving fake EOS tokens and then
executing the further code logic.

Fake EOS Notice. It is an advanced variant of fake-transfer vulnerabilities
with a well-implemented apply function. However, two parameters of the
transfer function, i.e., to and _self, have not been checked for equality.
Specifically, the attacker was able to create an intermediate contract that
forwarding the incoming transfer action to the victim contract. As such, the
victim contract would be misled into believing it was receiving EOS tokens.

##  III Approach

In this section, we will first introduce the overall design of the proposed
EVulHunter system, and then depict each component individually.

Figure 2: Overview of EVulHunter.

###  III-A Overview

Figure 2 illustrates the framework of the proposed system, which accepts the
EOSIO WASM code as the input. EVulHunter is composed of three modules,
including CFG Builder, WASM Simulator and Detector Engine.

We first build the CFGs based on
Octopus444https://github.com/quoscient/octopus, an open source project that is
able to parse EOSIO WASM code. Specifically, the CFGs of the apply function
will be generated in our CFG builder. After that, the entire CFGs will go
through the Detector Engine module. To facilitate the analyses, specific Basic
Blocks (BBs) might be selected to be executed by the WASM Simulator on demand,
where the stack the linear memory would be rebuilt to serve the purpose.

We will describe WASM Simulator and Detector Engine in the following
subsections, III-B and III-C respectively.

###  III-B WASM Simulator

WASM Simulator is designed to be a simple, concise and generic Virtual Machine
(VM) to support further analyses with high extensibility. Basically, this VM
keeps a Stack and Memory structure, which will be modified during tracing
instructions of the WASM code. To facilitate the analysis, several special
patterns have been observed and summarized, including _self, from, to and
strings in a format of 32-bit encoding integer, so that the simulator is able
to recover the semantic type information by mimicking the code execution. For
example, get_local 7 means a popped value from the stack will be assigned to a
variable named local7, even without the knowledge of the value type. After
simulated execution, it is possible to determine that local7 represents a
base32 encoded 64-bit integer and is exactly the “eosio.token” in the context.

Locating Indirect Call Functions. The WASM Simulator can be used to locate
indirect call functions in the apply function, which is the first challenge
mentioned in Section I. Specifically, the link between the apply function and
the target function is implemented with the indirect call mechanism by using
an API named execute_action() in all the versions of CDT. Although different
CDT may have different implementations, the invariant feature is that the
indirect call function pointer is always the last parameter of the
execute_action() function. As the top item in the stack is always the index
number of the indirect call function when calling execute_action, we can use
this feature to determine the indexes of the functions in the corresponding
WASM code.

###  III-C Detector Engine

All detectors are wrapped in the Detector Engine module, which may conduct
deep analyses by interacting with the WASM Simulator when necessary. We have
implemented two detectors for fake EOS transfer and fake transfer notice
vulnerabilities respectively, as follows:

  * •

In case of fake EOS transfer, if none of the functions pre-defined by the
developer is likely to be invoked when the incoming action is not authorized
by eosio.token, which further suggests that the transfer function won’t be
invoked as well, therefore it is not vulnerable to fake EOS transfer
vulnerability; otherwise it is tagged as a vulnerable contract.

  * •

In case of fake transfer notice, if there exists a comparison between to and
_self in the transfer function, it is not vulnerable to fake transfer notice
vulnerability; otherwise it is tagged as a vulnerable contract.

Covering Variants. We have to handle two types of variants. Firstly, variants
derived from the the evolutions of CDT. Secondly, variants come from one
specific verion of CDT, by which even the same logic might be served by
different groups of WASM instructions. We have tried our best to summarize and
cover cases as many as possible, if not all, to make the detectors robust and
complete. All variants being analyzed are summarized in Fig. 3, including
comparison patterns, comparison pairs and corresponding descriptions of
elements being compared shown in elliptical boxes. There are 3 comparison
patterns and 2 comparison pairs respectively, implying 6 combinations in
total.

Figure 3: Summary of Covering Variants.

##  IV Evaluation

TABLE I: The Distribution of our benchmark. | # Vulnerable | # Non-vulnerable | # Total
---|---|---|---
Fake EOS Transfer | 75 | 109 | 184
Fake EOS Notice | 141 | 54 | 195
Total | 159 | 82 | 241
TABLE II: The Overall Experiment Result. | True Positive | False Positive | True Negative | False Negative | Precision | Recall | Accuracy
---|---|---|---|---|---|---|---
Fake EOS Transfer | 75 | 26 | 83 | 0 | 74.26% | 100% | 85.87%
Fake EOS Notice | 141 | 0 | 54 | 0 | 100% | 100% | 100%
Total | 216 | 26 | 137 | 0 | 89.26% | 100% | 93.14%

###  IV-A Benchmark

It appear that no available benchmarks on vulnerable EOS smart contracts in
our community could be used for evaluation. Thus, we propose to manually craft
a benchmark from the known reported vulnerable
Dapps555https://github.com/peckshield/EOS/tree/master/known_dapp_attacks. For
example, it is reported that EOSBet contract was attacked by exploiting the
fake EOS transfer vulnerability in September 2018, and it suffered from the
fake EOS notice vulnerability in October 2018. In this way, we seek to
identify the corresponding vulnerable versions of the reported smart
contracts, and then collect the patched ones with no vulnerabilities.

At last, we have collected 241 EOS smart contracts in total (with 159
vulnerable ones), 184 of them were used as the benchmark of fake EOS transfer
vulnerability, and 195 of them were used to evaluate the fake EOS notice
vulnerability detection. The distribution of our benchmark is shown in Table
I. Note that one smart contract could have both fake EOS transfer and fake EOS
notice vulnerabilities.

###  IV-B Detection Result

Overall Result. Table II shows the overall detection result. For the fake EOS
transfer vulnerability, our tool could achieve an overall accuracy of 86%,
while for the fake EOS notice vulnerability, we could achieve an overall
accuracy of 100%. This result suggests that EVulHunter is able to detect the
fake notice vulnerabilities with high precision and recall. No false negatives
were found in our evaluation, while 26 false positives were found when
detecting fake EOS transfer vulnerability.

False Positives. We have conducted a manual investigation for all the 26 false
positives, which were finally attributed to the history versions of contract
“eosbetdice11”. Specifically, besides the account eosio.token, these contracts
acknowledge the legality of an extra account named eosbettokens 666Probably
related to tokens they issued to support their own business., which breaks our
rules introduced in III-C.

In EOSIO platform, any account can deploy a token contract, however, it is not
trivial for us to collect all possible legal accounts in an anonymous audit.
To alleviate this issue, our system could be easily extended to support a
whitelist mechanism which allows developers to customize extra legal accounts
by themselves.

###  IV-C Performance

We further evaluate the performance of EVulHunter, e.g., the time cost to
detect fake transfer vulnerabilities in a given smart contracts. For the 242
smart contracts in our benchmark, we have recorded the detection time of each
contract, and then analyzed the correlation between the size of smart contract
and the detection time.

As shown in Figure 4, it is obvious that time consumption of vulnerability
detection increased linearly with the increasing of the size of smart
contract, although the increasing rate is very small. In general, the
detection time rages from 1 second to 3 seconds. This result suggests that
EVulHunter is efficient to flag the vulnerabilities in EOS smart contracts,
which could be easily scale to thousands of smart contracts.

Figure 4: The Performance Evaluation of EVulHunter.

##  V Related Work

The security vulnerabilities of smart contracts have attracted great attention
of our research community. However, almost all of the previous studies focus
on analyzing the vulnerabilities in the Ethereum smart contracts [9, 7, 10,
11, 8, 12], as a large number of smart contracts in the Ethereum ecosystem are
open-sourced, and a number of analyzing tools could be leveraged. For example,
Liu et al. [9] presented an analyzer aimed at reentrancy, the type of
vulnerability which cause ‘TheDAO’ hard fork on Ethereum. Atzei et al. [7]
covered more security vulnerabilities and showed a series of attack which
allow attacker to steal money from smart contract. Luu et al. [10] built a
sysbolic execution tool to find potential security bugs. Tikhomirov et al.
[11] provided a static analysis tool to detect problematic language
constructs. Our work is the first attempt to detect fake transfer
vulnerabilities for EOSIO’s smart contracts. We believe our efforts could shed
some light on the detection of vulnerabilities in the EOSIO’s smart contracts.

##  VI Conclusion

In this work, we present the first systematic attempt to automatically detect
fake-transfer vulnerabilities from EOSIO WASM code. Specifically, it first
traverses the WASM code and constructs the corresponding Control Flow Graphs
(CFGs), and then detects the existence of vulnerabilities based on predefined
patterns. The evaluation results demonstrate the system performance, including
speed and accuracy.

There are a number of future lines of work we will explore. For example, our
detection simply relies on constant patterns, which may not be suitable to
model the characteristics of other types of vulnerabilities. As such, some
advanced program analysis techniques (e.g., symbolic execution) could be used
to achieve more accurate results. To this end, we may have to revise the
design of our system to support advanced analyses. Despite the limitations, we
still believe our efforts and observations could positively contribute to the
community.

## References

  * [1] Eos developer documentation, https://www.eosdocs.io/dappdevelopment/tokencontract/
  * [2] Eosio, https://eos.io/
  * [3] Ethereum, https://www.ethereum.org/
  * [4] Official developer portal of eosio, https://developers.eos.io/eosio-cpp/docs/abi
  * [5] Webassembly, https://en.wikipedia.org/wiki/WebAssembly
  * [6] Wiki of eosio, https://en.wikipedia.org/wiki/WebAssembly
  * [7] Nicola Atzei, Massimo Bartoletti, and Tiziana Cimoli.  A survey of attacks on ethereum smart contracts (sok).  In International Conference on Principles of Security and Trust, pages 164–186. Springer, 2017.
  * [8] Bo Jiang, Ye Liu, and WK Chan.  Contractfuzzer: Fuzzing smart contracts for vulnerability detection.  In Proceedings of the 33rd ACM/IEEE International Conference on Automated Software Engineering, pages 259–269. ACM, 2018.
  * [9] Chao Liu, Han Liu, Zhao Cao, Zhong Chen, Bangdao Chen, and Bill Roscoe.  Reguard: finding reentrancy bugs in smart contracts.  In Proceedings of the 40th International Conference on Software Engineering: Companion Proceeedings, pages 65–68. ACM, 2018.
  * [10] Loi Luu, Duc-Hiep Chu, Hrishi Olickel, Prateek Saxena, and Aquinas Hobor.  Making smart contracts smarter.  In Proceedings of CCS 2016, pages 254–269. ACM, 2016.
  * [11] Sergei Tikhomirov, Ekaterina Voskresenskaya, Ivan Ivanitskiy, Ramil Takhaviev, Evgeny Marchenko, and Yaroslav Alexandrov.  Smartcheck: Static analysis of ethereum smart contracts.  In WETSEB 2018, pages 9–16. IEEE, 2018.
  * [12] Petar Tsankov, Andrei Dan, Dana Drachsler-Cohen, Arthur Gervais, Florian Buenzli, and Martin Vechev.  Securify: Practical security analysis of smart contracts.  In Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security, pages 67–82. ACM, 2018.