#  Extreme ATM skew in a local volatility model with discontinuity: joint
density approach

Alexander Gairat111Gearquant. Email address<EMAIL_ADDRESS>and Vadim
Shcherbakov222 Royal Holloway, University of London. Email address:
<EMAIL_ADDRESS>
###### Abstract

This paper concerns a local volatility model in which volatility takes two
possible values, and the specific value depends on whether the underlying
price is above or below a given threshold value. The model is known, and a
number of results have been obtained for it. In particular, a power law
behaviour of the implied volatility skew has been established in the case when
the threshold is taken at the money. This result as well as some others have
been obtained by techniques based on the Laplace transform. The purpose of
this paper is to demonstrate how to obtain similar results by another method.
The proposed alternative approach is based on the natural relationship of the
model with Skew Brownian motion and consists of the systematic use of the
joint distribution of this stochastic process and some of its functionals.

Keywords: Local volatility model, Skew Brownian motion, implied volatility, at
the money skew

##  1 Introduction

This paper concerns a local volatility model (LVM), in which volatility takes
only two possible values. If the underlying price is larger or equal to some
threshold value $R$, then volatility is equal to
$\sigma_{\scriptscriptstyle+}$, and if the underlying price is less than $R$,
then volatility is equal to $\sigma_{\scriptscriptstyle-}$, where
$\sigma_{\scriptscriptstyle+}$ and $\sigma_{\scriptscriptstyle-}$ are given
positive constants. In what follows we refer to this model as the two-valued
LVM. If $\sigma_{\scriptscriptstyle+}=\sigma_{\scriptscriptstyle-}=\sigma$,
then the model is just the classic log-normal model with constant volatility
$\sigma$, under which the celebrated Black-Scholes (BS) formula for the price
of a European option has been obtained. The two-valued LVM is well known, and
a number of results for the model is available (e.g. see [3], [6], [7], [8]
and references therein).

This paper is motivated by the study of the two-valued LVM in [8]. Let us
briefly describe the main results of that paper. First, pricing formulas for
European options under the two-valued LVM were obtained there in the cases
when $S_{0}=R$ and $R=K$, where $S_{0}$ is the spot and $K$ is the strike. The
other results in [8] concern the analysis of the implied volatility surface in
the case $S_{0}=R$, i.e. the threshold is taken at the money (ATM). In
particular, it was shown in that paper that if the strike is
$K=e^{\gamma\sqrt{T}}$, where $T$ is time to maturity, then implied volatility
$\sigma_{{\small{\rm BS}}}(T,K)=\sigma_{{\small{\rm
BS}}}(T,e^{\gamma\sqrt{T}})$ converges, as $T\to 0$, to a smooth function
$\sigma_{{\small{\rm BS}}}(\gamma)$ of $\gamma$, and Taylor’s expansion of the
second order for this limit function was explicitly computed. An exact formula
for the ATM implied volatility skew was obtained, and it was shown that the
skew explodes as $T^{-1/2}$. This is in agreement with the empirical behaviour
of the skew (see [8] for a detailed discussion of this phenomenon).

The results in [8] were obtained by using the Laplace transform and related
research techniques. The purpose of the present paper is to demonstrate how to
obtain some of these results by another method. This alternative method is
based on the natural relationship of the two-valued LVM with Skew Brownian
motion (SBM). The latter is a continuous time Markov process obtained from
standard Brownian motion (BM) by independently choosing with certain fixed
probabilities the signs of the excursions away from the origin. In the case
when these probabilities are equal to $1/2$, then the process is standard BM.
It turns out that if the underlying price follows the two-valued LVM, then the
natural logarithm of the price divided by volatility is a special case of SBM
with a two-valued drift (to be explained). Our approach consists of using the
joint distribution of this process and its functionals, such as the local time
at the origin, the last visit to the origin and the occupation time. The
distribution was obtained in [3], where it was applied to option valuation
under the two-valued LVM. Using this distribution allowed to streamline some
computations in a special case of SBM in [4]. In the present paper, we give
another example of the application of the joint distribution.

First of all, we use the distribution to obtain pricing formulas for European
options in the case $S_{0}=R=1$ (the assumption $R=1$ is just a technical one,
as the general case of $R$ can be readily reduced to this one by dividing the
underlying price by $R$). Note that the option pricing in the general case of
the initial underlying price was considered in [3] (see Section 3 below for
details). However, the case $S_{0}=R=1$ (i.e. the initial underlying price is
at the discontinuity threshold) was not explicitly mentioned in that paper.
Although pricing formulas in this case can be obtained from more general
formulas obtained in [3] by passing to the limit $S_{0}\to 1$, we derive them
here directly by using the aforementioned joint distribution (as we did in [3]
in the general case). This allows us to once again demonstrate the proposed
method. Besides, the corresponding computations are simplified in the case
when the discontinuity threshold is taken at the money. In addition, we show
how to obtain the option prices in this case by combining the joint
distribution with the well-known Dupire’s forward equation.

Our next application of the joint distribution of SBM and its functionals is
the use of the distribution for obtaining the approximation of option prices
in the two-valued LVM by the corresponding BS prices (BS approximation). The
BS approximation is then used for the analysis of implied volatility. Briefly,
the BS-approximation is as follows. Consider a European option with maturity
$T$ and strike $K\geq 1$. Let $\mathsf{C}_{\rm lvm}(K,T)$ be the option price
under the two-valued LVM and let $\mathsf{C}_{\rm
BS}(\sigma_{\scriptscriptstyle+},K,T)$ be the BS price of the same option,
when volatility is equal to $\sigma_{\scriptscriptstyle+}$. Then
$|\mathsf{C}_{\rm lvm}(K,T)-2p\mathsf{C}_{\rm
BS}(\sigma_{\scriptscriptstyle+},K,T)|\leq cT$ for all sufficiently small $T$,
where
$p=\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}$.
A similar approximation holds for prices of European put options. It should be
noted that the BS approximation is not immediately visible from the final
pricing formulas. However, it can be readily obtained, if the option price is
written in terms of an integral of the joint density of the underlying price
process and its functionals. The BS approximation was already noted in [3].
However, it was not properly appreciated in that paper. Despite its
simplicity, the BS approximation turns out to be very precise and useful. In
particular, it allows to obtain some important results concerning the short
term behaviour of implied volatility by rather elementary methods.

The rest of the paper is organised as follows. In Section 2 we formally define
the two-valued LVM and state the key result concerning the aforementioned
joint distribution. Option valuation and the BS approximation for option
prices, is discussed in Section 3. Section 4 concerns the analysis of the
implied volatility surface. Proofs of some results are given in Section 5.

##  2 The model

Start with some notations. Let $(\Omega,{\cal F},\mathsf{P})$ be a probability
space on which all random variables under consideration are defined. The
expectation with respect to the probability measure $\mathsf{P}$ will be
denoted by $\mathsf{E}$. Throughout $(W_{t},\,t\geq 0)$ denotes standard
Brownian motion (BM), and ${\bf 1}_{A}$ denotes the indicator function of a
set or an event $A$.

Before we proceed further, note that without loss of generality we assume
throughout that the threshold value $R$ of the underlying price, where the
volatility changes its value, is $R=1$.

The two-valued LVM that was briefly described in the introduction the
underlying price is given by a stochastic process $S_{t}=(S_{t},\,t\geq 0)$
that follows the equation

$dS_{t}=\sigma(\log(S_{t}))S_{t}dW_{t},$ (1)

with the function $\sigma$ given by

$\sigma(x)=\sigma_{\scriptscriptstyle+}{\bf 1}_{\\{x\geq
0\\}}+\sigma_{\scriptscriptstyle-}{\bf 1}_{\\{x<0\\}},$ (2)

were $\sigma_{\scriptscriptstyle+}>0$ and $\sigma_{\scriptscriptstyle-}>0$ are
given constants.

Further, consider the process $X_{t}=(X_{t},\,t\geq 0)$ defined by

$X_{t}=\begin{cases}\frac{\log(S_{t})}{\sigma_{\scriptscriptstyle+}}&\text{for}\quad
S_{t}\geq 1,\\\
\frac{\log(S_{t})}{\sigma_{\scriptscriptstyle-}}&\text{for}\quad
S_{t}<1.\end{cases}$ (3)

By [3, Lemma 1]), the process $X_{t}$ follows the equation

$dX_{t}=m(X_{t})dt+(p-q)dL_{t}+dW_{t},$ (4)

where

$\displaystyle m(x)$
$\displaystyle=-\frac{\sigma_{\scriptscriptstyle+}}{2}{\bf 1}_{\\{x\geq
0\\}}-\frac{\sigma_{\scriptscriptstyle-}}{2}{\bf 1}_{\\{x<0\\}},$ (5)
$\displaystyle p$
$\displaystyle=\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}},$
(6) $\displaystyle q$
$\displaystyle=1-p=\frac{\sigma_{\scriptscriptstyle+}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}},$
(7)

and

$L_{T}=L^{(X)}_{T}=\lim\limits_{\varepsilon\to
0}\frac{1}{2\varepsilon}\int\limits_{0}^{T}1_{\\{-\varepsilon\leq
X_{u}\leq\varepsilon\\}}du$ (8)

is the local time of the process at the origin.

The process $X_{t}$ is a special case of Skew Brownian motion (SBM) with a
two-valued drift. Recall that SBM (without drift) is obtained from standard BM
by independently choosing with certain probabilities the signs of the
excursions away from the origin. An excursion is chosen to be positive with a
fixed probability $p$ and negative with probability $q=1-p$ (if these
probabilities are equal to $\frac{1}{2}$, then the process is standard BM).
The process $X_{t}$ is SBM with probabilities $p$ and $q$ given by (6) and (7)
respectively, and the two-valued drift (5) (see Section B).

As we already mentioned in the introduction, a key ingredient in our analysis
of the two-valued LVM is the use of the joint distribution of the process
$X_{t}$ and some of its functionals. These functionals include the local time
of the process and the following quantities. Namely, given $T>0$ let

$\displaystyle\tau_{0}$
$\displaystyle=\min\left\\{t\in[0,T]:X_{t}=0\right\\},$ (9)
$\displaystyle\tau$ $\displaystyle=\max\\{t\in(0,T]:X_{t}=0\\},$ (10)

be the first and the last visits to the origin respectively (on the interval
$[0,T]$), and let

$\displaystyle V$ $\displaystyle=\int\limits_{\tau_{0}}^{\tau}1_{\\{X_{t}\geq
0\\}}dt$ (11)

be the occupation time of the non-negative half line on the interval
$[\tau_{0},\tau]$.

The distribution of interest is given in Theorem 1 below. Note that the
theorem is a special case of a more general result for SBM obtained in [3]
(see Theorem 7 in Section 7).

######  Theorem 1.

Let $(X_{t},\,t\geq 0)$ be the process defined in (3). Let $L_{T}$, $\tau$ and
$V$ be quantities of this process defined by equations (8), (10) and (11)
respectively, and let $X_{0}=0$. Then, the joint density of random variables
$\tau,V,X_{T}$ and $L_{T}$ is given by

$f_{T}(t,v,x,l)=2a(x)h(v,lp)h(t-v,lq)h(T-t,x)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}v+\sigma_{\scriptscriptstyle-}^{2}(t-v)+\sigma^{2}(x)(T-t)}{8}-\frac{\sigma(x)}{2}x},$
(12)

for $0\leq v\leq t\leq T$ and $l\geq 0$, where

$a(x)=p{\bf 1}_{\\{x\geq 0\\}}+q{\bf 1}_{\\{x<0\\}},$ (13)

probabilities $p$ and $q$ are defined in (6) and (7) respectively, and

$h(s,y)=\frac{|y|}{\sqrt{2\pi
s^{3}}}e^{-\frac{y^{2}}{2s}},\,\,y\in{\mathbb{R}},s\in{\mathbb{R}}_{+},$ (14)

is the density of the first passage time to zero of the standard BM starting
at $y$.

######  Remark 1.

Note that it is convenient to rewrite the density in terms of the following
variables $u=t-v$ and $s=T-t$ for $t\in[0,T]$. It is easy to see that if
$X_{0}=0$ and $x\geq 0$, then the variable $u$ is the occupation time of the
negative half-line, and the variable $v+s$ is the total occupation time of the
positive half-line, and if $X_{0}=0$ and $x<0$, then the occupation time of
the negative half-line is $u+s$, and the total occupation time of the positive
half-line is $v$. In these terms we have that

$\begin{split}f_{T}(t,v,x,l)&=\widetilde{f}_{T}(v,u,s,x,l)=2a(x)h(v,lp)h(u,lq)h(s,x)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}v+\sigma_{\scriptscriptstyle-}^{2}u+\sigma^{2}(x)s}{2}-\frac{\sigma(x)}{2}x},\end{split}$
(15)

for $(v,u,s):v,u,s\geq 0,\,v+u+s=T$ and $l\geq 0$.

######  Theorem 2.

Let $p(0,x,T)$ be the probability density function of $X_{T}$ given that
$X_{0}=0$. Then

$\displaystyle p(0,x,T)$
$\displaystyle=2a(x)e^{-\frac{\sigma(x)}{2}x}\int_{0}^{T}\phi(T-s)h(s,x)e^{-\frac{\sigma^{2}(x)}{8}s}ds$
$\displaystyle=\begin{cases}2pe^{-\frac{\sigma_{\scriptscriptstyle+}}{2}x}\int_{0}^{T}\phi(T-s)h(s,x)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}ds&\quad\text{for}\quad
x\geq 0,\\\
2qe^{-\frac{\sigma_{\scriptscriptstyle-}}{2}x}\int_{0}^{T}\phi(T-s)h(s,x)e^{-\frac{\sigma_{\scriptscriptstyle-}^{2}}{8}s}ds&\quad\text{for}\quad
x<0,\\\ \end{cases}$

where

$\begin{split}\phi(t)&=\frac{1}{\sqrt{2\pi
t}(\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-})}\left(\sigma_{\scriptscriptstyle+}e^{-\frac{\sigma_{\scriptscriptstyle-}^{2}}{8}t}-\sigma_{\scriptscriptstyle-}e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}t}\right)\\\
&+\frac{1}{2}\frac{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}{(\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-})}\left(\mathcal{N}\left(\frac{\sqrt{t}\sigma_{\scriptscriptstyle-}}{2}\right)-\mathcal{N}\left(\frac{\sqrt{t}\sigma_{\scriptscriptstyle+}}{2}\right)\right)\end{split}$
(16)

and

$\mathcal{N}(z)=\int_{-\infty}^{z}\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}dx\quad\text{for}\quad
z\in{\mathbb{R}}.$ (17)

Theorem 2 is proved in Section 5.

######  Remark 2.

The function $\phi$ defined in (16) has the following probabilistic sense.
Recall that the function $h(s,x)$ is the density of the first passage time to
$0$ of the standard BM starting at $x$. The distribution of the first passage
time converges, as $x\to 0$, to the distribution concentrated at $0$. It
follows from this fact and Theorem 2 that

$\phi(T)=\frac{1}{2p}\lim\limits_{x\downarrow
0}p(0,x,T)=\frac{1}{2q}\lim\limits_{x\uparrow 0}p(0,x,T).$

##  3 Option valuation

Let us briefly recall results of [3] concerning the option valuation under the
two-valued LVM. Pricing formulas were obtained in that paper for knock-in call
options in the following cases: $R=1,S_{0}\geq 1,K>1$ and $R=1,S_{0}<1,K>1$.
By combining these results with the Black-Scholes prices for knock-out call
options one can obtain prices of European call options for other values of
$S_{0}$ and $K$. The pricing formula in the case $R=1,S_{0}\geq 1,K>1$ in [3]
is given in terms of a single integral, where an integrand is analytically
expressed in terms of the cumulative distribution function (cdf) of the
standard normal distribution $N(0,1$). In the special case $R=S_{0}=1,K>1$ one
gets the price of the European option call option. This case is considered in
Theorem 3 below. The pricing formula in the case $R=1,S_{0}<1,K>1$ in [3] is
also given in terms of a single integral with an integrand analytically
expressed in terms of the standard normal distribution and a bivariate normal
distribution.

###  3.1 Pricing formulas

Let

$\psi(a,s,k):=\begin{cases}\int_{k}^{\infty}e^{ax}h(s,x)dx=\frac{1}{\sqrt{2\pi
s}}e^{ak-\frac{k^{2}}{2s}}+ae^{\frac{a^{2}}{2}s}\mathcal{N}\left(\frac{as-k}{\sqrt{s}}\right)&\text{
for }k\geq 0,\\\ \int_{-\infty}^{k}e^{ax}h(s,x)dx=\frac{1}{\sqrt{2\pi
s}}e^{ak-\frac{k^{2}}{2s}}-ae^{\frac{a^{2}}{2}s}\mathcal{N}\left(\frac{k-as}{\sqrt{s}}\right)&\text{
for }k<0.\end{cases}$ (18)

Note that the equation for the function $\psi$ can be rewritten as follows

$\psi(a,s,k)=\frac{1}{\sqrt{2\pi s}}e^{ak-\frac{k^{2}}{2s}}+\rm sgn(k)\cdot
a\cdot e^{\frac{a^{2}}{2}s}\mathcal{N}\left(\rm
sgn(k)\frac{as-k}{\sqrt{s}}\right)\quad\text{for}\quad k\in{\mathbb{R}},$ (19)

where

$\rm sgn(k)=\begin{cases}1&\text{for }k\geq 0,\\\ -1&\text{for
}k<0.\end{cases}$

Define

$F\left(T,a,k\right)=\int_{0}^{T}\phi(T-s)\psi(a,s,k)e^{-\frac{a^{2}}{2}s}ds,$
(20)

where $\phi$ is the function defined in (16).

Consider European options with strike $K$ and time to expiry $T$. Let
$\mathsf{C}_{\rm lvm}(S_{0},K,T)$ and $\mathsf{P}_{\rm lvm}(S_{0},K,T)$ be the
price of the call option and the put option respectively, when the initial
value of the underlying price is $S_{0}$.

######  Theorem 3.

If $S_{0}=1$ and $K\geq 1$, then

$\mathsf{C}_{\rm
lvm}(1,K,T)=2p\left(F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},k\right)-e^{\sigma_{\scriptscriptstyle+}k}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},k\right)\right),$
(21)

where $k=\log(K)/\sigma_{\scriptscriptstyle+}$.

If $S_{0}=1$ and $K<1$, then

$\mathsf{P}_{\rm
lvm}(1,K,T)=2q\left(e^{\sigma_{\scriptscriptstyle-}k}F\left(T,-\frac{\sigma_{\scriptscriptstyle-}}{2},k\right)-F\left(T,\frac{\sigma_{\scriptscriptstyle-}}{2},k\right)\right),$
(22)

where $k=\log(K)/\sigma_{\scriptscriptstyle-}$.

By Theorem 3 we immediately get the following equation for the ATM

$\begin{split}\mathsf{V}_{\rm atm}(T)&=\mathsf{C}_{\rm
lvm}(1,1,T)=2p\left(F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)-F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)\right)\\\
&=\mathsf{P}_{\rm
lvm}(1,1,T)=2q\left(F\left(T,-\frac{\sigma_{\scriptscriptstyle-}}{2},0\right)-F\left(T,\frac{\sigma_{\scriptscriptstyle-}}{2},0\right)\right),\end{split}$
(23)

which still involves the function $F$. Corollary 1 below shows that the above
equation for the ATM price can be simplified in such a way that the ATM price
is analytically expressed in terms of the error function

${\rm Erf}(x)=\frac{2}{\sqrt{\pi}}\int_{0}^{x}e^{-t^{2}}dt.$

######  Corollary 1 (ATM price).

If $S_{0}=K=1$, then

$\begin{split}\mathsf{V}_{\rm
atm}(T)&=\frac{\sigma_{\scriptscriptstyle-}^{2}\sigma_{\scriptscriptstyle+}^{2}}{4(\sigma_{\scriptscriptstyle-}^{2}-\sigma_{\scriptscriptstyle+}^{2})}\left[\frac{\sqrt{8T}}{\sigma_{\scriptscriptstyle+}\sqrt{\pi}}e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}T}-\frac{\sqrt{8T}}{\sigma_{\scriptscriptstyle-}\sqrt{\pi}}e^{-\frac{\sigma_{\scriptscriptstyle-}^{2}}{8}T}\right.\\\
&\left.+\left(\frac{4}{\sigma_{\scriptscriptstyle+}^{2}}+T\right){\rm
Erf}\left(\frac{\sigma_{\scriptscriptstyle+}\sqrt{T}}{\sqrt{8}}\right)-\left(\frac{4}{\sigma_{\scriptscriptstyle-}^{2}}+T\right){\rm
Erf}\left(\frac{\sigma_{\scriptscriptstyle-}\sqrt{T}}{\sqrt{8}}\right)\right]\end{split}$
(24)

Proofs of both Theorem 3 and Corollary 1 are given in Section 5.

######  Remark 3.

It should be noted that the option pricing formulas in Theorem 3 differ from
those that are given in [8]. In both papers the option price is given by a
single integral, but the corresponding integrands differ. It might be of
interest to investigate the relationship between the two variants of. At the
same time, the ATM price in Corollary 1 is exactly the same as the one in [8].
In addition, note that the ATM price can be rewritten as follows

$\mathsf{V}_{\rm
atm}(T)=\frac{\sigma_{\scriptscriptstyle-}^{2}\sigma_{\scriptscriptstyle+}^{2}}{4(\sigma_{\scriptscriptstyle-}^{2}-\sigma_{\scriptscriptstyle+}^{2})}(I(\sigma_{\scriptscriptstyle+},T)-I(\sigma_{\scriptscriptstyle-},T)),$

where

$I(z,T)=\frac{\sqrt{8T}}{z\sqrt{\pi}}e^{-\frac{z^{2}}{8}T}+\left(\frac{4}{z^{2}}+T\right){\rm
Erf}\left(\frac{z\sqrt{T}}{\sqrt{8}}\right).$

###  3.2 Option prices and Dupire’s forward equation

In this section we provide (in the Theorem 4 below) another representation of
option prices. This representation is based on the Dupire’s forward equation
([1], [2]), which we are going to recall.

If the underlying price follows the local volatility model

$dS_{t}=\sigma(t,S_{t})S_{t}dW_{t},$

then the price $\mathsf{C}_{\rm lvm}(K,T)$ of a European call option with
strike $K$ and time to expiry $T$ satisfies the equation (the forward
equation)

$\frac{\partial\mathsf{C}_{\rm lvm}(K,T)}{\partial
T}=\frac{1}{2}K^{2}\sigma^{2}(T,K)\frac{\partial^{2}\mathsf{C}_{\rm
lvm}(K,T)}{\partial K^{2}}.$ (25)

It follows from the forward equation that

$\displaystyle\mathsf{C}_{\rm lvm}\left(K,T\right)-\left(S_{0}-K\right)_{+}$
$\displaystyle=\int_{0}^{T}\frac{\partial}{\partial t}\mathsf{C}_{\rm
lvm}\left(K,t\right)dt$
$\displaystyle=\frac{K^{2}}{2}\int_{0}^{T}\sigma^{2}(t,K)\frac{\partial^{2}}{\partial
K^{2}}\mathsf{C}_{\rm lvm}\left(K,t\right)dt$
$\displaystyle=\frac{K^{2}}{2}\int_{0}^{T}\sigma^{2}(t,K)\frac{\partial^{2}}{\partial
K^{2}}\mathsf{E}\left(\max\left(S_{t}-K,0\right)|S_{0}\right)dt$
$\displaystyle=\frac{K^{2}}{2}\int_{0}^{T}\sigma^{2}(t,K)\mathsf{E}\left(\delta\left(S_{t}-K\right)|S_{0}\right)dt,$

where $\delta(\cdot)$ is the delta-function. Noting that
$\mathsf{E}\left(\delta\left(S_{t}-K\right)|S_{0}\right)=p_{S}\left(S_{0},K,t\right)$,
where $p_{S}\left(S_{0},\cdot,t\right)$, is the probability density function
of $S_{t}$ given $S_{0}$, we arrive to the equation

$\mathsf{C}_{\rm
lvm}\left(K,T\right)-\left(S_{0}-K\right)_{+}=\frac{K^{2}}{2}\int_{0}^{T}\sigma^{2}(t,K)p_{S}\left(S_{0},K,t\right)dt,$
(26)

which we are going to use in the proof of Theorem 4 below.

######  Theorem 4.

$\mathsf{C}_{\rm lvm}(1,K,T)=\sqrt{K}\int_{0}^{T}\mathsf{V}_{\rm
atm}(T-s)h(s,\log(K)/\sigma_{\scriptscriptstyle+})e^{-\frac{1}{8}\sigma_{\scriptscriptstyle+}^{2}s}ds\quad\text{for}\quad
K\geq 1,$ (27)

and

$\mathsf{P}_{\rm lvm}(1,K,T)=\frac{1}{\sqrt{K}}\int_{0}^{T}\mathsf{V}_{\rm
atm}(T-s)h(s,-\log(K)/\sigma_{\scriptscriptstyle-})e^{-\frac{1}{8}\sigma_{\scriptscriptstyle-}^{2}s}ds\quad\text{for}\quad
K<1,$ (28)

where in both cases $\mathsf{V}_{\rm atm}$ is the ATM price.

The proof of Theorem 4 is given in Section 5.

###  3.3 Black-Scholes approximation

In this section we obtain the approximation of option prices $\mathsf{C}_{\rm
lvm}(1,K,T)$ and $\mathsf{P}_{\rm lvm}(1,K,T)$ by the corresponding BS prices.

Denote by $\mathsf{C}_{\rm BS}(\sigma,S_{0},K,T)$ and $\mathsf{P}_{\rm
BS}(\sigma,S_{0},K,T)$ the BS prices of European call option and European put
option respectively with strike $K$ and time to maturity $T$, when volatility
of the underlying asset is equal to $\sigma$.

######  Theorem 5.

If $K\geq 1$, then

$|\mathsf{C}_{\rm lvm}(1,K,T)-2p\mathsf{C}_{\rm
BS}(\sigma_{\scriptscriptstyle+},1,K,T)|\leq cT,$

and if $K<1$, then

$|\mathsf{P}_{\rm lvm}(1,K,T)-2q\mathsf{P}_{\rm
BS}(\sigma_{\scriptscriptstyle-},1,K,T)|\leq cT$

for some constant $c$ and all sufficiently small $T$.

We already mentioned in the introduction, that despite its simplicity, the BS
approximation is surprisingly accurate. In particular, this approximation
allows to obtain some key results concerning the behaviour of implied
volatility in the case of short term maturities.

##  4 Implied volatility

In this section we use the BS approximation to obtain some results concerning
implied volatility that were obtained in [8]. Implied volatility
$\sigma_{{\small{\rm BS}}}(T,k)$ is considered as a function of maturity $T$
and log-moneyness $k=\log(K/S_{0})$. Note that in the case $S_{0}=1$ we have
that $k=\log(K)$.

Start with a remark that is almost verbatim to Remark 3.4 in [8] concerning
the ATM implied volatility $\sigma_{{\small{\rm BS}}}(T,0)$. By definition,
$\sigma_{{\small{\rm BS}}}(T,0)$ is the solution of the equation

$\mathsf{C}_{\rm BS}(\sigma_{{\small{\rm BS}}}(T,0),1,1,T)=\mathsf{V}_{\rm
atm}(T).$

Given volatility $\sigma$we have that

$\mathsf{C}_{\rm
BS}(\sigma,1,1,T)=\mathcal{N}\left(\frac{\sigma\sqrt{T}}{2}\right)-\mathcal{N}\left(-\frac{\sigma\sqrt{T}}{2}\right)={\rm
Erf}\left(\frac{\sigma\sqrt{T}}{\sqrt{8}}\right).$

Therefore,

$\sigma_{{\small{\rm BS}}}(T,0)=\frac{\sqrt{8}}{\sqrt{T}}{\rm
Erf}^{-1}(\mathsf{V}_{\rm atm}(T)).$ (29)

Further, using that

${\rm
Erf}^{-1}(z)=\frac{\sqrt{\pi}}{2}\left(z+\frac{\pi}{12}z^{3}\right)+o(z^{3}),\quad\text{as}\quad
z\to 0,$

gives the short term expansion for the ATM implied volatility

$\sigma_{{\small{\rm
BS}}}\left(T,0\right)=2\frac{\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\left(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}\right)}-\frac{\left(\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}\right)^{2}\left(\sigma_{\scriptscriptstyle-}-\sigma_{\scriptscriptstyle+}\right)^{2}}{12\left(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}\right)^{3}}T+o(T),\quad\text{as}\quad
T\to 0,$ (30)

and, hence,

$\sigma_{{\small{\rm
BS}}}\left(0,0\right)=2\frac{\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\left(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}\right)}.$

###  4.1 Implied volatility in the central limit regime

In this section we consider the short term behaviour of implied volatility in
the case when $K=e^{\gamma\sqrt{T}}$ (i.e. in the central limit regime). This
case was considered in Theorem 3.1 in [8]. One of the results of the theorem
is the following equation

$\begin{split}\lim_{T\to 0}\sigma_{{\small{\rm
BS}}}(T,\gamma\sqrt{T})&=\frac{2\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}+\frac{\sqrt{\pi}}{\sqrt{2}}\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\gamma+\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}\left(\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{2\left(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}\right)}-\rm
sgn(\gamma)\right)\gamma^{2}\\\ &+o(\gamma^{2}),\quad\text{as}\quad\gamma\to
0.\end{split}$ (31)

In other words, the implied volatility $\sigma_{{\small{\rm
BS}}}(T,\gamma\sqrt{T})$ can be approximated for short term maturities $T$ by
a quadratic function of $\gamma$, which can be computed explicitly. Below we
compute this quadratic function by using the BS approximation derived in
Theorem 5.

For definiteness assume that $\gamma>0$ (i.e. $K=e^{\gamma\sqrt{T}}\geq 1$)
and use our result for prices of call options. Let $\mathsf{C}_{\rm
lvm}(1,e^{\gamma\sqrt{T}},T)$ be the call price in the two valued LVM (Theorem
3). The equation for the implied, volatility $\sigma_{{\small{\rm BS}}}$ is

$\mathsf{C}_{\rm BS}(\sigma_{{\small{\rm
BS}}},1,e^{\gamma\sqrt{T}},T)=\mathsf{C}_{\rm lvm}(1,e^{\gamma\sqrt{T}},T),$
(32)

where the left hand side is the BS price of the call option with maturity $T$
and the strike $K=e^{\gamma\sqrt{T}}$. Then

$\mathsf{C}_{\rm
BS}(\sigma,1,e^{\gamma\sqrt{T}},T)=\mathcal{N}(d_{1})-e^{\gamma\sqrt{T}}\mathcal{N}(d_{0}),$

where $d_{1}=-\frac{\gamma}{\sigma}+\frac{\sigma\sqrt{T}}{2}$ and
$d_{0}=-\frac{\gamma}{\sigma}-\frac{\sigma\sqrt{T}}{2}$. By Taylor’s theorem
we have that

$\displaystyle\mathsf{C}_{\rm BS}(\sigma,1,e^{\gamma\sqrt{T}},T)$
$\displaystyle=\left(\frac{\sigma}{\sqrt{2\pi}}-\frac{\gamma}{2}+\frac{\gamma^{2}}{\sigma\sqrt{2\pi}}\right)\sqrt{T}+o(\sqrt{T}),\quad\text{as}\quad
T\to 0.$ (33)

Further, combining the BS approximation (Theorem 5) for the right hand side of
(32) with (33) gives that

$\mathsf{C}_{\rm lvm}(1,e^{\gamma\sqrt{T}},T)=2p\mathsf{C}_{\rm
BS}(\sigma_{\scriptscriptstyle+},1,e^{\gamma\sqrt{T}},T)=2p\left(\frac{\sigma_{\scriptscriptstyle+}}{\sqrt{2\pi}}-\frac{\gamma}{2}+\frac{\gamma^{2}}{\sigma_{\scriptscriptstyle+}\sqrt{2\pi}}\right)\sqrt{T}+o(\sqrt{T}).$
(34)

Replacing both the left hand and the right hand sides of equation (32) by
their approximations (provided by equations (33) and (34) respectively) we
obtain the following equation

$\left(\frac{\sigma}{\sqrt{2\pi}}-\frac{\gamma}{2}+\frac{\gamma^{2}}{\sigma\sqrt{2\pi}}\right)\sqrt{T}=2p\left(\frac{\sigma_{\scriptscriptstyle+}}{\sqrt{2\pi}}-\frac{\gamma}{2}+\frac{\gamma^{2}}{\sigma_{\scriptscriptstyle+}\sqrt{2\pi}}\right)\sqrt{T}+o(\sqrt{T}),$
(35)

which means that under the assumptions made the implied volatility
$\sigma_{{\small{\rm BS}}}(\gamma,e^{\gamma\sqrt{T}})$ converges to a limit,
as $T\to 0$, and, moreover, this limit can be estimated by the solution of the
equation

$\frac{\sigma}{\sqrt{2\pi}}-\frac{\gamma}{2}+\frac{\gamma^{2}}{\sigma\sqrt{2\pi}}=p\sigma_{\scriptscriptstyle+}\sqrt{\frac{2}{\pi}}-p\gamma+\frac{p}{\sigma_{\scriptscriptstyle+}}\sqrt{\frac{2}{\pi}}\gamma^{2}.$
(36)

It is easy to see that equation (36) is basically a quadratic equation for the
unknown $\sigma$ with coefficients analytically depending on $\gamma$. This
implies that the solution is an analytic function $\sigma(\gamma)$ of
$\gamma$. Consider Taylor’s expansion of the second order for this function at
$\gamma=0$, that is

$\sigma(\gamma)=\left(c_{0}+c_{1}\gamma+\frac{1}{2}c_{2}\gamma^{2}\right)+o\left(\gamma^{2}\right),\quad\text{as}\quad\gamma\to
0,$ (37)

where $c_{0}$, $c_{1}$ and $c_{2}$ denote the values at $0$ of the function
itself, its 1st and 2nd derivatives respectively. Using this expansion for
approximating the left hand side of (36) gives the following equation

$\begin{split}\frac{c_{0}+c_{1}\gamma+\frac{1}{2}c_{2}\gamma^{2}}{\sqrt{2\pi}}&-\frac{\gamma}{2}+\frac{\gamma^{2}}{\sqrt{2\pi}}\frac{1}{(c_{0}+c_{1}\gamma+\frac{1}{2}c_{2}\gamma^{2})}\\\
&=\frac{c_{0}}{\sqrt{2\pi}}+\left(\frac{c_{1}}{\sqrt{2\pi}}-\frac{1}{2}\right)\gamma+\frac{1}{\sqrt{2\pi}}\left(\frac{c_{2}}{2}+\frac{1}{c_{0}}\right)\gamma^{2}+o(\gamma^{2}),\end{split}$
(38)

which, in turn, implies that

$\frac{c_{0}}{\sqrt{2\pi}}+\left(\frac{c_{1}}{\sqrt{2\pi}}-\frac{1}{2}\right)\gamma+\frac{1}{\sqrt{2\pi}}\left(\frac{\sigma_{2}}{2}+\frac{1}{c_{0}}\right)\gamma^{2}=p\sigma_{\scriptscriptstyle+}\sqrt{\frac{2}{\pi}}-p\gamma+\frac{p}{\sigma_{\scriptscriptstyle+}}\sqrt{\frac{2}{\pi}}\gamma^{2}+o(\gamma^{2}).$
(39)

Equating coefficients at $\gamma^{i},\,i=0,1,2$ in (39) we obtain that

$c_{0}=2p\sigma_{\scriptscriptstyle+},\quad
c_{1}=\sqrt{\frac{2}{\pi}}(1-2p)\quad\text{and}\quad
c_{2}=\frac{4p^{2}-1}{p\sigma_{\scriptscriptstyle+}},$

and, hence,

$\begin{split}\sigma(\gamma)&=2p\sigma_{\scriptscriptstyle+}+\frac{\sqrt{2}}{\sqrt{\pi}}(1-2p)\gamma+\frac{1}{2}\frac{4p^{2}-1}{p\sigma_{\scriptscriptstyle+}}\gamma^{2}+o(\gamma^{2})\\\
&=\frac{2\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}+\frac{\sqrt{\pi}}{\sqrt{2}}\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\gamma+\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}\left(\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{2\left(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}\right)}-1\right)\gamma^{2}+o(\gamma^{2}).\end{split}$
(40)

Alternatively, one can use the put price and repeat the above argument in the
case when $\gamma<0$, i.e. $K=e^{\gamma\sqrt{T}}<1$, to obtain that

$\begin{split}\sigma(\gamma)&=2q\sigma_{\scriptscriptstyle-}+\frac{\sqrt{2}}{\sqrt{\pi}}(2q-1)\gamma+\frac{1}{2}\frac{4q^{2}-1}{q\sigma_{\scriptscriptstyle-}}\gamma^{2}+o(\gamma^{2})\\\
&=\frac{2\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}+\frac{\sqrt{\pi}}{\sqrt{2}}\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\gamma+\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}\left(\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{2\left(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}\right)}+1\right)\gamma^{2}+o(\gamma^{2}).\end{split}$
(41)

Finally, note that (40) and (41) are special cases of (31) depending on the
sign of $\gamma$.

###  4.2 ATM implied volatility skew

Recall that $\mathsf{C}_{\rm BS}(\sigma,1,e^{k},T)$ is the BS price (with
volatility $\sigma$) of the cal option with the log-strike $k=\log(K)$ and
maturity $T$. As before, let $\mathsf{C}_{\rm lvm}(1,e^{k},T)$ be the call
price of the same option in the two-valued LVM. Given the log-strike $k$ and
maturity $T$ the corresponding implied volatility $\sigma_{{\small{\rm
BS}}}(T,k)$ is the solution of the following equation for $\sigma$

$\mathsf{C}_{\rm BS}(\sigma,1,e^{k},T)=\mathsf{C}_{\rm lvm}(1,e^{k},T).$

Denote $\partial_{k}=\frac{\partial}{\partial k}$.

######  Theorem 6.

The ATM skew is given by

$\partial_{k}\sigma_{{\small{\rm
BS}}}(T,0)=\frac{1}{\sqrt{T}}\sqrt{\frac{\pi}{2}}e^{\frac{T\sigma_{{\small{\rm
BS}}}^{2}(T,0)}{8}}\left(1-\frac{2\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\left(F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)+F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)\right)\right)$
(42)

and

$\partial_{k}\sigma_{{\small{\rm
BS}}}(T,0)=\frac{1}{\sqrt{T}}\sqrt{\frac{\pi}{2}}\frac{\sigma_{\scriptscriptstyle+}-\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}+o\left(\sqrt{T}\right),\quad\text{as}\quad
T\to 0.$ (43)

######  Remark 4.

Equation (42) is similar to the equation for the ATM skew obtained in Theorem
3.5 in [8]. In particular, the factor
$\frac{1}{\sqrt{T}}\sqrt{\frac{\pi}{2}}e^{\frac{T\sigma_{{\small{\rm
BS}}}^{2}(T,0)}{8}}$ is exactly the same as the one in that theorem. However,
the term
$\left(F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)+F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)\right)$
differs from the similar term in [8]. This difference is caused by the same
reason as the difference in the pricing formulas Note that the short term
asymptotic behaviour of the ATM skew given by (43) is exactly the same as in
[8] (e.g. see Remark 3.2 in that paper).

Before proceeding to the proof of Theorem 6, we prove below two auxiliary
statements. The first one is Proposition 1 that provides an asymptotic result
for the function $F$ defined in (20). The second auxiliary statement is
Proposition 2 that concerns the derivative of the call price with respect to
the log-strike.

######  Proposition 1.

For any fixed $a$ the following holds

$F\left(T,a,0\right)=\frac{a}{\sqrt{2\pi}}\sqrt{T}+\frac{1}{2}+o(\sqrt{T}),\quad\text{as}\quad
T\to 0.$

###### Proof of Proposition 1.

Observe that

$\displaystyle\psi(a,t,0)$ $\displaystyle=\frac{1}{\sqrt{2\pi
t}}+\frac{a}{2}+O\left(\sqrt{t}\right),\quad\text{as}\quad t\to 0,$ (44)
$\displaystyle\phi(t)$ $\displaystyle=\frac{1}{\sqrt{2\pi
t}}+\left(\frac{\sigma_{\scriptscriptstyle+}+\sigma_{\scriptscriptstyle-}}{8\sqrt{2\pi}}-\frac{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}{4}\right)\sqrt{t}+o\left(\sqrt{t}\right),\quad\text{as}\quad
t\to 0.$ (45)

Therefore,

$\displaystyle F\left(T,a,0\right)$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{0}^{T}\frac{e^{-\frac{a^{2}}{2}s}}{\sqrt{T-s}}\left(\frac{a}{2}+\frac{1}{\sqrt{2\pi
s}}\right)ds+o(\sqrt{T})$
$\displaystyle\frac{a}{2\sqrt{2\pi}}\int_{0}^{T}\frac{1}{\sqrt{T-s}}ds+\frac{1}{2\pi}\int_{0}^{T}\frac{1}{\sqrt{s(T-s)}}ds+o(\sqrt{T})$
$\displaystyle=\frac{a}{\sqrt{2\pi}}\sqrt{T}+\frac{1}{2}+o(\sqrt{T}),\quad\text{as}\quad
T\to 0.$

The proof of the proposition is finished. ∎

######  Proposition 2.

We have that

$\partial_{k}\mathsf{C}_{\rm
lvm}(1,e^{k},T)=-2\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},\frac{k}{\sigma_{\scriptscriptstyle+}}\right)e^{k},$
(46)

and

$\begin{split}\partial_{k}\mathsf{C}_{\rm
lvm}(1,1,T)&:=\partial_{k}\mathsf{C}_{\rm
lvm}(1,e^{k},T)|_{k=0}=-2\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)\\\
&=-\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}+\frac{1}{\sqrt{2\pi}}\frac{\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\sqrt{T}+o\left(\sqrt{T}\right),\quad\text{as}\quad
T\to 0.\end{split}$ (47)

###### Proof of Proposition 2.

By Theorem 3 we have that

$\mathsf{C}_{\rm
lvm}\left(1,e^{k},T\right)=2p\left(F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},\frac{k}{\sigma_{\scriptscriptstyle+}}\right)-e^{k}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},\frac{k}{\sigma_{\scriptscriptstyle+}}\right)\right),$

where

$F\left(T,a,\frac{k}{\sigma_{\scriptscriptstyle+}}\right)=\int_{0}^{T}\phi(T-s)\psi\left(a,s,\frac{k}{\sigma_{\scriptscriptstyle+}}\right)e^{-\frac{a^{2}}{2}s}ds,$

and functions $\phi$ and $\psi$ are given by (16) and (18) respectively.
Observe that $\partial_{k}\psi(a,s,k)=0$ for all $a,s$ and $k$. Therefore,

$\displaystyle\partial_{k}\mathsf{C}_{\rm lvm}(1,e^{k},T)$
$\displaystyle=-2p\cdot
e^{k}\int_{0}^{T}\phi(T-s)\psi\left(-\frac{\sigma_{\scriptscriptstyle+}}{2},s,\frac{k}{\sigma_{\scriptscriptstyle+}}\right)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}ds$
$\displaystyle=-2\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},\frac{k}{\sigma_{\scriptscriptstyle+}}\right)e^{k},$

as claimed in (46), and, hence, we have

$\partial_{k}\mathsf{C}_{\rm
lvm}(1,1,T)=-2\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right),$

i.e. the first equation in (47). Applying Proposition 1 with
$a=-\frac{\sigma_{\scriptscriptstyle+}}{2}$ gives the second equation in (47),
the short term asymptotics of $\partial_{k}\mathsf{C}_{\rm lvm}(1,1,T)$, as
claimed. ∎

###### Proof of Theorem 6.

Similarly to $\partial_{k}=\frac{\partial}{\partial k}$, denote
$\partial_{\sigma}=\frac{\partial}{\partial\sigma}$. By the chain rule, we
have that

$\partial_{k}\sigma_{{\small{\rm BS}}}(T,k)=\frac{\partial_{k}\mathsf{C}_{\rm
lvm}(1,e^{k},T)-\partial_{k}\mathsf{C}_{\rm BS}(\sigma_{{\small{\rm
BS}}}(T,k),1,e^{k},T)}{\partial_{\sigma}\mathsf{C}_{\rm
BS}(\sigma_{{\small{\rm BS}}}(T,k),1,e^{k},T)}.$ (48)

Further, observe that

$\displaystyle\partial_{\sigma}\mathsf{C}_{\rm BS}(\sigma_{{\small{\rm
BS}}}(T,0),1,1,T)$ $\displaystyle=\frac{e^{-\frac{\sigma_{{\small{\rm
BS}}}^{2}(T,0)}{8}T}}{\sqrt{2\pi}}\sqrt{T},$ (49)
$\displaystyle\partial_{k}\mathsf{C}_{\rm BS}(\sigma_{{\small{\rm
BS}}}(T,0),1,1,T)$ $\displaystyle=\frac{1}{2}\left(-1+{\rm
Erf}\left(\frac{\sqrt{T}\sigma_{{\small{\rm
BS}}}(T,0)}{2\sqrt{2}}\right)\right),$ (50)

and

${\rm Erf}\left(\frac{\sqrt{T}\sigma_{{\small{\rm
BS}}}(T,0)}{2\sqrt{2}}\right)=\mathsf{C}_{\rm BS}\left(\sigma_{{\small{\rm
BS}}}(T,0),T\right)=\mathsf{V}_{\rm atm}(T).$ (51)

Using (49), (50) and (51) we can rewrite (48) in the case $k=0$ as follows

$\partial_{k}\sigma_{{\small{\rm
BS}}}(T,0)=\frac{\sqrt{2\pi}e^{\frac{\sigma_{{\small{\rm
BS}}}^{2}(0,T)}{8}T}}{\sqrt{T}}\left(\partial_{k}\mathsf{C}_{\rm
lvm}(1,1,T)+\frac{1}{2}-\frac{\mathsf{V}_{\rm atm}(T)}{2}\right)$ (52)

Next, by equation 23,

$\mathsf{V}_{\rm
atm}(T)=\frac{2\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\left(F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)-F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)\right)$

and, by Proposition 2,

$\partial_{k}\mathsf{C}_{\rm
lvm}(1,1,T)=-\frac{2\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right).$

Therefore,

$\partial_{k}\mathsf{C}_{\rm lvm}(1,1,T)+\frac{1}{2}-\frac{\mathsf{V}_{\rm
atm}(T)}{2}=\frac{1}{2}-\frac{\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\left(F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)+F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},0\right)\right)$

and, getting back to (52), we obtain, after simple algebra, that
$\sigma_{{\small{\rm BS}}}(T,0)$ is equal to (42), as claimed.

Finally, equation (43) follows from (42) and Proposition 1 (we skip details).

∎

##  5 Proofs of Theorems 2–5 and Corollary 1

###### Proof of Theorem 2.

By equation (15) we have that

$\begin{split}p(0,x,T)&=2a(x)e^{-\frac{\sigma(x)x}{2}}\int_{u+v+s=T}\left(\int_{0}^{\infty}h(v,lp)h(u,lq)dl\right)h(s,x)e^{-\frac{\sigma_{\rm
p}^{2}v+\sigma_{\rm m}^{2}u+\sigma^{2}(x)s}{8}}dsdv.\end{split}$ (53)

Recall the following two equations that were used in the proof of [3, Theorem
3, Part 1)], namely

$\int\limits_{0}^{\infty}h(v,lp)h(u,lq)dl=\frac{pq}{2\sqrt{2\pi}\left(p^{2}u+q^{2}v\right)^{3/2}},$
(54)

and

$\int_{0}^{w}\frac{pq}{\sqrt{2\pi}\left(p^{2}(w-v)+q^{2}v\right)^{3/2}}e^{-\frac{\sigma_{\rm
p}^{2}v+\sigma^{2}_{\rm m}(w-v)}{8}}dv=\phi(w)\quad\text{for}\quad w>0,$ (55)

where the function $\phi$ is defined in (16). Using these equations in (53)
gives the claimed equation for the density of $X_{T}$, that is

$p(0,x,T)=2a(x)e^{-\frac{\sigma(x)}{2}x}\int_{0}^{T}\phi(T-s)h(s,x)e^{-\frac{\sigma^{2}(x)}{8}s}ds\quad\text{for}\quad
x\in{\mathbb{R}}.$

∎

###### Proof of Theorem 3.

Since $S_{0}=1$ and $K\geq 1$, we have that
$X_{0}=\log(S_{0})/\sigma_{\scriptscriptstyle+}=0$ and
$k=\log(K)/\sigma_{\scriptscriptstyle+}\geq 0$ respectively. By Theorem 2) we
have that

$\displaystyle\mathsf{C}_{\rm lvm}(1,K,T)$
$\displaystyle=\int_{k}^{\infty}\left(e^{\sigma_{\scriptscriptstyle+}x}-e^{\sigma_{\scriptscriptstyle+}k}\right)p(0,x,T)dx$
$\displaystyle=2p\int_{k}^{\infty}\left(e^{\sigma_{\scriptscriptstyle+}x}-e^{\sigma_{\scriptscriptstyle+}k}\right)e^{-\frac{\sigma_{\scriptscriptstyle+}}{2}x}\left(\int_{0}^{T}\phi(T-s)h(s,x)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}ds\right)dx$
$\displaystyle=2p\int_{0}^{T}\phi(T-s)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{2}s}\left(\int_{k}^{\infty}\left[e^{\frac{\sigma_{\scriptscriptstyle+}}{2}x}-e^{\sigma_{\scriptscriptstyle+}k}e^{-\frac{\sigma_{\scriptscriptstyle+}}{2}x}\right]h(s,x)dx\right)ds$
$\displaystyle=2p\int_{0}^{T}\phi(T-s)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}\left[\psi\left(\frac{\sigma_{\scriptscriptstyle+}}{2},s,k\right)-e^{\sigma_{\scriptscriptstyle+}k}\psi\left(-\frac{\sigma_{\scriptscriptstyle+}}{2},s,k\right)\right]ds$
$\displaystyle=2p\left(F\left(T,\frac{\sigma_{\scriptscriptstyle+}}{2},k\right)-e^{\sigma_{\scriptscriptstyle+}k}F\left(T,-\frac{\sigma_{\scriptscriptstyle+}}{2},k\right)\right),$

as claimed.

Equation (22) for the put price $\mathsf{P}_{\rm lvm}(1,K,T)$ can be obtained
in a similar way (we skip the details). ∎

###### Proof of Corollary 1.

Note first that

$\displaystyle\psi(a,s,0)-\psi(-a,s,0)$
$\displaystyle=ae^{\frac{a^{2}}{2}s}\quad\text{for all}\quad a.$ (56)

Combining (56) with equation (23) we obtain that

$\displaystyle\mathsf{V}_{\rm atm}(T)$ $\displaystyle=\mathsf{C}_{\rm
lvm}(1,1,T)$
$\displaystyle=2p\int_{0}^{T}\phi(T-s)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}\left(\psi\left(\frac{\sigma_{\scriptscriptstyle+}}{2},s,0\right)-\psi\left(-\frac{\sigma_{\scriptscriptstyle+}}{2},s,0\right)\right)ds$
$\displaystyle=p\sigma_{\scriptscriptstyle+}\int_{0}^{T}\phi(T-s)ds=p\sigma_{\scriptscriptstyle+}\int_{0}^{T}\phi(s)ds=\frac{\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+}}\int_{0}^{T}\phi(s)ds.$

Further, a direct computation gives that

$\begin{split}\int_{0}^{T}\phi(T-s)ds&=\frac{\sqrt{T}}{\sqrt{2\pi}}\frac{1}{(\sigma_{\scriptscriptstyle-}-\sigma_{\scriptscriptstyle+})}\left(\sigma_{\scriptscriptstyle-}e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}T}-\sigma_{\scriptscriptstyle+}e^{-\frac{\sigma_{\scriptscriptstyle-}^{2}}{8}T}\right)\\\
&+\frac{\sigma_{\scriptscriptstyle-}(4+\sigma_{\scriptscriptstyle+}^{2}T)}{4\sigma_{\scriptscriptstyle+}(\sigma_{\scriptscriptstyle-}-\sigma_{\scriptscriptstyle+})}{\rm
Erf}\left(\sigma_{\scriptscriptstyle+}\sqrt{\frac{T}{8}}\right)-\frac{\sigma_{\scriptscriptstyle+}(4+\sigma_{\scriptscriptstyle-}^{2}T)}{4\sigma_{\scriptscriptstyle-}(\sigma_{\scriptscriptstyle-}-\sigma_{\scriptscriptstyle+})}{\rm
Erf}\left(\sigma_{\scriptscriptstyle-}\sqrt{\frac{T}{8}}\right).\end{split}$
(57)

Combining the above and simplifying gives the ATM price, as claimed. ∎

###### Proof of Theorem 4.

Recall that in the two-valued LVM
$\sigma(t,K)=\sigma_{\scriptscriptstyle+}{\bf 1}_{\\{K\geq
1\\}}+\sigma_{\scriptscriptstyle-}{\bf 1}_{\\{K<1\\}}$. In the case $S_{0}=1$
equation (26) becomes

$\displaystyle\mathsf{C}_{\rm lvm}\left(1,K,T\right)$
$\displaystyle=\frac{K^{2}\sigma_{\scriptscriptstyle+}^{2}}{2}\int_{0}^{T}p_{S}\left(S_{0},K,t\right)dt\quad\text{for}\quad
K\geq 1.$

Noting that
$\displaystyle{p_{S}\left(1,K,t\right)=\frac{p(0,\log(K)\sigma_{\scriptscriptstyle+}^{-1},t)}{K\sigma_{\scriptscriptstyle+}}}$
for $K\geq 1$, where $p(0,\cdot,t)$ is the density of $X_{t}$ given that
$X_{0}=0$, we get that

$\mathsf{C}_{\rm
lvm}\left(1,K,T\right)=\frac{K\sigma_{\scriptscriptstyle+}}{2}\int_{0}^{T}p\left(0,\log(K)\sigma_{\scriptscriptstyle+}^{-1},t\right)dt\quad\text{for}\quad
K\geq 1.$

Using the equation for $p(0,k,t)$ from Lemma 2, we obtain

$\displaystyle\mathsf{C}_{\rm lvm}(1,K,T)$
$\displaystyle=\frac{K\sigma_{\scriptscriptstyle+}}{2}\int_{0}^{T}p(0,\log(K)\sigma_{\scriptscriptstyle+}^{-1},t)dt$
$\displaystyle=\sqrt{K}\sigma_{\scriptscriptstyle+}p\int_{0}^{T}\int_{0}^{t}\phi(t-s)h(s,\log(K)\sigma_{\scriptscriptstyle+}^{-1})e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}dsdt$
$\displaystyle=\sqrt{K}\int_{0}^{T}h(s,\log(K)\sigma_{\scriptscriptstyle+}^{-1})e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}\left[\frac{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}{\sigma_{\scriptscriptstyle+}+\sigma_{\scriptscriptstyle-}}\int_{0}^{T-s}\phi(t)dt\right]ds$
$\displaystyle=\sqrt{K}\int_{0}^{T}h\left(s,\log(K)\sigma_{\scriptscriptstyle+}^{-1}\right)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}s}\mathsf{V}_{\rm
atm}(T-s)ds\quad\text{for}\quad K\geq 1,$

where $\mathsf{V}_{\rm atm}(T-s)$ is the ATM price (Theorem 1) for maturity
$T-s$.

∎

###### Proof of Theorem 5.

We obtain the BS approximation only for the call option price $\mathsf{C}_{\rm
lvm}(1,K,T)$, as the case of the put option is similar.

Using equation (53) for the density of $X_{T}$ gives the following equation
for the price of the call option

$\mathsf{C}_{\rm
lvm}(1,K,T)=2pA(\sigma_{\scriptscriptstyle+},\sigma_{\scriptscriptstyle-}),$

where

$\displaystyle A(\sigma_{\scriptscriptstyle+},\sigma_{\scriptscriptstyle-})$
$\displaystyle=\int_{k}^{\infty}\left(e^{\sigma_{\scriptscriptstyle+}x}-e^{\sigma_{\scriptscriptstyle+}k}\right)e^{-\frac{\sigma_{\scriptscriptstyle+}}{2}x}\int_{0}^{\infty}\int_{0}^{T}\int_{0}^{T-s}h(v,lp)h(T-v-s,lq)h(s,x)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}v+\sigma_{\scriptscriptstyle-}^{2}u+\sigma_{\scriptscriptstyle+}^{2}s}{8}}dvdsdldx.$

For the BS price of the same call option in the case when volatility is equal
to $\sigma_{\scriptscriptstyle+}$ we have that

$\displaystyle\mathsf{C}_{\rm
BS}(\sigma_{\scriptscriptstyle+},1,K,T)=A(\sigma_{\scriptscriptstyle+},\sigma_{\scriptscriptstyle+})$
$\displaystyle=\int_{k}^{\infty}\left(e^{\sigma_{\scriptscriptstyle+}x}-e^{\sigma_{\scriptscriptstyle+}k}\right)e^{-\frac{\sigma_{\scriptscriptstyle+}}{2}x}\int_{0}^{\infty}\int_{0}^{T}\int_{0}^{T-s}h(v,lp)h(T-v-s,lq)h(s,x)e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}T}{8}}dvdsdldx.$

Observe that

$\displaystyle\left|e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}v+\sigma_{\scriptscriptstyle-}^{2}u+\sigma_{\scriptscriptstyle+}^{2}s}{8}}-e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}T}{8}}\right|$
$\displaystyle=\left|1-e^{-\frac{1}{8}(\sigma_{\scriptscriptstyle-}^{2}-\sigma_{\scriptscriptstyle+}^{2})u}\right|e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}T}{8}}\quad\text{for
all}\quad T>0,$

since $v+u+s=T$. Therefore,

$\displaystyle\left|e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}v+\sigma_{\scriptscriptstyle-}^{2}u+\sigma_{\scriptscriptstyle+}^{2}s}{8}}-e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}T}{8}}\right|$
$\displaystyle\leq\frac{\left|\sigma_{\scriptscriptstyle-}^{2}-\sigma_{\scriptscriptstyle+}^{2}\right|T}{8}e^{-\frac{\sigma_{\scriptscriptstyle+}^{2}}{8}T}+o(T)\quad\text{as}\quad
T\to 0,$

and, hence,
$|A(\sigma_{\scriptscriptstyle+},\sigma_{\scriptscriptstyle-})-A(\sigma_{\scriptscriptstyle+},\sigma_{\scriptscriptstyle+})|\leq
C_{1}T$ for some constant $C_{1}$ and all sufficiently small $T$, which in
turn implies that

$|\mathsf{C}_{\rm lvm}(1,K,T)-2p\mathsf{C}_{\rm
BS}(\sigma_{\scriptscriptstyle+},1,K,T)|\leq cT$

for some constant $c$ and all sufficiently small $T$. ∎

######  Remark 5.

It should be noted that BS approximation is not quite visible from the final
form of option pricing formulas established in Theorem 3. The key step in the
proof of Theorem 5 is the use of equation (53) for the density of $X_{T}$
instead of the final equation for the density given in Theorem (2).

##  6 Conclusion

In this paper we consider the LVM with discontinuity, in which volatility
takes two possible values. A particular value depends on the position of the
underlying price with respect to a given threshold. The model is well known,
and a number of results have been obtained for the model in recent years. In
particular, explicit pricing formulas for European options have been obtained
in [8] in the case when the threshold is taken at the money. These formulas
have been then used to establish that the skew explodes as $T^{-1/2}$, as
maturity $T\to 0$, which accurately reproduces the well-known power law
behaviour of the skew observed in real data. The research method in [8] is
based on the Laplace transform and related techniques. In the present paper we
propose another approach to obtaining the same results. The approach is based
on the natural relationship of the LVM under consideration with SBM and
consists of using the joint distribution of SBM and some of its functionals.
This distribution was obtained in [3], where it was also applied to option
valuation under the two-valued LVM in the general case. First of all, using
this method streamlines the computation of option prices. Secondly, the method
allows to obtain the simple approximation of these prices by the corresponding
BS prices in the case of short maturities. Despite its simplicity, this
approximation turns out to be rather precise and allows to establish the
aforementioned asymptotic behaviour of the skew.

##  7 Appendix. SBM with two-valued drift

The process $X_{t}$ defined in (3) is a special case of the stochastic process
$Z_{t}=(Z_{t},\,t\geq 0)$ defined as a strong solution of the equation

$dZ_{t}=m(Z_{t})dt+(p-q)dL^{(Z)}_{t}+dW_{t},$ (58)

where

$m(z)=m_{1}{\bf 1}_{\\{z\geq 0\\}}+m_{2}{\bf 1}_{\\{z<0\\}},$ (59)

$p\geq 0$ and $q\geq 0$ are given constants, such that $p+q=1$, $L_{t}^{(Z)}$
is the local time of the process $Z_{t}$ at the origin (defined similarly to
(8)), and $(W_{t},\,t\geq 0)$ is standard BM. The existence and uniqueness of
the strong solution of the equation is well known (e.g. see [5] and references
therein). In the special case $m_{1}=m_{2}=0$, the process $Z_{t}$ is SBM with
parameter $p\in[0,1]$ (e.g., see [5] and references therein). By analogy, the
process $Z_{t}$ can be called SBM with the two-valued drift $m$. Theorem 1 is
a special case of the following theorem.

######  Theorem 7 ([3], Theorem 2).

Let $\tau^{(Z)},V^{(Z)}$ and $L^{(Z)}_{T}$ be the last visit to the origin,
the occupation time and the local time at the origin of the process $Z_{t}$.
If $Z_{0}=0$, then the joint density of random variables
$\tau^{(Z)},V^{(Z)},Z_{T}$ and $L^{(Z)}_{T}$ is is given by

$\begin{split}f_{T}(t,v,x,l)&=2a(x)h(v,lp)h(t-v,lq)h(T-t,x)e^{-\frac{m_{1}^{2}v+m^{2}_{2}(t-v)+m^{2}(x)(T-t)}{2}-l(pm_{1}-qm_{2})+m(x)x},\\\
&\quad\text{for}\quad 0\leq v\leq t\leq T\quad\text{and}\quad l\geq
0.\end{split}$ (60)

for $(v,u,s):v,u,s\geq 0,\,v+u+s=T$ and $l\geq 0$.

######  Example 1.

In the special case $m_{1}=m_{2}=m$ and $p=1/2$ equation (12) becomes

$f_{T}(t,v,x,l)=h(v,l/2)h(t-v,l/2)h(T-t,x)e^{-\frac{m^{2}T}{2}+mx},$ (61)

for $0\leq v\leq t\leq T$ and $l\geq 0$, which is the joint density of random
variables $\tau^{(Z)},V^{(Z)},Z_{T}$ and $L^{(Z)}_{T}$ of the process
$Z_{t}=mt+W_{t}$ in the case when $Z_{0}=0$.

######  Remark 6.

SDE (4) is a special case of (58) with the drift specified by values
$m_{1}=-\frac{\sigma_{\scriptscriptstyle+}}{2}$ and
$m_{2}=-\frac{\sigma_{\scriptscriptstyle-}}{2}$ and probabilities $p$ and $q$
are given by (6) and (7) respectively. In this case we have that
$pm_{1}-qm_{2}=-\frac{\sigma_{\scriptscriptstyle-}\sigma_{\scriptscriptstyle+}}{2(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+})}+\frac{\sigma_{\scriptscriptstyle+}\sigma_{\scriptscriptstyle-}}{2(\sigma_{\scriptscriptstyle-}+\sigma_{\scriptscriptstyle+})}=0,$
which reduces density (60) to (12).

## References

  * [1] Dupire, B. (1994). Pricing with a smile. Risk, 7, pp. 18–20.
  * [2] Dupire, B. (1997). Pricing and hedging with smiles, in Mathematics of Derivative Securities. M. Dempster and S. Pliska, eds., Cambridge University Press, pp. 103–111.
  * [3] Gairat, A. and Shcherbakov, V. (2017). Density of Skew Brownian motion and its functionals with application in finance. Mathematical Finance, 27(4), pp. 1069–1088.
  * [4] Gairat, A. and Shcherbakov, V. (2022). Skew Brownian motion with dry friction: joint density approach. Statistics and Probability Letters, 187, 109511.
  * [5] Lejay, A. (2006). On the constructions of the skew Brownian motion. Probability Surveys, 3, pp. 413–466.
  * [6] Lipton, A. (2018). Oscillating Bachelier and Black–Scholes formulas. In: Lipton, A. (ed.) Financial Engineering–Selected Works of A. Lipton, pp. 371–394. World Scientific, Singapore.
  * [7] Lipton, A. and Sepp, A. (2011). Filling the gaps. Risk Magazine, pp. 66–71.
  * [8] Pigato, P. (2019). Extreme at the money skew in a local volatility model. Finance and Stochastics, 23(4), pp. 827–859.


# DeepClouds.ai: Deep learning enabled computationally cheap direct numerical
simulations

Moumita Bhowmik
Indian Institute of Tropical Meteorology
Pune, India
<EMAIL_ADDRESS>
Manmeet Singh
Indian Institute of Tropical Meteorology
Pune, India
<EMAIL_ADDRESS>
Suryachandra Rao
Indian Institute of Tropical Meteorology
Pune, India
<EMAIL_ADDRESS>
Souvik Paul
Hewlett Packard Enterprise
Pune, India
<EMAIL_ADDRESS>

###### Abstract

Simulation of turbulent flows, especially at the edges of clouds in the
atmosphere, is an inherently challenging task. Hitherto, the best possible
computational method to perform such experiments is the Direct Numerical
Simulation (DNS). DNS involves solving non-linear partial differential
equations for fluid flows, also known as Navier-Stokes equations, on
discretized grid boxes in a three-dimensional space. It is a valuable paradigm
that has guided the numerical weather prediction models to compute rainfall
formation. However, DNS cannot be performed for large domains of practical
utility to the weather forecast community. Here, we introduce DeepClouds.ai, a
3D-UNET that simulates the outputs of a rising cloud DNS experiment. The
problem of increasing the domain size in DNS is addressed by mapping an inner
3D cube to the complete 3D cube from the output of the DNS discretized grid
simulation. Our approach effectively captures turbulent flow dynamics without
having to solve the complex dynamical core. The baseline shows that the deep
learning-based simulation is comparable to the partial-differential equation-
based model as measured by various score metrics. This framework can be used
to further the science of turbulence and cloud flows by enabling simulations
over large physical domains in the atmosphere. It would lead to cascading
societal benefits by improved weather predictions via advanced
parameterization schemes.

##  1 Introduction

One of the cornerstones of atmospheric science is the study of cloud
microphysics. Cloud prediction is crucial for understanding atmospheric
turbulence and climate model uncertainties. Almost two-thirds of the Earth’s
surface is covered by clouds [1]. By reflecting solar radiation, transporting,
and releasing energy, they have an impact on the Earth-atmosphere system. They
also play a crucial role in regulating the atmospheric hydrological cycle,
modulating atmospheric chemistry, and impacting extreme weather events through
atmospheric moisture redistribution. These phenomena have their origins in
cloud microphysical processes, which have been investigated extensively in
recent decades [2, 3, 4, 5, 6, 7, 8, 9]. The microphysical characteristics of
cloud droplets are strongly influenced by the turbulent flow surrounding them.
Direct numerical simulation (DNS) is a straightforward approach for
disentangling many feedbacks between small-scale turbulence and droplet
dynamics by resolving the turbulence down to the smallest eddy size. Three-
dimensional (3D) DNS combines the Eulerian description of turbulent flow with
the Lagrangian development of an ensemble of cloud water droplets [10].

###  1.1 Critical analysis

DNS can directly solve Navier-Stokes equations without any turbulence model.
However, it has some significant weaknesses, including the need for a very
fine grid in DNS to capture the small-scale turbulent flow and the generation
of big data for turbulent flow. As a result of these factors, DNS is a
computationally expensive technique [11, 12]. In a nutshell, due to existing
computational resource constraints, DNS is not a feasible option for numerical
weather prediction and climate modelling [13, 14, 15]. Thus, there exists a
gap between the smallest scale cloud microphysics and large scale atmospheric
flow [16]. These limitations necessitate the need for alternate approaches to
simulate small-scale flows in the atmosphere.

###  1.2 Deep learning for DNS

Deep learning has revolutionized various fields, particularly geosciences [17,
18, 19] within the framework of computer vision. DNS typically involves
computations over three-dimensional grid boxes similar to the 3D images in the
computer vision community. Recently, a number of studies have attempted the
use of deep neural networks for solving turbulent flows in fluids [20, 21, 22,
23, 24, 25, 26, 27, 28, 29, 30] . Introducing deep learning in turbulence will
help understand the characteristics of clouds surrounding flow patterns. The
knowledge of cloud flow is highly essential for the development activities in
the atmospheric research community. Accurate representation of the smallest-
scale microphysical processes in weather prediction models by deep learning
would enable furthering the skill of numerical weather prediction (Fig. 1).
Computationally cheap DNS offered by deep learning would aid the development
of better parameterization techniques.

###  1.3 Aim & contribution

The purpose of this study is to enable the modelling of small-scale processes
across broad domains. Solving Navier-Stokes equations with small-scale
processes at spatial scales applicable to the weather and climate science
community is not feasible at present due to the computational costs of DNS. We
introduce DeepClouds.ai, which can perform the mapping of a smaller grid
box(A) within the DNS domain to a larger region (B) which encapsulates the
smaller region(A). Thus, A lies within B, and our deep learning model attempts
to estimate B from A. The mapping of DNS outputs from a smaller three-
dimensional cube to a larger three-dimensional cube is a difficult
undertaking.

Previous studies [20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30] have attempted
to simulate turbulence using data-driven approaches. However, none of them
uses a three-dimensional representation of the flow as it occurs in reality.
Further, none of the existing implementations have attempted the complex,
nonlinear, and computationally difficult task of simulating three-dimensional
atmospheric cloud flows. We use a three-dimensional convolutional neural
network, known as a 3D-UNET, to simulate cloud flows over 2x domain from a
region of size 1x. This is a test case of our newly built framework, which can
be applied to different atmospheric simulations. Because our work is a one-of-
a-kind attempt in this direction, it can be used as a baseline to benchmark
datasets created in the future using DeepClouds.ai.

Figure 1: Schematic showing the recursive simulations of large DNS domains
using DeepClouds.ai. The experiments over larger regions with small scale
micro-physics would improve precipitation forecasts in the numerical weather
prediction models.

##  2 Model Description

###  2.1 Flow Equations

In this section, we describe the fluid flow equations for turbulent flow
fields that form the core of DNS simulations used for training the deep
learning based model. The basic equations for turbulent flow fields are as
follows:

$\nabla.{\bf U}=0,~{}$ (1) $\displaystyle\frac{\partial{\bf U}}{\partial
t}+({\bf U}.\nabla){\bf U}=-\frac{1}{\rho_{0}}\nabla p$ $\displaystyle+$
$\displaystyle\nu\nabla^{2}{\bf U}$
$\displaystyle+Be_{z}+f_{LS}(x,t),~{}~{}~{}$ $\frac{\partial{\bf T}}{\partial
t}+{\bf U}.\nabla{\bf T}=k\nabla^{2}{\bf
T}+\frac{L}{c_{p}}C_{d}-\frac{g_{z}}{c_{p}}w_{z},~{}$ (3) $\frac{\partial
q_{v}}{\partial t}+{\bf U}.\nabla q_{v}=D\nabla^{2}q_{v}-C_{d}.~{}$ (4)

Here, U denotes flow velocity, T is flow temperature, and ${q_{v}}$ is vapor
mixing ratio. $B$ and $f_{LS}$ denote buoyancy term, and forcing term,
respectively. $L$ is latent heat of evaporation, $C_{p}$ is specific heat, and
$C_{d}$ is the condensation rate. A detailed description about these
parameters are given in [10, 31].

Pressure variation (P) is calculated from Hydrostatic equation:

$\frac{d{\bf P}}{dt}=-\frac{g{\bf P}w_{z}}{R_{t}{\bf T}}~{}$ (5)

Here, $R_{t}$ is specific gas constant of moist air, and the vector
$\textbf{g}=(0,0,-g)$ includes the gravitational acceleration
$g(=9.81m/s^{2})$. $w_{z}(=0.5m/s)$ is a constant upward velocity of rising
domain.

The condensation rate ($C_{d}$) per unit mass inside each grid cell is
determined by,

$C_{d}(x,t)=\frac{1}{m_{a}}\frac{dm_{l}}{dt}=\frac{4\pi\rho_{l}}{\rho_{0}a^{3}}\sum_{\beta=1}^{\Delta}S(X_{\beta},t)r(\textbf{X},t)~{}$
(6)

Here, $m_{a}$ is the mass of air per grid cell, the sum collects the droplets
inside each grid cells of size $a^{3}$ that surround the grid point $x$ and
$S$ is supersaturation. More details of (Equation 6) can be found in [10, 31].

Table 1: List of constants, reference values and initial simulation parameters Domain Size | $(25.6cm)^{3}$
---|---
Number of Grids | 256 $\times$ 256 $\times$ 256
Droplets number concentration | 343 $(cm)^{3}$
Initial Temperature | 14.30 C
Initial Pressure | 765 hPa
Relative Humidity of environment | 65$\%$
Initial Height | 500 m
Taylor microscale Reynolds number ($Re_{\lambda}$) | 59
Dissipation rate ($\epsilon$) | 1.5 $\times$ $10^{-5}m^{2}s^{-1}$
Kolmogrove Time ($\tau_{\eta}$) | 122 ms

###  2.2 Droplets Equations

The liquid water component is modeled as a Lagrangian ensemble of $N_{p}$
pointlike droplets in the DNS domain. The droplet motion is governed by drag
force and gravity, with the form:

$\frac{d{\bf X}(t)}{dt}={\bf V}({\bf X},t),~{}$ (7) $\frac{d{\bf V}({\bf
X},t)}{dt}=\frac{{\bf U}({\bf X},t)-{\bf V}({\bf X},t)}{\tau_{p}}+{\bf g},~{}$
(8) $\displaystyle\frac{dr({\bf X},t)}{dt}$ $\displaystyle=$
$\displaystyle\frac{K}{r({\bf X},t)+\xi}(S({\bf X},t)+1-(1+$
$\displaystyle\frac{a{r_{d}}^{3}}{{r({\bf
X},t)}^{3}+b{r_{d}}^{3})})exp(\frac{B}{r({\bf X},t)})),~{}$

where, $U(x,t)$ is the flow velocity at droplet’s position $X(t)$, $V(X,t)$ is
droplet’s velocity, $r({\bf X},t)$ is the droplet’s radius, and $\tau_{p}$
($=2\rho_{l}r^{2}/(9\rho_{o}\nu)$ is a finite particle response time. A detail
description about a, b, $K$, $\xi$, $r_{d}$ and $B$ can be found in appendix B
of [32].

###  2.3 Initial Setup

At the beginning of the DNS, cloud droplets of different sizes, containing dry
aerosols are placed randomly in one-third of the domain, and the rest of the
domain is filled with clear air as shown in (Fig. 2). These droplets can move
freely with turbulent flow and always remain inside the domain. In other
words, the domain is treated as adiabatic. The domain starts rising at about
500m from the cloud base. A 3D box of size $(25.6cm)^{3}$ is used for the
direct numerical simulation. The flow inside the domain is periodic in all
three (x,y,z)-directions. The details of the other initial parameters are
given in Table 1.

Figure 2: Conceptual diagram of rising cloud domain simulated by the
Eulearian-Langrangian DNS. Figure 3: Schematic of data-driven simulations of
clear-cloudy turbulent flows. The inner 3D cube of a direct numerical
simulation model is mapped to the larger 3D cube using a UNET algorithm. The
cartoon shows an imaginative view of the clear-cloudy turbulent flow
simulation that is used as a test case in this study.

Figure 4: Evolution of mean squared error loss while training DeepClouds.ai.

.

##  3 The dataset

We generate a labelled training dataset that consists of the precursor as a
smaller 3D cube extracted from a DNS simulation. The label is a larger domain
at the same time step in the DNS experiment.

###  3.1 DNS simulations for generating labelled data

We ran the DNS simulation corresponding to a 256 x 256 x 256 discretized cube
over a domain of size $(25.6cm)^{3}$ for 600 seconds. The number of CPU
processors required to perform the DNS simulation equals 512. It took a
cumulative simulation time of around 30 hours. DNS is a fully parallel model
and was run on an HPC system having computational capacity of 4.0 PFlops. It
includes 3315 dual-socket compute nodes, each with Intel Xeon Broadwell E52695
v4 18-core processors (2.10GHz). Pratyush is a Cray XC40 supercomputer with
liquid cooling. This HPC system has a total memory capacity of 414 TB, each
node has 128GB DDR4 (2667 MHz) memory. All nodes are interconnected using the
high speed, low latency Cray Aries Network with Dragonfly network topology,
which uses per packet adaptive routing.

The DNS model is integrated with a time-step of 0.0005 seconds, and the data
for turbulent flow inside the rising domain is generated every 0.2 seconds.
This dataset corresponds to instantaneous fields comprising four variables.
These fields are the velocity-(u,v,w) components of the flow fields and vapor
mixing ratio ($q_{v}$) at every grid points of the 3D box. We archive a single
file for each unique time step of the model because the size of DNS data is
enormous, comprising of big 3D matrices for the four variables. Thus, the
total number of time steps for extracting DNS output is equal to the number of
saved files. The number of samples for the training datasets comprises the
total number of files saved for the DNS simulation. Hence, there are 3000
samples consisting of 3D cubical matrices available from DNS simulation for
generating labelled datasets.

##  4 Methodology

From the dataset generated by DNS simulation, precursors or inputs
corresponding to the wind and mixing ratio on the inner three-dimensional cube
are extracted as shown in (Fig. 3). The targets or labels are the entire
three-dimensional DNS simulation for the different variables. All the
variables are trained simultaneously using a 3D-UNET to map from the inner
domain of DNS to the entire DNS cube. We first divide the entire dataset into
training (2000), validation(250) and testing(750). The training data is used
to compute the maximum (max) and minimum (min) values across all the 2000
samples for each variable. These max-min values are then used to normalize the
training, validation and test datasets.

The inner cube of the DNS simulation, which corresponds to a size 128 x 128 x
128, is artificially made to 256 x 256 x 256 using three-dimensional cubic
interpolation to facilitate it as an input to the 3D-UNET. The 3D U-Net
architecture is quite similar to the U-Net in terms of its design. It consists
of both the encoding and decoding components. The encoding and decoding parts
have a kernel size of 5 x 5 x 5 with max-pooling corresponding to 2 x 2 x 2.
Strides of two in each dimension are used and the model consists of a
cumulative 1,794,385 parameters. The learning rate used for training is 0.0001
and the output activation function is ’hyperbolic tangent’. 3D-UNET emerged
from the image segmentation community and has shown promising results in
emulating DNS simulations.

We train the 3D-UNET with DNS outputs corresponding to a smaller region as the
precursors and the entire domain as the labels. Learning is performed in
regression mode using multiple 3D convolutions and deconvolutions. TensorFlow
is used for performing the data-driven optimization with the rectified linear
units (RELU) employed for non-linear activation. We monitor the validation
loss and use early stopping while saving the best model on the validation
dataset at train time. The loss is computed as mean squared error and the
training is performed iteratively on each sample due to the I/O constraints of
the large dataset. The server used for running the jobs is Apollo 6500 Gen10
Plus System comprising of 8 NVIDIA A100 Tensor Core GPUs (40GB). The processor
used is Dual AMD EPYC 7402 Processors (24core $@$ 2.8GHz) with 1TB Memory (16
$\times$ 64GB DIMMS) – 3200 MT/s DDR4 Memory. The GPU is able to fit only one
sample per training epoch, and hence the choice of batch size is driven by our
computational resources. It takes around 3 days to complete the training task
on a single NVIDIA A100 GPU.

While the training time may seem large, the benefit of DeepClouds.ai is
realized during testing. We were able to generate the test predictions for 750
samples in under 1 hour using a single NVIDIA A100 GPU, typically taking
around 1 day for the DNS multi-CPU simulations. For the testing, central cubes
of size 128 x 128 x 128 are extracted from the DNS testing data that
DeepClouds.ai has not seen before. Similar to training, the test inputs are
first enlarged to a size 256 x 256 x 256 using cubic interpolation and the
trained 3D-UNET model performs predictions on these inputs. The test
predictions are compared with the outputs from 256 x 256 x 256 discretized
cube of the DNS simulations. Thus the 3D-UNET learns to perform three-
dimensional spatial enhancement by a factor of 2x.

##  5 Results & Discussion

Score | Value
---|---
Root Mean Squared Error | 5.56e-9
Peak-Signal-to-noise-ratio | 131
Linear regression coefficient | 0.67
Table 2: Various metrics comparing the deep learning-based predictions with
the direct numerical simulation outputs. Figure 5: Visualization of vapor
mixing ratio ($q_{v}$) profiles in 3D domain obtained from DNS and 3D-UNET
model simulation at t = 420sec on the test dataset Figure 6: Visualization of
(a) u-component and (b) v-component of 3D velocity profiles obtained from DNS
and 3D-UNET model simulation at t = 400sec on the test dataset

The efficiency of DeepClouds.ai to enlarge the inner 3D region of a cubical
cloud DNS simulation for the mixing ratio is shown in Table 2. As compared to
the absolute numbers from (Fig. 4), the RMSE is 3-4 orders less than the
absolute values of mixing ratio. The data-driven simulations of cloud flows
have a high PSNR of 131 dB, as shown by Table 2. The learned deep learning
model is trained on DNS data of turbulent flow in an upward increasing cloudy
domain, and predictions are made for an unseen test data set.

The vapors mixing ratio ($q_{v}$) is defined as the mass of water vapor to the
mass of dry air. Fig. 4 shows the $q_{v}$ profiles from DNS (Fig. 4a) and the
learned model (Fig. 4b). It is important to note that since the domain is
adiabatic, the mixing ratio will be conserved inside the domain. This
indicates that when the water vapor content in the flow increases, so does the
value of the mixing ratio, and vice versa. The amount of moisture in the flow
is critical for following the growth trajectories of cloud droplets. If the
water vapor amount in the domain increases with time, it indicates that the
cloud droplets are evaporating. Whereas if the water vapor amount in the
domain decreases, it is indicative of the expanding droplets by condensation.
If this process continues, the droplets will form raindrops, and precipitation
will occur. Fig. 4a demonstrates that the mixing ratio is highest around the
domain’s edges and lower near the center in DNS simulations. DeepClouds.ai
(3D-UNET) is able to well capture the same pattern in the predicted
$q_{v}$-profiles (see Fig. 4b). The temperature difference between the flow
and the cloud droplets is the cause of high $q_{v}$ near the boundary.
Droplets evaporate as they approach the domain’s boundaries, and the mass of
water vapor in the flow increases.

Fig. 3 shows a visual comparison of 3D turbulent flow velocity in streamwise
(i.e. u-component) (Fig. 3 (a-b)) and transverse (v-components) (Fig. 3 (c-d))
directions at instantaneous DNS simulation timestep of t= 400 sec. This
visualization aims to examine the 3D turbulence intensity within the cloudy
domain as it rises upward in time. Fig. 3a shows that the flow profile has a
low variance for the u-component relative to the v-component of the 3D
velocity profiles from the DNS and 3D-UNET. It is clear that the 3D-UNET
performs a reasonable job at representing the smoothness of flow in the
streamwise direction.

In comparison to streamwise turbulent flow, transverse turbulent flow is
substantially rougher (Fig. 3c). This means that the flow in this direction
isn’t evenly distributed across the domain. The trained model adequately
captures the flow’s non-uniformity along the v-components (see Fig. 3d). This
non-uniformity is caused by the domain’s upward movement as well as the
existence of cloud droplets within the domain. The surrounding temperature and
pressure of the environment change with each time step as the domain rises,
resulting in a velocity differential in the turbulent flow. Cloud droplets can
travel freely inside the domain, and as a result, they begin to interact with
one another. These interactions produce impediments for the flow field,
resulting in velocity component non-uniformity. It can be seen from both
figures (Fig. 4 and Fig. 3) that our deep learning model performs reasonably
well relative to the direct numerical simulation.

##  6 Conclusions

In this work, we present DeepClouds.ai which is a data-driven approach to
perform direct numerical simulations (DNS) of large domain small-scale
processes in the atmosphere. DNS is the only possible option at present to
simulate the complex, nonlinear turbulent small-scale flow fields in the
atmosphere, particularly that over the physical clouds. It lacks in its
ability to be carried out over large domains due to the high computational
costs. Even the biggest supercomputer with hundreds of thousands of CPU cores
would struggle to perform a DNS over the spatial scales of a real atmospheric
cloud. The inability to perform such simulations limits the present numerical
weather prediction models in their capacity to form raindrops at small-scales,
i.e., the scales at which the processes would occur in reality. Thus, there is
a grave need to perform DNS simulations for large domains.

DeepClouds.ai is a framework that built using the outputs from DNS simulations
enabling the recursive large domain small-scale data-driven simulations. To
the best of our knowledge, this is a first attempt towards this step. Our
results show realistic simulations of upward rising clear-cloudy turbulent
flow. For the deep learning of DNS simulation, we use a 3D-UNET. The datasets
and the framework will act as a baseline for benchmarking deep learning
enabled direct numerical simulations. DeepClouds.ai serves as a link between
the data science community and physicists working on cloud flows. It has the
potential to benefit a variety of societal concerns by enabling improved
skills in numerical weather prediction systems.

Future research will focus on using various architectures of computer vision
algorithms to improve on the baseline established in this study. For enhanced
data-driven DNS, we would additionally incorporate temporal data in addition
to spatial data. The datasets and codes will be made publicly available to
contribute in the advancement of cloud physics and atmospheric science
research.

## References

  * [1] M. B. Baker.  Cloud microphysics and climate.  Scinece, 276, 1997.
  * [2] P.A. Vaillancourt, M.K. Yau, and W.W. Grabowski.  Microscopic approach to cloud droplet growth by condensation. part i: Model description and results without turbulence.  Journal of the Atmospheric Sciences, 58, 2001.
  * [3] P.A. Vaillancourt, M.K. Yau, P. Bartello, and W.W. Grabowski.  Microscopic approach to cloud droplet growth by condensation. part ii: Turbulence, clustering, and condensational growth.  Journal of the Atmospheric Sciences, 58, 2001.
  * [4] H.R. Pruppacher and J.D. Klett.  Microphysics of clouds and precipitation.  Springer, Dordrecht, 2010.
  * [5] R.R. Rogers and M.K. Yau.  A short course in cloud physics.  3rd ed. ed, International series in natural philosophy. Pergamon Press, Oxford; New York., 1989.
  * [6] A.S. Lanotte, A. Seminara, and F. Toschi.  Cloud droplet growth by condensation in homogeneous isotropic turbulence.  Journal of the Atmospheric Sciences, 66, 2009.
  * [7] G. Sardina, F. Picano, L. Brandt, and R. Caballero.  Continuous growth of droplet size variance due to condensation in turbulent clouds.  Journal of the Atmospheric Sciences, 115, 2015.
  * [8] T. Gotoh, T. Suehiro, and I. Saito.  Continuous growth of cloud droplets in cumulus cloud.  New Journal of Physics, 18, 2016.
  * [9] F. Yang, P. Kollias, R.A. Shaw, and A.M. Vogelmann.  Cloud droplet size distribution broadening during diffusional growth: ripening amplified by deactivation and reactivation.  Atmospheric Chemistry and Physics, 18, 2018.
  * [10] B. Kumar, F Janetzko, J. Schumacher, and R Shaw.  Extreme responses of a coupled scalar particle system during turbulent mixing.  New J. Phys., 14:15020, 2012.
  * [11] Lois Thomas, W. Wojciech Grabowski, and Bipin Kumar.  Diffusional growth of cloud droplets in homogeneous isotropic turbulence: Dns, scaled-up dns, and stochastic model.  Atmos. Chem. Phys., 20:9087–9100, 2020.
  * [12] Bipin Kumar, Moumita Bhowmik, and Ravi S. Nanjundiah.  Clouds, microphysical processes and small-scale simulations.  Current Science.
  * [13] C. G. Speziale.  Turbulent flows.  2000\.
  * [14] C. G. Speziale.  Analytical methods for the development of reynolds-stress closures in turbulence.  Annual Review of Fluid Mechanics, 23(1):107–157, 1991.
  * [15] V. Voller and F. Porte-Agel.  Moore’s law and numerical modeling.  Journal of Computational Physics, 179(2):698–703, 2002.
  * [16] B. Marcia Baker and Thomas Peter.  Small-scale cloud processes and climate.  Nature, 451:299––300, 2008.
  * [17] S. L. Brunton, B. R. Noack, and P Koumoutsakos.  Machine learning for fluid mechanics.  Annu. Rev. Fluid Mech, 52, 2020.
  * [18] S. Pandey and J. Schumacher.  Reservoir computing model of two-dimensional turbulent convection.  Phys. Rev. Fluids, 5, 2020.
  * [19] Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, Nuno Carvalhais, et al.  Deep learning and process understanding for data-driven earth system science.  Nature, 566(7743):195–204, 2019.
  * [20] Saakaar Bhatnagar, Yaser Afshar, Shaowu Pan, Karthik Duraisamy, and Shailendra Kaushik.  Prediction of aerodynamic flow fields using convolutional neural networks.  Computational Mechanics, 64(2):525–545, 2019.
  * [21] Nils Thuerey, Konstantin Wei$\beta$enow, Lukas Prantl, and Xiangyu Hu.  Deep learning methods for reynolds-averaged-navier-stokes simulations of airfoil flows.  AIAA Journal, 58(1):25–36, 2020.
  * [22] Nils Thuerey, Philipp Holl, Maximilian Mueller, Patrick Schnell, Felix Trost, and Kiwon Um.  Physics-based deep learning.  arXiv preprint arXiv:2109.05237, 2021.
  * [23] Andrea Beck, David Flad, and Claus-Dieter Munz.  Deep neural networks for data-driven less closure models.  Journal of Computational Physics, 398:108910, 2019.
  * [24] Chenyue Xie, Ke Li, Chao Ma, and Jianchun Wang.  Modeling subgrid-scale force and divergence of heat flux of compressible isotropic turbulence by artificial neural network.  Physical Review Fluids, 4(10):104605, 2019.
  * [25] Chenyue Xie, Jianchun Wang, Hui Li, Minping Wan, and Shiyi Chen.  Spatially multi-scale artificial neural network model for large eddy simulation of compressible isotropic turbulence.  AIP Advances, 10(1):015044, 2020.
  * [26] Kiwon Um, Robert Brand, Yun Fei, Philipp Holl, , and Nils Thuerey.  Learning from differentiable physics to interact with iterative pde-solvers.  In Advances in Neural Information Processing Systems, 2020.
  * [27] Dmitrii Kochkov, A. Jamie Smith, Ayya Alieva, Qing Wang, P. Michael Brenner, and Stephan Hoyer.  Machine learning accelerated computational fluid dynamics.  Proceedings of the National Academy of Sciences, 118(21), 2021.
  * [28] Romit Maulik, Omer San, Adil Rasheed, and Vedula Prakash.  Subgrid modelling for two dimensional turbulence using neural networks.  J. Atmos. Sci., 52:3620–3634, 1995.
  * [29] Florian Heyder, P. Juan Mellado, and J. Schumacher.  Two-dimensional convective boundary layer:numerical analysis and echo state network model.  arXiv preprint arXiv:2108.06195, 2021.
  * [30] Yu Cheng, Marco Giometto, Lin Ling Kauffmann, Pit and, Chen Cao, Cody Zupnick, Harold Li, Qi Li, Ryan Abernathey, and Pierre Gentine.  Deep learning for subgrid-scale turbulence modeling in large-eddy simulations of the atmospheric boundary layer.  arXiv preprint arXiv:1910.12125.
  * [31] B. Kumar, J. Schumacher, and R. A. Shaw.  Lagrangian mixing dynamics at the cloudy-clear air interface.  J. Atmos. Sci., 71(7):2564–2580, 2014.
  * [32] A. V. Korolev.  The influence of supersaturation fluctuations on droplet size spectra formation.  J. Atmos. Sci., 52:3620–3634, 1995.