# Semimartingale and continuous-time Markov chain approximation for rough
stochastic local volatility models ††thanks: The work was supported by
National Natural Science Foundation of China (Grant No. 12071373) and the
Fundamental Research Funds for the Central Universities China (JBK1805001). We
would like to thank Antoine Jacquier and Philipp Schoenbauer for discussions
leading to improvement in the paper.

Jingtang Ma, Wensheng Yang and Zhenyu Cui School of Economic Mathematics,
Southwestern University of Finance and Economics, Chengdu, 611130, P.R. China
(Email: mjt@swufe.edu.cn). School of Economic Mathematics, Southwestern
University of Finance and Economics, Chengdu, 611130, P.R. China (Email:
yangws@swufe.edu.cn). Corresponding author. School of Business, Stevens
Institute of Technology, Hoboken, New Jersey 07030, United States. (Email:
zcui6@stevens.edu).

###### Abstract

Rough volatility models have recently been empirically shown to provide a good
fit to historical volatility time series and implied volatility smiles of SPX
options. They are continuous-time stochastic volatility models, whose
volatility process is driven by a fractional Brownian motion with Hurst
parameter less than half. Due to the challenge that it is neither a
semimartingale nor a Markov process, there is no unified method that not only
applies to all rough volatility models, but also is computationally efficient.
This paper proposes a semimartingale and continuous-time Markov chain (CTMC)
approximation approach for the general class of rough stochastic local
volatility (RSLV) models. In particular, we introduce the perturbed stochastic
local volatility (PSLV) model as the semimartingale approximation for the RSLV
model and establish its existence, uniqueness and Markovian representation. We
propose a fast CTMC algorithm and prove its weak convergence. Numerical
experiments demonstrate the accuracy and high efficiency of the method in
pricing European, barrier and American options. Comparing with existing
literature, a significant reduction in the CPU time to arrive at the same
level of accuracy is observed.

JEL classification: C63, G13

Keywords: Continuous-time Markov chain, rough stochastic local volatility
models, semimartingale approximation, option pricing

##  1 Introduction

Recently, a new class of stochastic volatility model, named the rough
volatility model, was proposed in Gatheral et al. (2018), and has since then
generated significant amount of interests from both academia and industry. The
key insight of this model is to assume that the latent stochastic volatility
process is driven by a fractional Brownian motion, in contrast to a standard
Brownian motion (e.g. in traditional stochastic (local) volatility models).
The trajectories of the volatility process are continuous but exhibit
irregular path properties due to the fractional Brownian motion driver. From
empirical studies, Gatheral et al. (2018) find that the log-volatility
essentially behaves like a fractional Brownian motion with Hurst exponent $H$
of order $0.1$, at any reasonable time scale. Further empirical evidence has
been documented in Fukasawa et al. (2019), where the authors constructed a
quasi-likelihood estimator applied to realized volatility time series, and
confirmed that the Hurst parameter is much smaller than half, i.e. volatility
is indeed rough.

The rough volatility model has enjoyed huge success in reproducing many
stylized facts of historical volatility time series and implied volatility
smiles for SPX options. On one hand, rough volatility models provide
remarkably accurate fit to the shape of implied volatility smiles, and in
particular for at-the-money skew curves. They also reproduce stylized facts
for realized volatilities (El Euch et al., 2019; Livieri et al., 2018). On the
other hand, it is consistent with economic micro-structural models and
naturally emerges from economic agents’ behaviors, as shown in El Euch et al.
(2018). It also has intrinsic connections with Hawkes processes, see Jaisson
and Rosenbaum (2016); Dandapani et al. (2021).

The existence of the fractional kernel forces the variance process to leave
both the semimartingale and Markovian worlds, hence one important yet
challenging problem in the rough volatility research is to find an efficient
and accurate method to evaluate derivative prices, whereas closed-form
formulae are in general not available. On one hand, several analytical
approximation methods have been introduced and studied in Forde and Zhang
(2017); Guennoun et al. (2018); Forde et al. (2021b). On the other hand, the
Monte Carlo simulation of rough volatility model has been studied in Bayer et
al. (2016); Forde et al. (2021a), etc. Due to the memory in the volatility
process, the Monte Carlo simulation of rough volatility process is very time
consuming, and there is recent research on improving its efficiency, see
McCrickerd et al. (2018); Bayer et al. (2020); Ma and Wu (2021). In the
special class of affine rough volatility models, one is able to price options
through Fourier transform based methods by utilizing the characteristic
function via solving fractional Ricatti systems. In particular, El Euch and
Rosenbaum (2019) generalized the classical Heston model to the rough Heston
model and derive the characteristic function of the log asset price. See also
Abi Jaber and El Euch (2019b); Richard et al. (2021); Abi Jaber and El Euch
(2019a) for extensions. Note that some numerical challenges still remain for
solving the fractional Riccatti system in an efficient way, and see Callegaro
et al. (2021) for recent developments along that direction. In general,
without the affine structure, the characteristic function for the rough
volatility model is not available.

In general, the Monte Carlo method has low computational efficiency, and the
method of finding the characteristic function is not suitable for all rough
volatility models outside of the special class of affine rough volatility
models. Inspecting the existing literature111There is a website dedicated to
collecting the most up-to-date literature on rough volatility research:
https://sites.google.com/site/roughvol/home/risks-1 , to the best of authors’
knowledge, there is no method that is not only generally applicable to various
rough volatility models but also has good accuracy and computational
efficiency. This motivates us to search for a general method that is widely
applicable to the broad class of rough stochastic local volatility (RSLV)
models (see equation (1)). Inspired by the recent success on using the
continuous-time Markov chain (CTMC) method in derivatives pricing (see the
survey Cui et al. (2019) and references therein), this paper aims to extend
the applicability of the CTMC method to the realm of rough volatility models
for the first time.

Our method comprises of two main steps. The first step is a novel
semimartingale approximation to the RSLV model, and we obtain the “perturbed
stochastic local volatility (PSLV) model”, which is new to the literature.
This removes the singularity in the kernel function. The second step is the
CTMC approximation to the PSLV model, and we manage to obtain an explicit
formula involving matrix expressions. In the first step above, we provide a
new semimartingale approximation to the general class of RSLV models, which is
of independent theoretical interest. It is important to distinguish this
semimartingale approximation from a recent series of literature on the
Markovian approximation to rough volatility models, which started from Abi
Jaber and El Euch (2019c), Alfonsi and Kebaier, (2021), and see also Harms
(2019) for the case of fractional Brownian motion. The main idea there is to
represent the fractional process as an integral over a family of Ornstein-
Uhlenbeck processes, and then apply numerical discretization (i.e. quadrature)
to the integral. The final result is a $(n+2)$-dimensional stochastic
differential equation system, where $n$ is the number of grid points of
numerical integration. In contrast, our approach is based on the perturbation
idea, which was first introduced in Dung (2011) for the case of fractional
Brownian motions, and we extend it to the case of general RSLV models to
arrive at the PSLV model. The final result is a stochastic differential
equation system, and the stochastic volatility process in this system is a
semimartingale and Markov process. See Remark 2.4.

The contributions of this paper are three-fold:

  1. 1.

This paper extends the traditional stochastic local volatility (SLV) model to
the rough version, and names it the rough stochastic local volatility (RSLV)
model. Using the semimartingale approximation to the RSLV model, a new model
named the “perturbed stochastic local volatility” (PSLV) model, and its
stochastic differential form are obtained. In addition, this paper discusses
the existence, uniqueness, regularity, semimartingale property and the Markov
property of the PSLV model, and also prove that it converges weakly to the
original RSLV model.

  2. 2.

A novel CTMC approximation method is developed, and we express options prices
under the RSLV model in explicit matrix formulae. To the best of authors’
knowledge, this is the first CTMC algorithm designed for the RSLV model.
Theoretical convergence of this algorithm is established. In addition, a fast
algorithm (Algorithm 5.1) for pricing European and barrier option is given.
Compared with the traditional coupled two-dimensional CTMC method (Cui et al.,
2018), the new CTMC algorithm is decoupled. There is a significant improvement
in computer storage space and computing capacity (see Remark 5.1).

  3. 3.

Numerical examples demonstrate the accuracy and high efficiency of our method.
In particular, the method can deliver European and barrier options prices up
to 3 digits of accuracy in 0.18 seconds of CPU time. Note that the method is
universally fast and accurate across all RSLV models, and it is applicable not
only to path-independent options such as European call/put options, but also
to path-dependent options such as barrier options and American options. There
is a significant reduction in the CPU time to arrive at the same level of
accuracy, as compared to benchmark methods in the literature.

The remainder of this paper is organized as follows. Section 2 presents the
new PSLV model, studies its properties, proves its convergence to the RSLV
model and provides its Markovian representation. Section 3 gives the CTMC
approximation and establishes its weak convergence. Section 4 considers the
European, barrier and American options pricing problems under the RSLV model
and gives the explicit matrix expressions for their prices. Numerical
experiments are reported in Section 5. Finally, Section 6 concludes the paper.

##  2 Semimartingale approximation of rough stochastic local volatility models

###  2.1 Rough stochastic local volatility models

We consider the asset price $\\{S_{t}:t\in\mathbb{T}\\}$ with
$\mathbb{T}:=[0,T]$ which is defined on a filtered probability space
$(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$, where
$\mathbb{F}=\\{\mathcal{F}_{t}\\}_{t\geq 0}$ denotes the standard filtration
generated by a two-dimensional $\mathbb{F}-$Brownian motion $(B,B^{\perp})$
and $W=\rho B+\sqrt{1-\rho^{2}}B^{\perp}$ with constant correlation
$\rho\in(-1,1)$. We first recall the classical stochastic local volatility
(SLV) model as follows

$\hbox{SLV:}\quad\left\\{\begin{array}[]{ll}dS_{t}=\mu(S_{t},V_{t})dt+\varphi(V_{t})\nu(S_{t})dW_{t},\\\
dV_{t}=b(V_{t})dt+\sigma(V_{t})dB_{t},\end{array}\right.$

where $\mu(\cdot,\cdot):\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$,
$\varphi(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$,
$\nu(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$,
$b(\cdot):\mathbb{R}\rightarrow\mathbb{R}$,
$\sigma(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$. $W_{t}$ and $B_{t}$ are
two Brownian motions with a constant correlation $\rho\in(-1,1)$, that is,
$E[dW_{t}dB_{t}]=\rho dt$. The SLV model nests several representative models
in the literature as special cases, such as Heston (Heston, 1993), $4/2$
(Grasselli, 2017), Stein-Stein (Stein and Stein, 1991), $3/2$ (Lewis, 2000),
Hull-White (Hull and White, 1987), $\alpha$-Hypergeometric (Da Fonseca and
Martini, 2016), SABR222Stochastic Alpha Beta Rho (Hagan et al., 2002), Heston-
SABR (Van der Stoep et al., 2014), Quadratic SLV (Lipton, 2002), etc.

Recasting the classical SLV model into its rough correspondent, we have the
following definition of the rough stochastic local volatility (RSLV) model.

######  Definition 2.1 (rough stochastic local volatility model)

Under the risk-neutral measure, assume that the underlying asset price $S_{t}$
follows a rough stochastic local volatility (RSLV) model characterized by the
following two-dimensional diffusion system:

$\hbox{RSLV:}\quad\left\\{\begin{array}[]{ll}dS_{t}=\mu(S_{t},V_{t})dt+\varphi(V_{t})\nu(S_{t})dW_{t},\\\
V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left(b(V_{s})ds+\sigma(V_{s})dB_{s}\right),\end{array}\right.$
(1)

where $K(t,s):=\frac{(t-s)^{H-\frac{1}{2}}}{\Gamma(H+1/2)}$ is the fractional
kernel with the Hurst parameter $H\in(0,1/2)$.

In order to ensure the strong existence of continuous solutions to (1), the
following regularity assumption is necessary.

###### Assumption 2.1

$\mu(\cdot,\cdot):\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$,
$\nu(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$,
$b(\cdot):\mathbb{R}\rightarrow\mathbb{R}$,
$\varphi(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$,
$\sigma(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$ are all Lipschitz
continuous functions with linear growth.

###### Proposition 2.1

Under Assumption 2.1, the equation system (1) admits a unique strong
continuous solution. Moreover $V_{t}$ and $S_{t}$ satisfies

$\sup_{t\in\mathbb{T}}\mathbb{E}[|V_{t}|^{p}]<\infty,~{}~{}\sup_{t\in\mathbb{T}}\mathbb{E}[|S_{t}|^{p}]<\infty,~{}~{}p>0.$

In addition, $V$ and $S$ admit Hölder continuous paths on $\mathbb{T}$ of any
order strictly less than $H$.

Proof We refer to Abi Jaber and El Euch (2019a) for the proofs. $\quad\Box$

###### Remark 2.1

This paper focuses on kernel function of the following form
$K(t,s)=\frac{(t-s)^{H-\frac{1}{2}}}{\Gamma(H+1/2)}$. However, under certain
assumptions, the method in this paper can be applied directly to a more
general class of kernel functions. We refer to Abi Jaber and El Euch (2019a)
for discussions on the regularity condition of the kernel function.

It is well known that the fractional kernel forces the variance process to
leave both the semimartingale and Markovian worlds, which makes numerical
approximation procedures a difficult and challenging task in practice. Abi
Jaber and El Euch (2019c) use a multi-factor model to perform a Markovian
approximation to RSLV. We now show an alternative approach to approximate RSLV
by a semimartingale model through a perturbation idea.

###  2.2 Semimartingale approximation

Inspired by Dung (2011), who performs a semimartingale approximation to the
fractional Brownian motion, we approximate the fractional kernel $K(t,s)$ by a
perturbed kernel $K(t-\varepsilon,s)$ with $0<\varepsilon<<1$. This leads to
the following approximation process $V_{t}^{\varepsilon}$ of the variance
process $V_{t}$:

$V_{t}^{\varepsilon}:=V^{\varepsilon}_{0}+\int^{t}_{0}K(t+\varepsilon,s)\left(b(V^{\varepsilon}_{s})ds+\sigma(V^{\varepsilon}_{s})dB_{s}\right),~{}~{}V^{\varepsilon}_{0}=V_{0},$

where

$K(t+\varepsilon,s)=\frac{(t+\varepsilon-s)^{H-\frac{1}{2}}}{\Gamma(H+1/2)}.$

###### Remark 2.2

The original kernel function $K(t,s)$ is singular at the point $s=t$, and the
kernel function $K(t+\varepsilon,s)$ obtained by thesemimartingale
approximation is smooth for any $s\in[0,t]$.

The first lemma establishes the semimartingale property of the process
$V_{t}^{\varepsilon}$.

###### Lemma 2.1

For any $\varepsilon>0$, $V_{t}^{\varepsilon}$ is a
$\mathcal{F}_{t}-$semimartingale with decomposition

$V_{t}^{\varepsilon}=V^{\varepsilon}_{0}+\int^{t}_{0}\left(K(t+\varepsilon,s)b(V^{\varepsilon}_{s})+\psi^{\varepsilon}_{s}\right)ds+\int^{t}_{0}K(s+\varepsilon,s)\sigma(V^{\varepsilon}_{s})dB_{s},$
(2)

where

$\psi^{\varepsilon}_{s}:=\int^{s}_{0}\partial_{1}K(s+\varepsilon,u)\sigma(V^{\varepsilon}_{u})dB_{u},~{}~{}~{}\partial_{1}K(s+\varepsilon,u)=\frac{(H-\frac{1}{2})(s+\varepsilon-u)^{H-\frac{3}{2}}}{\Gamma(H+\frac{1}{2})}.$

Proof The lemma follows from a straightforward application of the stochastic
Fubini theorem (c.f. Theorem 2.2 in Veraar (2012)):

$\displaystyle\int^{t}_{0}\psi^{\varepsilon}_{s}ds$
$\displaystyle=\int^{t}_{0}\int^{s}_{0}\partial_{1}K(s+\varepsilon,u)\sigma(V^{\varepsilon}_{u})dB_{u}ds$
$\displaystyle=\int^{t}_{0}\left(\int^{t}_{u}\partial_{1}K(s+\varepsilon,u)ds\right)\sigma(V^{\varepsilon}_{u})dB_{u}$
$\displaystyle=\int^{t}_{0}\big{(}K(t+\varepsilon,u)-K(u+\varepsilon,u)\big{)}\sigma(V^{\varepsilon}_{u})dB_{u}$
$\displaystyle=\int^{t}_{0}K(t+\varepsilon,s)\sigma(V^{\varepsilon}_{s})dB_{s}-\int^{t}_{0}K(s+\varepsilon,s)\sigma(V^{\varepsilon}_{s})dB_{s}.$

This completes the proof. $\quad\Box$

The second lemma establishes the strong existence and uniqueness of
$V_{t}^{\varepsilon}$.

###### Lemma 2.2

Under Assumption 2.1, for any $\varepsilon>0$, there exists a unique strong
solution $V_{t}^{\varepsilon}$. Moreover $V^{\varepsilon}_{t}$ satisfies

$\sup_{t\in\mathbb{T}}\mathbb{E}[|V^{\varepsilon}_{t}|^{p}]<\infty,~{}~{}p>0,$

and admits Hölder continuous paths on $\mathbb{T}$ of any order strictly less
than $H$.

Proof We first prove the existence of $V^{\varepsilon}_{t}$. Thanks to the
continuity of $K(t+\varepsilon,s)$ and boundedness of
$\partial_{1}K(t+\varepsilon,s)=\frac{(H-1/2)(t+\varepsilon-s)^{H-3/2}}{\Gamma(H+1/2)}$,
the proof of Theorems 3.3 and 3.4 in Abi Jaber and El Euch (2019a) can be
directly applied to prove the existence of $V^{\varepsilon}_{t}$. Now we show
the pathwise uniqueness. Since $V_{t}^{\varepsilon}$ is a semimartingale as
shown in Lemma 2.1 and
$K(t+\varepsilon,t)=\frac{\varepsilon^{H-\frac{1}{2}}}{\Gamma(H+1/2)}<\infty$,
we can use a similar proof as that of Proposition B.3. in Abi Jaber and El
Euch (2019c) to establish the uniqueness of $V_{t}^{\varepsilon}$. This
completes the proof. $\quad\Box$

The next theorem proves the convergence of the semimartingale approximation.

###### Theorem 2.1

Under Assumption 2.1, the process $V^{\varepsilon}_{t}$ converges to $V_{t}$
in $\mathbf{L}^{2}(\Omega,\mathbb{T})$ as $\varepsilon$ tends to $0$,
uniformly in $t\in\mathbb{T}$.

Proof

$\displaystyle E[|V^{\varepsilon}_{t}-V_{t}|^{2}]$ $\displaystyle=$
$\displaystyle
E\Big{[}\Big{|}\int^{t}_{0}K(t+\varepsilon,s)\big{(}b(V^{\varepsilon}_{s})-b(V_{s})\big{)}ds+\int^{t}_{0}\big{(}K(t+\varepsilon,s)-K(t,s)\big{)}b(V_{s})ds$
$\displaystyle+$
$\displaystyle\int^{t}_{0}K(t+\varepsilon,s)\big{)}\big{(}\sigma(V^{\varepsilon}_{s})-\sigma(V_{s})\big{)}dB_{s}+\int^{t}_{0}\big{(}K(t+\varepsilon,s)-K(t,s)\big{)}\sigma(V_{s})dB_{s}\Big{|}^{2}\Big{]}.$

Recalling the power mean inequality: for $k\geq
1,x_{1},x_{2},\dots,x_{\ell}\geq 0,\ell\in\mathbb{N}_{+}$,
$\sum^{l}_{i=1}\frac{x_{i}}{\ell}\leq(\sum^{l}_{i=1}\frac{x^{k}_{i}}{\ell})^{1/k}$
and using Itô isometry and the Cauchy-Schwarz’s inequality, we have

$\displaystyle E[|V^{\varepsilon}_{t}-V_{t}|^{2}]$ $\displaystyle\leq$
$\displaystyle
4E\Big{[}\Big{(}\int^{t}_{0}K(t+\varepsilon,s)\big{(}b(V^{\varepsilon}_{s})-b(V_{s})\big{)}ds\Big{)}^{2}+\Big{(}\int^{t}_{0}\big{(}K(t+\varepsilon,s)-K(t,s)\big{)}b(V_{s})ds\Big{)}^{2}$
$\displaystyle+\int^{t}_{0}K^{2}(t+\varepsilon,s)\big{)}\big{(}\sigma(V^{\varepsilon}_{s})-\sigma(V_{s})\big{)}^{2}ds+\int^{t}_{0}\big{(}K(t+\varepsilon,s)-K(t,s)\big{)}^{2}\sigma^{2}(V_{s})ds\Big{]}$
$\displaystyle\leq$ $\displaystyle
4E\Big{[}\int^{t}_{0}K^{2}(t+\varepsilon,s)\Big{(}t\big{(}b(V^{\varepsilon}_{s})-b(V_{s})\big{)}^{2}+\big{(}\sigma(V^{\varepsilon}_{s})-\sigma(V_{s})\big{)}^{2}\Big{)}ds$
$\displaystyle+\int^{t}_{0}(K(t+\varepsilon,s)-K(t,s)\big{)}^{2}\big{(}tb^{2}(V_{s})+\sigma^{2}(V_{s})\big{)}ds\Big{]}.$

By the conditions that $b$, $\sigma$ are Lipschitz continuous with linear
growth, and the stochastic Fubini theorem, we have

$\displaystyle E[|V^{\varepsilon}_{t}-V_{t}|^{2}]$ $\displaystyle\leq
4\Big{(}C_{1}(1+t)\int^{t}_{0}K^{2}(t+\varepsilon,s)E[|V^{\varepsilon}_{s}-V_{s}|^{2}]ds$
$\displaystyle+\int^{t}_{0}(K(t+\varepsilon,s)-K(t,s)\big{)}^{2}C_{1}(1+t)\big{(}1+E[|V_{s}|]+E[V_{s}^{2}]\big{)}ds\Big{)}.$

Here and throughout this paper, we use $C_{i}$, $i=1,2,\dots$, to represent
positive constants. By Taylor expansion, there is

$K(t+\varepsilon,s)-K(t,s)=\frac{(t+\varepsilon-s)^{H-\frac{1}{2}}-(t-s)^{H-\frac{1}{2}}}{\Gamma(H+1/2)}=\varepsilon\frac{(H-\frac{1}{2})(t+\varepsilon-s)^{H-\frac{3}{2}}}{\Gamma(H+1/2)}+o(\varepsilon),$

and $E[|V_{s}|^{p}]<\infty$ for $p>0$ from Proposition 2.1, then we obtain

$\displaystyle\int^{t}_{0}\Big{[}(K(t+\varepsilon,s)-K(t,s)\big{)}^{2}C(1+t)\big{(}1+E[|V_{s}|]+E[V_{s}^{2}]\big{)}\Big{]}ds$
$\displaystyle\leq
C_{2}\int^{t}_{0}\Big{(}\varepsilon\frac{(H-\frac{1}{2})(t+\varepsilon-s)^{H-\frac{3}{2}}}{\Gamma(H+1/2)}+o(\varepsilon)\Big{)}ds=C_{2}\varepsilon^{H+\frac{1}{2}}+\mathcal{O}(\varepsilon).$

Finally, using Grönwall’s inequality leads to

$\displaystyle E[|V^{\varepsilon}_{t}-V_{t}|^{2}]$ $\displaystyle\leq
4\Big{[}C_{1}(1+t)\int^{t}_{0}K^{2}(t+\varepsilon,s)E[|V^{\varepsilon}_{s}-V_{s}|^{2}]ds+C_{2}\varepsilon^{H+\frac{1}{2}}+\mathcal{O}(\varepsilon)\Big{]}$
$\displaystyle\leq C_{3}\varepsilon^{H+\frac{1}{2}}.$

This completes the proof. $\quad\Box$

Given the existence, uniqueness, semimartingale and convergence properties of
$V_{t}^{\varepsilon}$, we now define the so-called perturbed stochastic local
volatility (PSLV) model $(S^{\varepsilon}_{t},V^{\varepsilon}_{t})$, which
serves as an approximation of $(S_{t},V_{t})$.

######  Definition 2.2 (perturbed stochastic local volatility models)

We define the following stochastic local volatility model
$(S^{\varepsilon}_{t},V^{\varepsilon}_{t})$ with perturbation parameters
$0<\varepsilon<<1$ as the unique strong solution of

$\hbox{PSLV:}\quad\left\\{\begin{array}[]{ll}dS^{\varepsilon}_{t}=\mu(S^{\varepsilon}_{t},V^{\varepsilon}_{t})dt+\varphi(V^{\varepsilon}_{t})\nu(S^{\varepsilon}_{t})dW_{t},\\\
V_{t}^{\varepsilon}=V^{\varepsilon}_{0}+\int^{t}_{0}K(t+\varepsilon,s)\left(b(V^{\varepsilon}_{s})ds+\sigma(V^{\varepsilon}_{s})dB_{s}\right),~{}~{}~{}V_{0}^{\varepsilon}=V_{0},\end{array}\right.$
(3)

under the same filtered probability space
$(\Omega,\mathcal{F},\mathbb{F},\mathbb{P})$ as $(S_{t},V_{t})$ defined by
(1).

###### Proposition 2.2

Under Assumption 2.1, the equation system (3) admits a unique strong
continuous solution.

Proof The strong existence and uniqueness of $V_{t}^{\varepsilon}$ is given by
Lemma 2.2. Moreover, $S_{t}^{\varepsilon}$ satisfies a stochastic differential
equation and there exists an unique strong solution to it under Assumption 2.1
(see e.g., Oksendal (2013)). This completes the proof. $\quad\Box$

As shown in Theorem 2.1, $V_{t}^{\varepsilon}$ converges to $V_{t}$ as
$\varepsilon$ tends to $0$. The next theorem shows the convergence of
$S^{\varepsilon}_{t}$ to $S_{t}$.

###### Theorem 2.2

Under Assumption 2.1, the process $S^{\varepsilon}_{t}$ converges to $S_{t}$
in $\mathbf{L}^{2}(\Omega,\mathbb{T})$ as $\varepsilon$ tends to $0$,
uniformly in $t\in\mathbb{T}$.

Proof The proof of this theorem is similar to Theorem 2.1. Specifically, using
power mean inequality, Itô isometry and the Cauchy-Schwarz’s inequality, we
obtain

$\displaystyle E[|S^{\varepsilon}_{t}-S_{t}|^{2}]$ $\displaystyle\leq
2E\Big{[}\Big{(}\int^{t}_{0}\big{(}\mu(S^{\varepsilon}_{s},V^{\varepsilon}_{s})-\mu(S_{s},V_{s})\big{)}ds\Big{)}^{2}+\int^{t}_{0}\big{(}\varphi(V^{\varepsilon}_{s})\nu(S^{\varepsilon}_{s})-\varphi(V_{s})\nu(S_{s})\big{)}^{2}ds\Big{]}$
$\displaystyle\leq
2E\Big{[}\int^{t}_{0}t\big{(}\mu(S^{\varepsilon}_{s},V^{\varepsilon}_{s})-\mu(S_{s},V_{s})\big{)}^{2}ds+\big{(}\varphi(V^{\varepsilon}_{s})\nu(S^{\varepsilon}_{s})-\varphi(V_{s})\nu(S_{s})\big{)}^{2}ds\Big{]}.$

Then using the condition that $\mu$, $\varphi$ and $\nu$ are Lipschitz
continuous with linear growth, and the stochastic Fubini theorem, we have

$\displaystyle E[|S^{\varepsilon}_{t}-S_{t}|^{2}]$ $\displaystyle\leq
C_{4}(1+t)\Big{(}\int^{t}_{0}E[|S^{\varepsilon}_{s}-S_{s}|^{2}]ds+\int^{t}_{0}E[|V^{\varepsilon}_{s}-V_{s}|^{2}]ds\Big{)}.$

Thanks to Theorem 2.1 and Grönwall’s inequality, we have

$\displaystyle E[|S^{\varepsilon}_{t}-S_{t}|^{2}]$ $\displaystyle\leq
C_{5}\varepsilon^{H+\frac{1}{2}}.$

This completes the proof. $\quad\Box$

###  2.3 Markovian representation of the PSLV model

In the PSLV model (3) obtained by semimartingale approximation, since the
singularity of the integral kernel $K$ at point $0$ is eliminated, we can
prove that $V^{\varepsilon}_{t}$ is Markovian through the following theorem.

###### Theorem 2.3

The stochastic process $V^{\varepsilon}_{t}$ is a Markov process, and
satisfies the following stochastic differential equation:

$dV_{t}^{\varepsilon}=K^{\varepsilon}b(V^{\varepsilon}_{t})dt+\varphi^{\varepsilon}_{t}dt+K^{\varepsilon}\sigma(V^{\varepsilon}_{t})dB_{t},$
(4)

where

$K^{\varepsilon}:=K(t+\varepsilon,t)$

and

$\varphi^{\varepsilon}_{t}:=\int^{t}_{0}\partial_{1}K(t+\varepsilon,u)\big{(}b(V^{\varepsilon}_{t})dt+\sigma(V^{\varepsilon}_{u})dB_{u}\big{)},~{}~{}~{}\partial_{1}K(t+\varepsilon,u)=\frac{(H-\frac{1}{2})(t+\varepsilon-u)^{H-\frac{3}{2}}}{\Gamma(H+\frac{1}{2})}.$

Proof Differentiating equation (2), (4) is obtained directly. Now we show the
Markovian property of $V^{\varepsilon}_{t}$. Integrating form $s$ to $t$, for
$0\leq s\leq t\leq T$, both sides of the equation (4) gives

$\displaystyle V_{t}^{\varepsilon}-V_{s}^{\varepsilon}$ $\displaystyle=$
$\displaystyle\int^{t}_{s}K(u+\varepsilon,u)\big{(}b(V^{\varepsilon}_{u})du+\sigma(V^{\varepsilon}_{u})dB_{u}\big{)}+\int^{t}_{s}\int^{u}_{0}\partial_{1}K(u+\varepsilon,\xi)\big{(}b(V^{\varepsilon}_{\xi})d\xi+\sigma(V^{\varepsilon}_{\xi})dB_{\xi}\big{)}du$
$\displaystyle=$
$\displaystyle\int^{t}_{s}K(u+\varepsilon,u)\big{(}b(V^{\varepsilon}_{u})du+\sigma(V^{\varepsilon}_{u})dB_{u}\big{)}+\int^{t}_{s}\int^{u}_{s}\partial_{1}K(u+\varepsilon,\xi)\big{(}b(V^{\varepsilon}_{\xi})d\xi+\sigma(V^{\varepsilon}_{\xi})dB_{\xi}\big{)}du$
$\displaystyle+\int^{t}_{s}\int^{s}_{0}\partial_{1}K(u+\varepsilon,\xi)\big{(}b(V^{\varepsilon}_{\xi})d\xi+\sigma(V^{\varepsilon}_{\xi})dB_{\xi}\big{)}du$
$\displaystyle=$
$\displaystyle\int^{t}_{s}K(u+\varepsilon,u)\big{(}b(V^{\varepsilon}_{u})du+\sigma(V^{\varepsilon}_{u})dB_{u}\big{)}+\int^{t}_{s}\int^{u}_{s}\partial_{1}K(u+\varepsilon,\xi)\big{(}b(V^{\varepsilon}_{\xi})d\xi+\sigma(V^{\varepsilon}_{\xi})dB_{\xi}\big{)}du$
$\displaystyle+\int^{t}_{s}\big{(}K(s+\varepsilon,\xi)-K(\varepsilon,\xi)\big{)}\big{(}b(V^{\varepsilon}_{\xi})d\xi+\sigma(V^{\varepsilon}_{\xi})dB_{\xi}\big{)}.$

Thanks to the zero-mean property of the Itô integral and Fubini theorem, there
is

$\displaystyle E[V^{\varepsilon}_{t}|\mathcal{F}_{s}]=$ $\displaystyle
E\bigg{[}V_{s}^{\varepsilon}+\int^{t}_{s}\bigg{(}K(u+\varepsilon,u)b(V^{\varepsilon}_{u})+\int^{u}_{s}\partial_{1}K(u+\varepsilon,\xi)b(V^{\varepsilon}_{\xi})d\xi$
$\displaystyle+\big{(}K(s+\varepsilon,u)-K(\varepsilon,u)\big{)}b(V^{\varepsilon}_{u})\bigg{)}du\bigg{|}\mathcal{F}_{s}\bigg{]}.$

According to Theorem 17.2.3 in Cohen and Elliott (2015), we have

$E[V^{\varepsilon}_{t}|\mathcal{F}_{s}]=E[V^{\varepsilon}_{t}|V^{\varepsilon}_{s}],~{}~{}\text{for}~{}~{}0\leq
s\leq t\leq T,$

and it follows that $V_{t}^{\varepsilon}$ is a Markov process for
$t\in\mathbb{T}$. This completes the proof. $\quad\Box$

###### Remark 2.3

It is well known that $V_{t}$ in the RSLV model (1) is not a Markov process
and does not have Itô differential expression. The reason is that the integral
kernel satisfies $K(t,t)=\infty$. However, after the semimartingale
approximation, $V^{\varepsilon}_{t}$ in the PSLV model (3), which has a smooth
kernel and
$K(t+\varepsilon,t)=\frac{\varepsilon^{H-\frac{1}{2}}}{\Gamma(H+1/2)}$ for a
fixed $0<\varepsilon<<1$, is a semimartingale, and it also has Ito
differential expressions and Markov property.

Although (4) shows the Itô differential form of $V^{\varepsilon}_{t}$, it is
not conducive to calculation and simulation because the term
$\psi^{\varepsilon}_{t}$ is still in the Itô integral form. Next, we consider
another differential expression. Inspired by Abi Jaber and El Euch (2019c),
the perturbed fractional kernel
$K(t+\varepsilon,s)=\frac{(t+\varepsilon-s)^{H-\frac{1}{2}}}{\Gamma(H+1/2)}$
can be written as a Laplace transform of a positive measure $m$:

$K(t+\varepsilon,s)=\int^{\infty}_{0}e^{-\gamma(t+\varepsilon-s)}m(d\gamma),~{}~{}~{}~{}m(d\gamma)=\frac{\gamma^{-H-\frac{1}{2}}d\gamma}{\Gamma(H+1/2)\Gamma(1/2-H)}.$

Then by the stochastic Fubini theorem, we obtain that

$V_{t}^{\varepsilon}=V^{\varepsilon}_{0}+\int^{\infty}_{0}e^{-\gamma\varepsilon}V_{t}^{\varepsilon,\gamma}m(d\gamma),$
(5)

where

$V_{t}^{\varepsilon,\gamma}:=\int^{t}_{0}e^{-\gamma(t-s)}\left(b(V^{\varepsilon}_{s})ds+\sigma(V^{\varepsilon}_{s})dB_{s}\right).$
(6)

###### Theorem 2.4

The PSLV model (3) can be expressed as the following system stochastic
differential equation:

$\left\\{\begin{array}[]{ll}dS^{\varepsilon}_{t}=\mu(S^{\varepsilon}_{t},V^{\varepsilon}_{t})dt+\varphi(V^{\varepsilon}_{t})\nu(S^{\varepsilon}_{t})dW_{t},\\\
dV^{\varepsilon}_{t}=-\left(\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
V_{t}^{\varepsilon,\gamma}m(d\gamma)\right)dt+K^{\varepsilon}b(V^{\varepsilon}_{t})dt+K^{\varepsilon}\sigma(V^{\varepsilon}_{t})dB_{t},\\\
dV_{t}^{\varepsilon,\gamma}=\left(-\gamma
V_{t}^{\varepsilon,\gamma}+b(V^{\varepsilon}_{t})\right)dt+\sigma(V^{\varepsilon}_{t})dB_{t},\end{array}\right.$
(7)

Proof By Itô’s lemma, for fixed $\gamma$, we have

$\displaystyle dV_{t}^{\varepsilon,\gamma}$
$\displaystyle=d\left(\int^{t}_{0}e^{-\gamma(t-s)}\left(b(V^{\varepsilon}_{s})ds+\sigma(V^{\varepsilon}_{s})dB_{s}\right)\right)$
$\displaystyle=-\gamma\left(\int^{t}_{0}e^{-\gamma(t-s)}\left(b(V^{\varepsilon}_{s})ds+\sigma(V_{s})dB_{s}\right)\right)dt+b(V^{\varepsilon}_{t})dt+\sigma(V^{\varepsilon}_{t})dB_{t}$
$\displaystyle=\left(-\gamma
V_{t}^{\varepsilon,\gamma}+b(V^{\varepsilon}_{t})\right)dt+\sigma(V^{\varepsilon}_{t})dB_{t}.$

From the definition of $V_{t}^{\varepsilon,\gamma}$ (6), we have

$\displaystyle E[|V_{t}^{\varepsilon,\gamma}|]$
$\displaystyle=E\left[\left|\int^{t}_{0}e^{-\gamma(t-s)}\left(b(V^{\varepsilon}_{s})ds+\sigma(V^{\varepsilon}_{s})dB_{s}\right)\right|\right]$
$\displaystyle\leq
E\left[\int^{t}_{0}e^{-\gamma(t-s)}\big{(}\left|b(V^{\varepsilon}_{s})\right|ds+\left|\sigma(V^{\varepsilon}_{s})\right|dB_{s}\big{)}\right].$

Since $b(\cdot):\mathbb{R}\rightarrow\mathbb{R}$ and
$\sigma(\cdot):\mathbb{R}\rightarrow\mathbb{R}_{+}$ are all Lipschitz
continuous functions with linear growth and
$\sup\limits_{t\in\mathbb{T}}\mathbb{E}[|V^{\varepsilon}_{t}|^{p}]<\infty,~{}~{}p>0$,
then by Fubini theorem, we have that

$\displaystyle E[|V_{t}^{\varepsilon,\gamma}|]$ $\displaystyle\leq
C\int^{t}_{0}e^{-\gamma(t-s)}ds=C\frac{1-e^{-\gamma t}}{\gamma}.$

We use the Lebesgue dominated convergence theorem to rewrite (5) as

$\displaystyle dV_{t}^{\varepsilon}$
$\displaystyle=\int^{\infty}_{0}e^{-\gamma\varepsilon}dV_{t}^{\varepsilon,\gamma}m(d\gamma)$
$\displaystyle=-\left(\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
V_{t}^{\varepsilon,\gamma}m(d\gamma)\right)dt+\int^{\infty}_{0}e^{-\gamma\varepsilon}m(d\gamma)b(V^{\varepsilon}_{t})dt+\int^{\infty}_{0}e^{-\gamma\varepsilon}m(d\gamma)\sigma(V^{\varepsilon}_{t})dB_{t}$
$\displaystyle=-\left(\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
V_{t}^{\varepsilon,\gamma}m(d\gamma)\right)dt+K^{\varepsilon}b(V^{\varepsilon}_{t})dt+K^{\varepsilon}\sigma(V^{\varepsilon}_{t})dB_{t}.$

Note that by comparing with the formula (4), it can be seen that the first
term in the right of above formula is actually equal to
$\varphi^{\varepsilon}_{t}$ in (4) by using Laplace transform to
$\partial_{1}K(t+\varepsilon,u)$. This completes the proof. $\quad\Box$

###### Remark 2.4

Based on the multifactor approximation, similar Itô differential expressions
with Markov properties can be obtained under the rough Heston model (see
formula (1.4) in Abi Jaber and El Euch (2019c)). It is worth noting that their
approximation method obtains an $(n+2)$-dimensional model, where $n$ is the
number of grid points of numerical integration and the multifactor model
converges to the original model as $n$ tends to infinity. In contrast, the
approximate process $V_{t}^{\varepsilon}$ with Markov property can be obtained
by the semimartingale approximation without numerical integration, thereby
avoiding the difficulties caused by multi-dimensional problems in simulation
and calculation.

Next we consider how to use the Markov chain approximation methods to solve
this stochastic differential equation system.

##  3 CTMC approximation

In this section, we use the CTMC method introduced in Mijatović and Pistorius
(2013) to approximate $(S^{\varepsilon}_{t},V^{\varepsilon}_{t})$ defined in
(3) by a continuous-time Markov chain. To simplify the analysis, we first
decouple the correlation between the two driving Brownian motions by
introducing an auxiliary process $X^{\varepsilon}_{t}$ in the following lemma.

###### Lemma 3.1

Define the functions $g(x):=\int_{\cdot}^{x}\frac{1}{\nu(u)}du$ and
$f(x):=\int_{\cdot}^{x}\frac{\varphi(u)}{K^{\varepsilon}\sigma(u)}du$. Then
the dynamics in (3) can be rewritten as

$\left\\{\begin{array}[]{ll}dX^{\varepsilon}_{t}=\theta(X^{\varepsilon}_{t},V^{\varepsilon}_{t})dt+\sqrt{1-\rho^{2}}\varphi(V^{\varepsilon}_{t})dB^{\bot}_{t},\\\
dV^{\varepsilon}_{t}=-\left(\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
V_{t}^{\varepsilon,\gamma}m(d\gamma)\right)dt+K^{\varepsilon}b(V^{\varepsilon}_{t})dt+K^{\varepsilon}\sigma(V^{\varepsilon}_{t})dB_{t},\\\
dV_{t}^{\varepsilon,\gamma}=\left(-\gamma
V_{t}^{\varepsilon,\gamma}+b(V^{\varepsilon}_{t})\right)dt+\sigma(V^{\varepsilon}_{t})dB_{t},\end{array}\right.$
(8)

where the dynamics of the auxiliary process
$X^{\varepsilon}_{t}:=g(S^{\varepsilon}_{t})-\rho f(v^{\varepsilon}_{t})$ and
the standard Brownian motion $B^{\bot}_{t}:=\frac{W_{t}-\rho
B_{t}}{\sqrt{1-\rho^{2}}}$ is independent from $B_{t}$ with a constant
correlation $\rho\in(-1,1)$. Here

$\displaystyle\theta(X^{\varepsilon}_{t},V^{\varepsilon}_{t})$
$\displaystyle:=\frac{\mu(S^{\varepsilon}_{t},V^{\varepsilon}_{t})}{\nu(S^{\varepsilon}_{t})}-\frac{\nu^{\prime}(S^{\varepsilon}_{t})}{2}\varphi^{2}(V^{\varepsilon}_{t})-\frac{\rho}{2}K^{\varepsilon}\big{(}\sigma(V^{\varepsilon}_{t})\varphi^{\prime}(V^{\varepsilon}_{t})-\sigma^{\prime}(V^{\varepsilon}_{t})\varphi(V^{\varepsilon}_{t})\big{)}$
(9) $\displaystyle-\rho\left(-\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
V_{t}^{\varepsilon,\gamma}m(d\gamma)+K^{\varepsilon}b(V^{\varepsilon}_{t})\right)\frac{\varphi(V^{\varepsilon}_{t})}{K^{\varepsilon}\sigma(V^{\varepsilon}_{t})}.$

Proof This proof follows similarly from Lemma 1 in Cui et al. (2018).
$\quad\Box$

###  3.1 The construction of the CTMC approximation

We first recall the basic setup of CTMC. Denote

$\displaystyle\mathcal{M}$
$\displaystyle:=\\{1,2,\ldots,M\\},~{}~{}~{}~{}\mathcal{M}^{\circ}:=\\{2,\ldots,M-1\\},$
$\displaystyle\mathcal{N}$
$\displaystyle:=\\{1,2,\ldots,N\\},~{}~{}~{}~{}~{}\mathcal{N}^{\circ}:=\\{2,\ldots,N-1\\}.$

Recall that a stochastic process $\alpha(t)$ taking values in the set
$\mathcal{M}$ of M possible states is a CTMC if the distribution of
$\alpha(t+\Delta t)$, conditioned on the current state and the past history up
to time $t$, depends only on the current state $\alpha(t)$. The transition
dynamics of $\alpha(t)$ are characterized by the rate matrix
${\bm{Q}}=(q_{ij})_{M\times M}$, whose elements $q_{ij}$ satisfy (i)
$q_{ii}\leq 0$, and $q_{ij}\geq 0$, if $i\neq j$, and (ii) $\sum_{j}q_{ij}=0$,
$\forall i\in\mathcal{M}$.

In terms of $q_{ij}$, for a time increment $\Delta t\geq 0$, the transition
probability matrix ${\bf P}(\Delta t)$ has the following matrix exponential
representation:

${\bf P}(\Delta t)=\exp({\bm{Q}}\cdot\Delta
t):=\sum_{i=0}^{\infty}\frac{({\bm{Q}}\cdot\Delta t)^{k}}{k!},$

with elements $p_{ij}(\Delta t):=\hbox{Prob}\left\\{\alpha(t+\Delta
t)=j\mid\alpha(t)=i\right\\}$.

We first derive the CTMC approximation
$\widetilde{V}^{\varepsilon}_{t}\in\\{v^{\varepsilon}_{1},v^{\varepsilon}_{2},\dots,v^{\varepsilon}_{M}\\}$
of the variance process $V^{\varepsilon}_{t}$ over a general non-uniform grid
$\\{v^{\varepsilon}_{i}\\}_{i=1}^{M}$ for $i\in\mathcal{M}$. Recall (5) and
let

$\widetilde{V}_{t}^{\varepsilon}=\widetilde{V}^{\varepsilon}_{0}+\int^{\infty}_{0}e^{-\gamma\varepsilon}\widetilde{V}_{t}^{\varepsilon,\gamma}m(d\gamma),$

where $\widetilde{V}^{\varepsilon}_{0}:=V^{\varepsilon}_{0}$,
$\widetilde{V}_{t}^{\varepsilon,\gamma}\in\\{v^{\varepsilon,\gamma}_{i}\\}_{i=1}^{M}$
with
$v^{\varepsilon,\gamma}_{i}=\frac{v^{\varepsilon}_{i}-V^{\varepsilon}_{0}}{\int^{\infty}_{0}e^{-\gamma\varepsilon}e^{-\gamma}m(d\gamma)}e^{-\gamma}=:\frac{v^{\varepsilon}_{i}-V^{\varepsilon}_{0}}{R}e^{-\gamma}$
and

$\hbox{Prob}\left\\{\widetilde{V}_{t}^{\varepsilon,\gamma}=v^{\varepsilon,\gamma}_{i}\big{|}\widetilde{V}_{t}^{\varepsilon}=v^{\varepsilon}_{i}\right\\}=1.$

###### Remark 3.1

The relationship between $V_{t}^{\varepsilon,\gamma}$ and
$V_{t}^{\varepsilon}$ is given by the equation (5). Note that the solution of
the integral equation (5) exists and is unique (see Abi Jaber and El Euch
(2019c)). Thus there is a one-to-one correspondence between
$V_{t}^{\varepsilon}$ and $V_{t}^{\varepsilon,\gamma}$. Recall that
$V_{t}^{\varepsilon}$ is Markov, hence we can use a CTMC to approximate it,
and the corresponding finite state space is ${v_{i}^{\varepsilon}}$. Note that
the intermediate auxiliary process $V_{t}^{\varepsilon,\gamma}$ is not Markov,
and we are not constructing a CTMC approximation to it. Hence
$v^{\varepsilon,\gamma}_{i}$ should not be interpreted as the grid
corresponding to a CTMC, but it is rather solved from the integral equation
(5) when a value $v_{i}^{\varepsilon}$ is substituted into that equation.
Moreover, the variable $\gamma$ in $V_{t}^{\varepsilon,\gamma}$ (defined by
(6)) appears in the exponential form, which is the reason why
$v^{\varepsilon,\gamma}_{i}$ is set to the form
$\frac{v^{\varepsilon}_{i}-V^{\varepsilon}_{0}}{R}e^{-\gamma}$. As a by-
product, for a given grid of the CTMC approximating $V_{t}^{\varepsilon}$, we
have a uniquely defined corresponding value for $V_{t}^{\varepsilon,\gamma}$.
Knowing this one-to-one link between a realization of the Markov process
$V_{t}^{\varepsilon}$ and the non-Markovian process
$V_{t}^{\varepsilon,\gamma}$ is important, and is crucial for the design of
the CTMC approximation to $V_{t}^{\varepsilon}$. Recall from (7) that the
drift term of $V_{t}^{\varepsilon}$ contains an integral with respect to
$V_{t}^{\varepsilon,\gamma}$. Based on the above link, when we carry out
moment matching to construct the CTMC approximation to $V_{t}^{\varepsilon}$,
we can actually express the drift term as an explicit function of
$v_{i}^{\varepsilon}$ and separate the integral with respect to $\gamma$ into
a separate constant term $R$. This is the key advantage of the CTMC method as
we avoid the discretization of the integral with respect to $\gamma$ through a
quadrature method, and this fact precisely leads to a dimension reduction.
Essentially, the process $V_{t}^{\varepsilon,\gamma}$ is just an intermediate
auxiliary process that is uniquely characterized through the integral equation
(5). It does not need to be Markov, and the property of this intermediate
process does not affect our construction of the CTMC approximation. We
construct the CTMC approximation only to the $X_{t}^{\varepsilon}$ and
$V_{t}^{\varepsilon}$, but not $V_{t}^{\varepsilon,\gamma}.$ To sum up, we use
a CTMC $\widetilde{V}^{\varepsilon}_{t}$ to approximate $V_{t}^{\varepsilon}$
based on its Markov property. We first establish the grid points
$v_{i}^{\varepsilon}$, and then use the integral equation (5) to solve
$v_{i}^{\varepsilon,\gamma}$ in terms of $v_{i}^{\varepsilon}$, and finally
substitute it into the drift term of the SDE (7) of $V_{t}^{\varepsilon}$ and
set up moment matching equations to obtain the generator matrix of
$\widetilde{V}^{\varepsilon}_{t}$.

According to Mijatović and Pistorius (2013) and the SDE of
$V_{t}^{\varepsilon}$ in (7), the elements of the tridiagonal generator matrix
$\mathbf{Q}=(q_{ij})_{M\times M}$ of $\widetilde{V}^{\varepsilon}_{t}$ for
$i\in\mathcal{M}^{\circ}$, $j\in\mathcal{M}$, are uniquely determined through
the following system of equations:

$\left[\begin{array}[]{ccc}1&1&1\\\ -h_{i}&0&h_{i+1}\\\
h^{2}_{i}&0&h^{2}_{i+1}\end{array}\right]\left[\begin{array}[]{c}q_{i,i-1}\\\
q_{i,i}\\\ q_{i,i+1}\end{array}\right]=\left[\begin{array}[]{c}0\\\
(v^{\varepsilon}_{i}-V^{\varepsilon}_{0})\widehat{R}+K^{\varepsilon}b(v^{\varepsilon}_{i})\\\
(K^{\varepsilon})^{2}\sigma^{2}(v^{\varepsilon}_{i})\end{array}\right],$ (10)

where $h_{i}=v^{\varepsilon}_{i}-v^{\varepsilon}_{i-1}$,
$\widehat{R}:=-\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
e^{-\gamma}m(d\gamma)/R$, $q_{1,j}=q_{M,j}=0$, $q_{i,j}=0$ for $|i-j|>1$,
$i,j\in\mathcal{M}$. Solving (10) gives

$\left\\{\begin{array}[]{l}q_{i,i-1}=\frac{(K^{\varepsilon})^{2}\sigma^{2}(v^{\varepsilon}_{i})-((v^{\varepsilon}_{i}-V^{\varepsilon}_{0})\widehat{R}+K^{\varepsilon}b(v^{\varepsilon}_{i}))h_{i+1}}{h_{i}(h_{i+1}+h_{i})},\\\
q_{i,i}=\frac{-(K^{\varepsilon})^{2}\sigma^{2}(v^{\varepsilon}_{i})+((v^{\varepsilon}_{i}-V^{\varepsilon}_{0})\widehat{R}+K^{\varepsilon}b(v^{\varepsilon}_{i}))(h_{i+1}-h_{i})}{h_{i+1}h_{i}},\\\
q_{i,i+1}=\frac{(K^{\varepsilon})^{2}\sigma^{2}(v^{\varepsilon}_{i})+((v^{\varepsilon}_{i}-V^{\varepsilon}_{0})\widehat{R}+K^{\varepsilon}b(v^{\varepsilon}_{i}))h_{i}}{h_{i+1}(h_{i+1}+h_{i})}.\end{array}\right.$
(11)

Thus the Markov process $V^{\varepsilon}_{t}$ from (8) is approximated by a
continuous-time Markov chain $\widetilde{V}^{\varepsilon}_{t}$ with the
generator matrix ${\bm{Q}}=(q_{i,j})_{M\times M}$, whose entries are given in
(11).

###### Remark 3.2

It is worth noting that the generator matrix ${\bm{Q}}=(q_{i,j})_{M\times M}$
of $\widetilde{V}^{\varepsilon}_{t}$ is well-defined and time-homogeneous.
Each element $q_{i,j}$ shown in (11) is explicitly expressed, where
$v^{\varepsilon}_{i}$ and $h_{i}$ are provided by grids design, $R$ and
$\widehat{R}$ are both constants, and
$K^{\varepsilon}=\varepsilon^{H-\frac{1}{2}}/\Gamma(H+1/2)$. Therefore, the
CTMC method is fully explicit and very computational friendly.

Next, we derive the CTMC approximation
$\widetilde{X}^{\varepsilon}_{t}\in\\{x^{\varepsilon}_{1},x^{\varepsilon}_{2},\dots,x^{\varepsilon}_{N}\\}$
of the auxiliary process $X^{\varepsilon}_{t}$ over a general non-uniform grid
$\\{x^{\varepsilon}_{i}\\}_{i=1}^{N}$, where
$x^{\varepsilon}_{i}=\sum^{i}_{\ell=1}\delta_{\ell}$ for $i\in\mathcal{N}$,
and $\delta_{\ell}$ is the grid interval. After approximating the variance
process $V^{\varepsilon}_{t}$ by the CTMC $\widetilde{V}^{\varepsilon}_{t}$,
the auxiliary process becomes a nonlinear regime-switching diffusion and its
parameters have $M$ states. In this way, we can use the technique introduced
in Cui et al. (2018) to approximate $X^{\varepsilon}_{t}$ by a continuous-time
Markov chain $\widetilde{X}^{\varepsilon}_{t}$. In particular, according to
(8), for each $\ell\in\mathcal{M}$, we define a second-layer Markov chain
approximation which is determined by the rate matrix
${\bf\Lambda}_{\ell}=(\lambda^{\ell}_{ij})_{N\times N}$ with

$\left[\begin{array}[]{ccc}1&1&1\\\ -\delta_{i}&0&\delta_{i+1}\\\
\delta_{i}^{2}&0&\delta_{i+1}^{2}\end{array}\right]\left[\begin{array}[]{c}\lambda^{\ell}_{i,i-1}\\\
\lambda^{\ell}_{i,i}\\\
\lambda^{\ell}_{i,i+1}\end{array}\right]=\left[\begin{array}[]{c}0\\\
\theta(x^{\varepsilon}_{i},v^{\varepsilon}_{l})\\\
(1-\rho^{2})\varphi^{2}(v^{\varepsilon}_{l})\end{array}\right],$ (12)

where $\lambda^{\ell}_{1i}=\lambda^{\ell}_{Ni}=0$ and $\lambda^{\ell}_{i,j}=0$
for $|i-j|>1$, $i,j\in\mathcal{N}$ and $\ell\in\mathcal{M}$. Solving (12)
gives

$\left\\{\begin{array}[]{l}\lambda^{\ell}_{i,i-1}=\frac{(1-\rho^{2})\varphi^{2}(v^{\varepsilon}_{\ell})-\theta(x^{\varepsilon}_{i},v^{\varepsilon}_{l})\delta_{i+1}}{\delta_{i}(\delta_{i+1}+\delta_{i})},\\\
\lambda^{\ell}_{i,i}=\frac{-(1-\rho^{2})\varphi^{2}(v^{\varepsilon}_{\ell})+\theta(x^{\varepsilon}_{i},v^{\varepsilon}_{l})(\delta_{i+1}-\delta_{i})}{\delta_{i+1}\delta_{i}},\\\
\lambda^{\ell}_{i,i+1}=\frac{(1-\rho^{2})\varphi^{2}(v^{\varepsilon}_{\ell})+\theta(x^{\varepsilon}_{i},v^{\varepsilon}_{l})\delta_{i}}{\delta_{i+1}(\delta_{i+1}+\delta_{i})}.\end{array}\right.$
(13)

Note that $\lambda^{\ell}_{i,j}$ is well-defined for $i,j\in\mathcal{N}$,
$\ell\in\mathcal{M}$, where

$\displaystyle\theta(x^{\varepsilon}_{i},v^{\varepsilon}_{\ell})$
$\displaystyle=\frac{\mu(s^{\varepsilon}_{i,\ell},v^{\varepsilon}_{\ell})}{\nu(s^{\varepsilon}_{i,\ell})}-\frac{\nu^{\prime}(s^{\varepsilon}_{i,\ell})}{2}\varphi^{2}(v^{\varepsilon}_{\ell})+\frac{\rho}{2}\big{(}\sigma(v^{\varepsilon}_{\ell})\varphi^{\prime}(v^{\varepsilon}_{\ell})-\sigma^{\prime}(v^{\varepsilon}_{\ell})\varphi(v^{\varepsilon}_{\ell})\big{)}$
$\displaystyle-\rho\left((v^{\varepsilon}_{\ell}-V^{\varepsilon}_{0})\widehat{R}+\frac{\varepsilon^{H-\frac{1}{2}}b(v^{\varepsilon}_{\ell})}{\Gamma(H+1/2)}\right)\frac{\varphi(v^{\varepsilon}_{\ell})\Gamma(H+1/2)}{\varepsilon^{H-\frac{1}{2}}\sigma(v^{\varepsilon}_{\ell})},$

with $s^{\varepsilon}_{i,\ell}=g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{\ell}))$.

According to Song et al. (2016),
$(\widetilde{X}^{\varepsilon}_{t},\widetilde{V}^{\varepsilon}_{t})$ can be
represented as a one-dimensional CTMC with a $NM\times NM$ transition rate
matrix:

$\displaystyle{\bf\Lambda}$
$\displaystyle=\left(\begin{array}[]{cccc}q_{11}{\bf
I}_{N}+{\bf\Lambda}_{1}&q_{12}{\bf I}_{N}&\cdots&q_{1M}{\bf I}_{N}\\\
q_{21}{\bf I}_{N}&q_{22}{\bf I}_{N}+{\bf\Lambda}_{2}&\cdots&q_{2M}{\bf
I}_{N}\\\ \vdots&\vdots&\ddots&\vdots\\\ q_{M1}{\bf I}_{N}&q_{M2}{\bf
I}_{N}&\cdots&q_{MM}{\bf I}_{N}+{\bf\Lambda}_{M}\end{array}\right),$ (14)

where ${\bf I}_{N}$ is the $N\times N$ identity matrix, $q_{i,j}$ and
${\bf\Lambda}_{l}=(\lambda^{l}_{i,j})_{N\times N}$ are respectively given by
(11) and (13). Recall (8) and denote

$\widetilde{S}^{\varepsilon}_{t}:=g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{t}+\rho
f(\widetilde{V}^{\varepsilon}_{t})\big{)}.$

Then for any continuous function $\phi$, with the setting that
$g(s^{\varepsilon})-\rho
f(v^{\varepsilon})=x^{\varepsilon}=x^{\varepsilon}_{i}$,
$v^{\varepsilon}=v^{\varepsilon}_{l}$,

$\displaystyle\mathbb{E}\left[\phi(S^{\varepsilon}_{T},V^{\varepsilon}_{T})|S^{\varepsilon}_{t}=s^{\varepsilon},V^{\varepsilon}_{t}=v^{\varepsilon}\right]$
(15) $\displaystyle=$
$\displaystyle\mathbb{E}\left[\phi(g^{-1}\big{(}X^{\varepsilon}_{t}+\rho
f(V^{\varepsilon}_{t})\big{)},V^{\varepsilon}_{T})|X^{\varepsilon}_{t}=x^{\varepsilon},V^{\varepsilon}_{t}=v^{\varepsilon}\right]$
$\displaystyle\approx$
$\displaystyle\mathbb{E}\left[\phi(g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{t}+\rho
f(\widetilde{V}^{\varepsilon}_{t})\big{)},\widetilde{V}^{\varepsilon}_{T})|\widetilde{X}^{\varepsilon}_{t}=x^{\varepsilon}_{i},\widetilde{V}^{\varepsilon}_{t}=v^{\varepsilon}_{l}\right]$
$\displaystyle=$
$\displaystyle\mathbf{e}_{i,l}\cdot\exp(\mathbf{\Lambda}(T-t))\cdot\mathbf{\Phi},$

where ${\bf e}_{i,l}$ is a $1\times MN$ vector with all entries equal to $0$
except that the $(l-1)N+i$ entry is equal to $1$, ${\bf\Lambda}$ is given by
(14), and the payoff vector ${\bf\Phi}$ an $MN\times 1$ vectors with elements
${\bf\Phi}_{(l-1)N+i}=\phi(g^{-1}\big{(}x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l})\big{)},v^{\varepsilon}_{l})$, for $l\in\mathcal{M}$,
$i\in\mathcal{N}$.

###  3.2 Convergence analysis for the CTMC approximation

Denote $\mathbb{S}:=[s_{\min},s_{\max}]$ for
$-\infty<s_{\min}<s_{\max}<\infty$. Let
$\mathbb{S}^{\varepsilon}:=[s^{\varepsilon}_{\min},s^{\varepsilon}_{\max}]$,
$\mathbb{V}^{\varepsilon}:=[v^{\varepsilon}_{\min},v^{\varepsilon}_{\max}]$,
$\mathbb{X}^{\varepsilon}:=[x^{\varepsilon}_{\min},x^{\varepsilon}_{\max}]$
and $\mathbb{T}:=[0,T]$ be the range of the asset process
$S^{\varepsilon}_{t}$, volatility process $V^{\varepsilon}_{t}$, auxiliary
process $X^{\varepsilon}_{t}$ and time $t$, where
$-\infty<s^{\varepsilon}_{\min}<s^{\varepsilon}_{\max}<\infty$,
$-\infty<v^{\varepsilon}_{\min}<v^{\varepsilon}_{\max}<\infty$,
$x^{\varepsilon}_{\min}=g(s^{\varepsilon}_{\min})+\rho
f(v^{\varepsilon}_{\min})$ and
$x^{\varepsilon}_{\max}=g(s^{\varepsilon}_{\max})+\rho
f(v^{\varepsilon}_{\max})$. The convergence analysis of the CTMC approximation
is based on the following lemma.

###### Lemma 3.1

Let $S_{t}\in\mathbb{S}$ for $t\in\mathbb{T}$ be a Feller process whose
infinitesimal generator is given by

$\mathcal{L}G(s):=\lim_{t\to 0}\frac{\mathbb{E}[G(S_{t})|S_{0}=s]-G(s)}{t}.$

Let $\widetilde{S}_{t}^{n}\in{\mathbb{S}}^{n}$ be the continuous-time Markov
chain with the generator ${\mathcal{L}}^{n}$,
${\mathbb{S}}^{n}\subseteq{\mathbb{S}}$ and when n tends to infinity,
${\mathbb{S}}^{n}={\mathbb{S}}$. Assume that for each
$G\in\mathcal{C}^{2}({\mathbb{S}})$,

$\lim_{n\to\infty}\max_{s\in{\mathbb{S}}^{n}}|\mathcal{L}G(s)-{\mathcal{L}}^{n}G(s)|=0.$
(16)

Then $\widetilde{S}_{t}^{n}$ converges weakly to $S_{t}$ as $n$ goes to
infinity.

Proof According to Section 8.7 in Durrett (2018) and Theorem 10.1.1 in Kushner
and Dupuis (2001), it can be deduced from condition (16) that
$\widetilde{S}_{t}^{n}$ is tight. Then by Theorem 4.2.11 in Ethier (2009), we
obtain the weak convergence of $\widetilde{S}_{t}^{n}$ to $S_{t}$. This
completes the proof. $\quad\Box$

The first theorem gives the weak convergence of
$\widetilde{V}^{\varepsilon}_{t}$ to $V^{\varepsilon}_{t}$.

###### Theorem 3.1

Assuming that the grid interval $h$ satisfies: $h_{\ell}=\mathcal{O}(M^{-1})$,
$|h_{\ell}-h_{\ell+1}|\leq\mathcal{O}(M^{-2})$ and Assumption 2.1 holds, the
continuous-time Markov chain $\widetilde{V}^{\varepsilon}_{t}$ with the
generator ${\bm{Q}}=(q_{i,j})_{M\times M}$ converges weakly to
${V}^{\varepsilon}_{t}$, as $M$ goes to infinity.

Proof Theorem 2.4 shows the Markov property of $V_{t}^{\varepsilon}$.
According to the continuity of Riemann integral, $V_{t}^{\varepsilon}$ in (3)
has continuous simple path, therefore is a Feller process. Here we assume that
$\hbox{Prob}\left\\{V_{t}^{\varepsilon,\gamma}=v^{\varepsilon,\gamma}_{\ell}\big{|}V_{t}^{\varepsilon}=v^{\varepsilon}_{\ell}\right\\}=1.$
Moreover, for $G\in\mathcal{C}^{2}({\mathbb{V}}^{\varepsilon})$,
$\ell\in\mathcal{M}^{\circ}$,

$\displaystyle|\mathcal{L}G(v^{\varepsilon}_{\ell})-{\bm{Q}}G(v^{\varepsilon}_{\ell})|=$
$\displaystyle\Big{|}\frac{1}{2}(K^{\varepsilon})^{2}\sigma^{2}(v^{\varepsilon}_{\ell})G^{\prime\prime}(v^{\varepsilon}_{l})+\Big{(}(v^{\varepsilon}_{\ell}-V^{\varepsilon}_{0})\widehat{R}+K^{\varepsilon}b(v^{\varepsilon}_{\ell})\Big{)}G^{\prime}(v^{\varepsilon}_{\ell})$
$\displaystyle-\left(q_{l,l-1}G(v^{\varepsilon}_{l-1})+q_{l,l}G(v^{\varepsilon}_{\ell})+q_{l,l+1}G(v^{\varepsilon}_{l+1})\right)\Big{|}.$

Then we use (11), the triangle inequality, the Lipschitz continuity and the
linear growth of $b$ and $\sigma$ to give

$\displaystyle|\mathcal{L}G(v^{\varepsilon}_{\ell})-{\bm{Q}}G(v^{\varepsilon}_{\ell})|$
(17) $\displaystyle\leq$ $\displaystyle
C\Bigg{[}\left|G^{\prime\prime}(v^{\varepsilon}_{\ell})-\frac{2}{(h_{l+1}+h_{l})}\left(\frac{G(v^{\varepsilon}_{l+1})-G(v^{\varepsilon}_{\ell})}{h_{l+1}}-\frac{G(v^{\varepsilon}_{l})-G(v^{\varepsilon}_{l-1})}{h_{l}}\right)\right|$
$\displaystyle+\bigg{|}G^{\prime}(v^{\varepsilon}_{\ell})-\frac{1}{(h_{l+1}+h_{l})}\left(\frac{h_{l}(G(v^{\varepsilon}_{l+1})-G(v^{\varepsilon}_{\ell}))}{h_{l+1}}+\frac{h_{l+1}(G(v^{\varepsilon}_{l})-G(v^{\varepsilon}_{l-1}))}{h_{l}}\right)\bigg{|}\Bigg{]}.$

By Taylor expansion, we have that

$\displaystyle G(v^{\varepsilon}_{l+1})$
$\displaystyle=G(v^{\varepsilon}_{\ell})+G^{\prime}(v^{\varepsilon}_{\ell})h_{l+1}+G^{\prime\prime}(v^{\varepsilon}_{l})\frac{h^{2}_{l+1}}{2}+o(h^{2}_{l+1}),$
(18) $\displaystyle G(v^{\varepsilon}_{l-1})$
$\displaystyle=G(v^{\varepsilon}_{\ell})-G^{\prime}(v^{\varepsilon}_{\ell})h_{l}+G^{\prime\prime}(v^{\varepsilon}_{\ell})\frac{h^{2}_{l}}{2}+o(h_{l}^{2}).$

Plugging (18) into (LABEL:ineq:L-Q) and using the condition
$h_{\ell}=\mathcal{O}(M^{-1})$,
$|h_{\ell}-h_{\ell+1}|\leq\mathcal{O}(M^{-2})$, we obtain

$|\mathcal{L}G(v^{\varepsilon}_{\ell})-{\bm{Q}}^{\gamma}G(v^{\varepsilon}_{\ell})|\leq
o(1).$

Then by Lemma 3.1, $\widetilde{V}^{\varepsilon}_{t}$ converges weakly to
${V}^{\varepsilon}_{t}$ as $M$ goes to infinity. This completes the proof.
$\quad\Box$

The next theorem shows the weak convergence of
$\widetilde{X}^{\varepsilon}_{t}$ to $X^{\varepsilon}_{t}$.

###### Theorem 3.2

Assuming that Assumption 2.1 holds, the grid interval $h$, $\delta$ satisfy:
$h_{\ell}=\mathcal{O}(M^{-1})$,
$|h_{\ell}-h_{\ell+1}|\leq\mathcal{O}(M^{-2})$, and
$\delta_{i}=\mathcal{O}(N^{-1})$,
$|\delta_{i}-\delta_{i+1}|\leq\mathcal{O}(N^{-2})$, the continuous-time Markov
chain $\widetilde{X}^{\varepsilon}_{t}$ converges weakly to
${X}^{\varepsilon}_{t}$, as $M$, $N$ go to infinity.

Proof It is easy to verify that under Assumption 2.1,
$\theta(x^{\varepsilon},v^{\varepsilon})$, which is the drift term of
${X}^{\varepsilon}_{t}$ and defined by (9), is continuous and bounded for
$(x^{\varepsilon},v^{\varepsilon})\in\mathbb{X}^{\varepsilon}\times\mathbb{V}^{\varepsilon}$.
Then by similar arguments as Theorem 3.1, we can obtain the desired
conclusion. This completes the proof. $\quad\Box$

##  4 Option pricing under RSLV model

In this section, we use the semimartingale and CTMC approximation techniques
to give the explicit approximate expression of option prices under the RSLV
model.

Vanilla option prices for $(S_{T},V_{T})$ under the RSLV model (1) is given
by:

$\mathbb{E}\left[e^{-rT}\phi(S_{T},V_{T})|S_{0}=s,V_{0}=v\right],$

where

$\phi(S_{T},V_{T})=\left\\{\begin{aligned}
&(S_{T}-D)^{+}~{}~{}~{}~{}for~{}a~{}call,\\\
&(D-S_{T})^{+}~{}~{}~{}~{}for~{}a~{}put,\end{aligned}\right.$

and $r$ is the risk-free interest rate, $D$ is the strike price, and $T$ is
the maturity. After the semimartingale approximation of $(S_{T},V_{T})$ by
$(S^{\varepsilon}_{T},V^{\varepsilon}_{T})$, and the CTMC approximation by
$(\widetilde{S}^{\varepsilon}_{T},\widetilde{V}^{\varepsilon}_{T})$, the
European options prices under the RSLV model have the following approximate
formula.

######  Algorithm 4.1 (European options)

Given that $s=s^{\varepsilon}$, $v=v^{\varepsilon}=v^{\varepsilon}_{l}$ and
$x^{\varepsilon}_{i}=g(s^{\varepsilon})+\rho f(v^{\varepsilon}_{l})$, the
European option price under the RSLV model can be approximately calculated by

$\displaystyle\mathbb{E}\left[e^{-rT}\phi(S_{T},V_{T})|S_{0}=s,V_{0}=v\right]$
$\displaystyle\approx$
$\displaystyle\mathbb{E}\left[e^{-rT}\phi(S^{\varepsilon}_{T},V^{\varepsilon}_{T})|S^{\varepsilon}_{0}=s^{\varepsilon},V^{\varepsilon}_{0}=v^{\varepsilon}\right]$
$\displaystyle\approx$
$\displaystyle\mathbb{E}\left[e^{-rT}\phi(g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{t}+\rho
f(\widetilde{V}^{\varepsilon}_{t})\big{)},\widetilde{V}^{\varepsilon}_{T})|\widetilde{X}^{\varepsilon}_{0}=x^{\varepsilon}_{i},\widetilde{V}^{\varepsilon}_{0}=v^{\varepsilon}_{l}\right]$
$\displaystyle=$ $\displaystyle e^{-rT}\cdot{\bf
e}_{i,l}\cdot\exp({\bf\Lambda}T)\cdot{\bf\Phi^{(1)}}.$

Here ${\bf e}_{i,l}$ is a $1\times NM$ vector with all entries equal to $0$
except that the $(l-1)N+i$ entry is equal to $1$, ${\bf\Lambda}$ is given by
(14), and ${\bf\Phi^{(1)}}$ is an $NM\times 1$ vector with elements for
$l\in\mathcal{M}$, $i\in\mathcal{N}$,

${\bf\Phi^{(1)}}_{(l-1)N+i}=\phi\big{(}g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\big{)}=\left\\{\begin{aligned}
&\left(g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))-D\right)^{+}~{}~{}~{}~{}for~{}a~{}call,\\\
&\left(D-g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\right)^{+}~{}~{}~{}~{}for~{}a~{}put.\end{aligned}\right.$

Similarly, we have the semimartingale and CTMC approximate formula of barrier
options prices under the RSLV model.

######  Algorithm 4.2 (barrier options)

Given that $s=s^{\varepsilon}$, $v=v^{\varepsilon}=v^{\varepsilon}_{l}$ and
$x^{\varepsilon}_{i}=g(s^{\varepsilon})+\rho f(v^{\varepsilon}_{l})$, the
barrier option price under the RSLV model with $0\leq L<U<\infty$, can be
approximately calculated by

$\displaystyle\mathbb{E}\left[e^{-rT}\phi(S_{T},V_{T})\mathds{1}_{\\{L<S_{T}<U\\}}|S_{0}=s,V_{0}=v\right]$
$\displaystyle\approx$
$\displaystyle\mathbb{E}\left[e^{-rT}\phi(S^{\varepsilon}_{T},V^{\varepsilon}_{T})\mathds{1}_{\\{L<S_{T}<U\\}}|S^{\varepsilon}_{0}=s^{\varepsilon},V^{\varepsilon}_{0}=v^{\varepsilon}\right]$
$\displaystyle\approx$
$\displaystyle\mathbb{E}\left[e^{-rT}\phi(g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{T}+\rho
f(\widetilde{V}^{\varepsilon}_{T})\big{)},\widetilde{V}^{\varepsilon}_{T})\mathds{1}_{\\{L<g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{T}+\rho
f(\widetilde{V}^{\varepsilon}_{T})\big{)}<U\\}}|\widetilde{X}^{\varepsilon}_{0}=x^{\varepsilon}_{i},\widetilde{V}^{\varepsilon}_{0}=v^{\varepsilon}_{l}\right]$
$\displaystyle=$ $\displaystyle e^{-rT}\cdot{\bf
e}_{i,l}\cdot\exp({\bf\Lambda}T)\cdot{\bf\Phi^{(2)}}.$

Here ${\bf e}_{i,l}$ is a $1\times NM$ vector with all entries equal to $0$
except that the $(l-1)N+i$ entry is equal to $1$, ${\bf\Lambda}$ is given by
(14), and ${\bf\Phi^{(2)}}$ is an $NM\times 1$ vector with elements for
$l\in\mathcal{M}$, $i\in\mathcal{N}$,

${\bf\Phi^{(2)}}_{(l-1)N+i}=\left\\{\begin{aligned}
&\left(g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))-D\right)^{+}\mathds{1}_{\\{L<g^{-1}\big{(}x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l})\big{)}<U\\}}~{}~{}~{}~{}for~{}a~{}call,\\\
&\left(D-g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\right)^{+}\mathds{1}_{\\{L<g^{-1}\big{(}x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l})\big{)}<U\\}}~{}~{}~{}~{}for~{}a~{}put.\end{aligned}\right.$

Next we consider (finite-maturity) American options, whose prices are given by

$\max_{\tau\in\mathcal{T}}\mathbb{E}\left[e^{-r\tau}\phi(S_{\tau},V_{\tau})|S_{0}=s,V_{0}=v\right],$

where the set $\mathcal{T}$ comprises of the collection of
$\mathbb{F}$-stopping times taking values between $0$ and $T$. This means that
American options can be exercised at any time in $[0,T]$. $\mathcal{T}$ can be
approximated by a finite set of admissible exercise times
$\mathcal{T}^{n}:=\\{\tau_{i}\\}^{n}_{i=0}$, where $\tau_{i}=\frac{iT}{n}$,
where $n$ is the number of monitoring dates. Under admissible exercise times
set $\mathcal{T}^{n}$, the option is called the Bermudan option. By
semimartingale and CTMC approximations, the value of the American option under
the RSLV model can be approximately expressed as:

######  Algorithm 4.3 (American options)

Given that $s=s^{\varepsilon}$, $v=v^{\varepsilon}=v^{\varepsilon}_{l}$ and
$x^{\varepsilon}_{i}=g(s^{\varepsilon})+\rho f(v^{\varepsilon}_{l})$, the
Bermudan option price under RSLV model can be approximately calculated by

$\displaystyle\max_{\tau\in\mathcal{T}}\mathbb{E}\left[e^{-r\tau}\phi(S_{\tau},V_{\tau})|S_{0}=s,V_{0}=v\right]$
$\displaystyle\approx$
$\displaystyle\max_{\tau\in\mathcal{T}^{n}}\mathbb{E}\left[e^{-r\tau}\phi(S_{\tau},V_{\tau})|S_{0}=s,V_{0}=v\right]$
$\displaystyle\approx$
$\displaystyle\max_{\tau\in\mathcal{T}^{n}}\mathbb{E}\left[e^{-r\tau}\phi(S^{\varepsilon}_{\tau},V^{\varepsilon}_{\tau})|S^{\varepsilon}_{0}=s^{\varepsilon},V^{\varepsilon}_{0}=v^{\varepsilon}\right]$
$\displaystyle\approx$
$\displaystyle\max_{\tau\in\mathcal{T}^{n}}\mathbb{E}\left[e^{-r\tau}\phi(g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{\tau}+\rho
f(\widetilde{V}^{\varepsilon}_{\tau})\big{)},\widetilde{V}^{\varepsilon}_{\tau})|\widetilde{X}^{\varepsilon}_{0}=x^{\varepsilon}_{i},\widetilde{V}^{\varepsilon}_{0}=v^{\varepsilon}_{l}\right]=:B_{0},$

where $B_{0}={\bf e}_{i,l}\cdot{\bf B_{0}}$, and

$\left\\{\begin{aligned} {\bf B_{n}}&={\bf\Phi^{(1)}},\\\ {\bf
B_{i}}&=\max\left\\{e^{-rT/n}\exp({\bf\Lambda}T/n)\cdot{\bf
B_{i+1}},{\bf\Phi^{(1)}}\right\\},~{}i=n-1,n-2,\dots,0.\end{aligned}\right.$

Here $\bf{B_{i}}$ is a $NM\times 1$ vector, ${\bf e}_{i,l}$ a $1\times NM$
vector with all entries equal to $0$ except that the $(l-1)N+i$ entry is equal
to $1$, ${\bf\Lambda}$ is given by (14), and ${\bf\Phi^{(1)}}$ is an $NM\times
1$ vector with elements for $l\in\mathcal{M}$, $i\in\mathcal{N}$,

${\bf\Phi^{(1)}}_{(l-1)N+i}=\phi\big{(}g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\big{)}=\left\\{\begin{aligned}
&\left(g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))-D\right)^{+}~{}~{}~{}~{}for~{}a~{}call,\\\
&\left(D-g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\right)^{+}~{}~{}~{}~{}for~{}a~{}put.\end{aligned}\right.$

######  Remark 4.1 (Convergence analysis of options pricing)

Note that Theorem 2.2 and Theorem 3.2 show respectively the weak convergence
of the semimartingale approximation and the CTMC approximation. For European
and barrier options, since the functions $\phi$, $g$, $f$ are all continuous,
the options prices obtained by the semimartingale and CTMC approximations
converge to the original RSLV options prices, due to the continuous mapping
theorem (see, e.g., Kushner and Dupuis (2001)). For American options, the weak
convergence follows from the following three types of convergence: the
convergence of time step discretization (refer to Chapter $2.4$ in Quecke
(2007)), semimartingale approximation (due to continuous mapping theorem), and
the CTMC approximation (see Section 6 in Eriksson and Pistorius (2015)).

##  5 Numerical results

In this section, we extend the traditional stochastic local volatility model
(see Table 1) to rough stochastic local volatility model (see Table 2). Then a
series of examples for European, barrier and American options are used to
illustrate the accuracy and efficiency of the semimartingale and CTMC
approximation method introduced in Section 4. All numerical experiments are
carried out with Matlab R$2016$a on a Core i$7$ desktop with $16$GB RAM and
speed $3.60$ GHz.

Table 1 collects some popular stochastic local volatility models in the
literature. The traditional stochastic local volatility models are all driven
by standard Brownian motions. However, as confirmed by Gatheral et al. (2018),
for a very wide range of assets, historical fluctuation time series exhibit
rougher behavior than the Brownian motion. Therefore, it is very practical to
recast the traditional stochastic local volatility models in the rough
setting. Similar to the way El Euch and Rosenbaum (2019) deducing the rough
Heston model, we generalize all models listed in Table 1 to their rough
counterparts. The results are listed in Table 2. In order to numerically solve
the option pricing problem under RSLV models, we use the semimartingale
approximation introduced in Section 2.2 and auxiliary processes introduced in
Lemma 3.1 to approximately transform those models into PSLV models with two
independent Brownian motions. The results are listed in Table 3.

Table 1: Examples of the stochastic local volatility models.

Heston | $dS_{t}=(r-q)S_{t}dt+S_{t}\sqrt{V_{t}}dW_{t}$ | $\eta,\vartheta\in{\mathbb{R}}$
---|---|---
(Heston (1993)) | $dV_{t}=\eta(\vartheta-V_{t})dt+\sigma\sqrt{V_{t}}dB_{t}$ | $\sigma>0$
$4/2$ model | $dS_{t}=(r-q)S_{t}dt+S_{t}[a\sqrt{V_{t}}+b/\sqrt{V_{t}}]dW_{t}$ | $\eta,\vartheta,a,b\in{\mathbb{R}}$
(Grasselli (2017)) | $dV_{t}=\eta(\vartheta-V_{t})dt+\sigma\sqrt{V_{t}}dB_{t}$ | $\sigma>0$
$\alpha$-Hyper | $dS_{t}=(r-q)S_{t}dt+S_{t}\exp(V_{t})dW_{t}$ | $\eta\in{\mathbb{R}}$
(Da Fonseca and Martini (2016)) | $dV_{t}=(\eta-\vartheta\exp(aV_{t}))dt+\sigma dB_{t}$ | $\vartheta,a,\sigma>0$
SABR | $dS_{t}=V_{t}S_{t}^{\beta}dW_{t}$ | $\beta\in[0,1)$
(Hagan et al. (2002)) | $dV_{t}=\sigma V_{t}dB_{t}$ | $\sigma>0$
Heston-SABR | $dS_{t}=(r-q)S_{t}dt+\sqrt{V_{t}}S_{t}^{\beta}dW_{t}$ | $\beta\in[0,1)$
(Van der Stoep et al. (2014)) | $dV_{t}=\eta(\vartheta-V_{t})dt+\sigma\sqrt{V_{t}}dB_{t}$ | $\eta,\vartheta,\sigma>0$
Quadratic SLV | $dS_{t}=(r-q)S_{t}dt+\sqrt{V_{t}}(aS_{t}^{2}+bS_{t}+c)dW_{t}$ | $4ac>b^{2}$
(Lipton (2002)) | $dV_{t}=\eta(\vartheta-V_{t})dt+\sigma\sqrt{V_{t}}dB_{t}$ | $a,\eta,\vartheta,\sigma>0$

  * *

Here $r$ is the risk-free interest rate, $q$ is the dividend yield, $V_{t}$ is
the volatility and satisfies different stochastic differential equation for
different models.

Table 2: Examples of the rough stochastic local volatility models.

Rough Heston | $dS_{t}=(r-q)S_{t}dt+S_{t}\sqrt{V_{t}}dW_{t}$ | $\eta,\vartheta\in{\mathbb{R}}$
---|---|---
$V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left(\eta(\vartheta-V_{s})ds+\sigma\sqrt{V_{s}}dB_{s}\right)$ | $\sigma>0$
Rough $4/2$ model | $dS_{t}=(r-q)S_{t}dt+S_{t}[a\sqrt{V_{t}}+b/\sqrt{V_{t}}]dW_{t}$ | $\eta,\vartheta,a,b\in{\mathbb{R}}$
$V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left(\eta(\vartheta-V_{s})ds+\sigma\sqrt{V_{s}}dB_{s}\right)$ | $\sigma>0$
Rough $\alpha$-Hyper | $dS_{t}=(r-q)S_{t}dt+S_{t}\exp(V_{t})dW_{t}$ | $\eta\in{\mathbb{R}}$
$V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left((\eta-\vartheta\exp(aV_{s}))ds+\sigma dB_{s}\right)$ | $\vartheta,a,\sigma>0$
Rough SABR | $dS_{t}=V_{t}S_{t}^{\beta}dW_{t}$ | $\beta\in[0,1)$
$V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left(\sigma V_{s}dB_{s}\right)$ | $\sigma>0$
Rough Heston-SABR | $dS_{t}=(r-q)S_{t}dt+\sqrt{V_{t}}S_{t}^{\beta}dW_{t}$ | $\beta\in[0,1)$
$V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left(\eta(\vartheta-V_{s})ds+\sigma\sqrt{V_{s}}dB_{s}\right)$ | $\eta,\vartheta,\sigma>0$
Rough quadratic SLV | $dS_{t}=(r-q)S_{t}dt+\sqrt{V_{t}}(aS_{t}^{2}+bS_{t}+c)dW_{t}$ | $4ac>b^{2}$
$V_{t}=V_{0}+\int^{t}_{0}K(t,s)\left(\eta(\vartheta-V_{s})ds+\sigma\sqrt{V_{s}}dB_{s}\right)$ | $a,\eta,\vartheta,\sigma>0$

  * *

Here $r$ is the risk-free interest rate, $q$ is the dividend yield, $V_{t}$ is
the volatility process with initial value $V_{0}>0$, and
$K(t,s)=(t-s)^{H-\frac{1}{2}}/\Gamma(H+1/2)$ is the fractional kernel with the
Hurst parameter $H\in(0,1/2)$.

Table 3: Examples of dynamics and variance transforms for the perturbed
stochastic local volatility models.

Perturbed Heston | $\left\\{\begin{array}[]{ll}dX^{\varepsilon}_{t}=\left(r-q-\frac{V^{\varepsilon}_{t}}{2}-\frac{\rho\eta(\vartheta-V^{\varepsilon}_{t})}{\sigma}+\frac{\rho\chi(V_{t}^{\varepsilon,\gamma})}{K^{\varepsilon}\sigma}\right)dt+\sqrt{1-\rho^{2}}\sqrt{V^{\varepsilon}_{t}}dB^{\bot}_{t}\\\ dV^{\varepsilon}_{t}=-\chi(V_{t}^{\varepsilon,\gamma})dt+K^{\varepsilon}\left(\eta(\vartheta-V^{\varepsilon}_{t})dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\right)\\\ dV_{t}^{\varepsilon,\gamma}=(-\gamma V_{t}^{\varepsilon,\gamma}+\eta(\vartheta-V^{\varepsilon}_{t}))dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\\\ X^{\varepsilon}_{t}=\log S^{\varepsilon}_{t}-\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma}\end{array}\right.$
---|---
Perturbed $4/2$ model | $\left\\{\begin{array}[]{ll}dX^{\varepsilon}_{t}=\left(r-q+\frac{\rho}{\sigma}(\frac{\chi(V_{t}^{\varepsilon,\gamma})}{K^{\varepsilon}}-\eta(\vartheta-V^{\varepsilon}_{t}))(a+\frac{b}{V^{\varepsilon}_{t}})-\frac{\rho K^{\varepsilon}b\sigma}{2V^{\varepsilon}_{t}}\right)dt\\\ ~{}~{}~{}~{}~{}-\frac{1}{2}(a\sqrt{V^{\varepsilon}_{t}}+\frac{b}{\sqrt{V^{\varepsilon}_{t}}})^{2}dt+\sqrt{1-\rho^{2}}(a\sqrt{V^{\varepsilon}_{t}}+\frac{b}{\sqrt{V^{\varepsilon}_{t}}})dB^{\bot}_{t}\\\ dV^{\varepsilon}_{t}=-\chi(V_{t}^{\varepsilon,\gamma})dt+K^{\varepsilon}\left(\eta(\vartheta-V^{\varepsilon}_{t})dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\right)\\\ dV_{t}^{\varepsilon,\gamma}=\left(-\gamma V_{t}^{\varepsilon,\gamma}+\eta(\vartheta-V^{\varepsilon}_{t})\right)dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\\\ X_{t}=\log S^{\varepsilon}_{t}-\rho\frac{aV^{\varepsilon}_{t}+b\log V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma}\end{array}\right.$
Perturbed $\alpha$-Hyper | $\left\\{\begin{array}[]{l}dX^{\varepsilon}_{t}=(r-q-\frac{\rho\exp(V^{\varepsilon}_{t})}{\sigma}(\frac{\chi(V_{t}^{\varepsilon,\gamma})}{K^{\varepsilon}}-\eta+\vartheta\exp(aV^{\varepsilon}_{t}))dt\\\ ~{}~{}~{}~{}~{}+(\frac{\rho K^{\varepsilon}\sigma\exp(V^{\varepsilon}_{t})}{2}-\frac{\exp(2V^{\varepsilon}_{t})}{2})dt+\sqrt{1-\rho^{2}}\exp(V^{\varepsilon}_{t})dB^{\bot}_{t}\\\ dV^{\varepsilon}_{t}=-\chi(V_{t}^{\varepsilon,\gamma})dt+K^{\varepsilon}\left((\eta-\vartheta\exp(aV^{\varepsilon}_{t}))dt+\sigma dB_{t}\right)\\\ dV_{t}^{\varepsilon,\gamma}=\left(-\gamma V_{t}^{\varepsilon,\gamma}+\eta-\vartheta\exp(aV^{\varepsilon}_{t})\right)dt+\sigma dB_{t}\\\ X^{\varepsilon}_{t}=\log S^{\varepsilon}_{t}-\frac{\rho\exp(V^{\varepsilon}_{t})}{K^{\varepsilon}\sigma}\end{array}\right.$
Perturbed SABR | $\left\\{\begin{array}[]{ll}dX^{\varepsilon}_{t}=\left(-\frac{\beta v^{2}_{t}}{2(1-\beta)\left(X^{\varepsilon}_{t}+\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma}\right)}+\frac{\rho\chi(V_{t}^{\varepsilon,\gamma})}{K^{\varepsilon}\sigma}\right)dt+\sqrt{1-\rho^{2}}V^{\varepsilon}_{t}dB^{\bot}_{t}\\\ dV^{\varepsilon}_{t}=-\chi(V_{t}^{\varepsilon,\gamma})dt+K^{\varepsilon}\sigma V^{\varepsilon}_{t}dB_{t}\\\ dV_{t}^{\varepsilon,\gamma}=-\gamma V_{t}^{\varepsilon,\gamma}dt+\sigma V^{\varepsilon}_{t}dB_{t}\\\ X^{\varepsilon}_{t}=\frac{(S^{\varepsilon}_{t})^{1-\beta}}{1-\beta}-\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma}\end{array}\right.$
Perturbed Heston-SABR | $\left\\{\begin{array}[]{ll}dX^{\varepsilon}_{t}=\left((r-q)(1-\beta)(X^{\varepsilon}_{t}+\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma})-\frac{\beta V^{\varepsilon}_{t}}{2(1-\beta)(X^{\varepsilon}_{t}+\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma})}\right)dt\\\ ~{}~{}~{}~{}~{}+\frac{\rho}{\sigma}\left(\frac{\chi(V_{t}^{\varepsilon,\gamma})}{K^{\varepsilon}}-\eta(\vartheta-V^{\varepsilon}_{t})\right)dt+\sqrt{1-\rho^{2}}\sqrt{v_{t}}dB^{\bot}_{t}\\\ dV^{\varepsilon}_{t}=-\chi(V_{t}^{\varepsilon,\gamma})dt+K^{\varepsilon}\left(\eta(\vartheta-V^{\varepsilon}_{t})dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\right)\\\ dV_{t}^{\varepsilon,\gamma}=(-\gamma V_{t}^{\varepsilon,\gamma}+\eta(\vartheta-V^{\varepsilon}_{t}))dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\\\ X^{\varepsilon}_{t}=\frac{(S^{\varepsilon}_{t})^{1-\beta}}{1-\beta}-\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma}\end{array}\right.$
Perturbed quadratic SLV | $\left\\{\begin{array}[]{ll}dX^{\varepsilon}_{t}=\left(\frac{(r-q)S^{\varepsilon}_{t}}{a(S^{\varepsilon}_{t})^{2}+bS^{\varepsilon}_{t}+c}-\frac{2aS^{\varepsilon}_{t}+b}{2}V^{\varepsilon}_{t}+\frac{\rho}{\sigma}(\frac{\chi(V_{t}^{\varepsilon,\gamma})}{K^{\varepsilon}}-\eta(\vartheta-V^{\varepsilon}_{t}))\right)dt\\\ ~{}~{}~{}~{}~{}+\sqrt{1-\rho^{2}}\sqrt{v_{t}}dB^{\bot}_{t}\\\ dV^{\varepsilon}_{t}=-\chi(V_{t}^{\varepsilon,\gamma})dt+K^{\varepsilon}\left(\eta(\vartheta-V^{\varepsilon}_{t})dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\right)\\\ dV_{t}^{\varepsilon,\gamma}=(-\gamma V_{t}^{\varepsilon,\gamma}+\eta(\vartheta-V^{\varepsilon}_{t}))dt+\sigma\sqrt{V^{\varepsilon}_{t}}dB_{t}\\\ X^{\varepsilon}_{t}=\frac{2\arctan\left(\frac{2aS^{\varepsilon}_{t}+b}{\sqrt{4ac-b^{2}}}\right)}{\sqrt{4ac-b^{2}}}-\frac{\rho V^{\varepsilon}_{t}}{K^{\varepsilon}\sigma}\end{array}\right.$

  * *

Here $r$ is the risk-free interest rate, $q$ the dividend yield,
$X^{\varepsilon}_{t}$, $V^{\varepsilon}_{t}$, $V_{t}^{\varepsilon,\gamma}$ and
$S^{\varepsilon}_{t}$ defined by (8),
$\chi(V_{t}^{\varepsilon,\gamma}):=\int^{\infty}_{0}e^{-\gamma\varepsilon}\gamma
V_{t}^{\varepsilon,\gamma}m(d\gamma)$,
$K^{\varepsilon}=\varepsilon^{H-\frac{1}{2}}/\Gamma(H+1/2)$ the with the
perturbation parameter $0<\varepsilon<<1$ and the Hurst parameter
$H\in(0,1/2)$.

Next, we use the method introduced in Section 4 to numerically solve the
European and barrier options pricing under RSLV models. All models considered
in this paper share the same set of parameters:

$\displaystyle
S_{0}=10,~{}V_{0}=0.04,~{}T=1,~{}D=4,~{}\rho=-0.75,~{}r=0,~{}\sigma=0.8,~{}\beta=0.7,$
$\displaystyle\eta=4,~{}\vartheta=0.035,~{}a=0.02,~{}b=0.05,~{}c=1,~{}L=2,~{}U=15,~{}H=0.12.$

In the numerical implementation, we use the same boundary points
$v_{\min}=10^{-3}V_{0}$, $v_{\max}=4V_{0}$, $s_{\min}=10^{-3}S_{0}$,
$s_{\max}=4S_{0}$ and $x_{\min}=10^{-3}X_{0}$, $x_{\max}=4X_{0}$, with
$X_{0}=g(S_{0})-\rho f(v_{0})$ and choose the piecewise uniform grids.

######  Example 5.1 (European and barrier options)

The value of European and barrier options can be calculated by Algorithm 4.1
and Algorithm 4.2. Due to the independence between the auxiliary process
$X^{\varepsilon}_{t}$ and the volatility process $V^{\varepsilon}_{t}$, the
following fast algorithm can be used to speed up calculations.

######  Algorithm 5.1 (A fast 2 dimensional CTMC algorithm for European and
barrier options)

Under the setting that $g(s^{\varepsilon})-\rho
f(v^{\varepsilon})=x^{\varepsilon}=x^{\varepsilon}_{i}$,
$v^{\varepsilon}=v^{\varepsilon}_{\ell}$,

$\displaystyle\mathbb{E}\left[\phi(S^{\varepsilon}_{T},V^{\varepsilon}_{T})|S^{\varepsilon}_{t}=s^{\varepsilon},V^{\varepsilon}_{t}=v^{\varepsilon}\right]$
(19) $\displaystyle\approx$
$\displaystyle\mathbb{E}\left[\phi(g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{t}+\rho
f(\widetilde{V}^{\varepsilon}_{t})\big{)},\widetilde{V}^{\varepsilon}_{T})|\widetilde{X}^{\varepsilon}_{t}=x^{\varepsilon}_{i},\widetilde{V}^{\varepsilon}_{t}=v^{\varepsilon}_{\ell}\right]$
$\displaystyle=$
$\displaystyle\sum^{M}_{j=1}(e^{{\bm{Q}}T})_{\ell,j}\mathbb{E}\left[\phi(g^{-1}\big{(}\widetilde{X}^{\varepsilon}_{T}+\rho
f(\widetilde{V}^{\varepsilon}_{T})\big{)},v^{\varepsilon}_{j})|\widetilde{X}^{\varepsilon}_{t}=x^{\varepsilon}_{i}\right]$
$\displaystyle=$
$\displaystyle\sum^{M}_{j=1}(\exp({\bm{Q}}(T-t)))_{\ell,j}\times\mathbf{e}_{i}\cdot\exp(\mathbf{\Lambda}_{l}(T-t))\cdot\mathbf{\Phi}_{\ell},$

where ${\bf e}_{i}$ is a $1\times N$ vector with all entries equal to $0$
except that the $i$th entry is equal to $1$, ${\bm{Q}}$ is given by (11),
${\bf\Lambda}_{\ell}$ is given by (13), and the payoff vector ${\bf\Phi}_{l}$
is a $N\times 1$ vectors with elements
$({\bf\Phi}_{\ell})_{i}=\phi(g^{-1}\big{(}x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{\ell})\big{)},v^{\varepsilon}_{\ell})$, for
$\ell\in\mathcal{M}$, $i\in\mathcal{N}$. For European option, we have

$\phi\big{(}g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\big{)}=\left\\{\begin{aligned}
&\left(g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))-D\right)^{+}~{}~{}~{}~{}for~{}a~{}call,\\\
&\left(D-g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\right)^{+}~{}~{}~{}~{}for~{}a~{}put.\end{aligned}\right.$

and for barrier option, there is

$\phi\big{(}g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\big{)}=\left\\{\begin{aligned}
&\left(g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))-D\right)^{+}\mathds{1}_{\\{L<g^{-1}\big{(}x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l})\big{)}<U\\}}~{}~{}~{}~{}for~{}a~{}call,\\\
&\left(D-g^{-1}(x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l}))\right)^{+}\mathds{1}_{\\{L<g^{-1}\big{(}x^{\varepsilon}_{i}+\rho
f(v^{\varepsilon}_{l})\big{)}<U\\}}~{}~{}~{}~{}for~{}a~{}put.\end{aligned}\right.$

###### Remark 5.1

When applying the CTMC method to a two-dimensional problem, the two-
dimensional probability transition matrix is usually converted into a large
matrix (see Song et al. (2016) and Xi et al. (2019)). In this way, we need to
calculate the matrix exponential where the exponent is a $NM\times NM$ matrix.
This is time consuming, and requires a lot of computer storage space for
intermediate inputs. Since the two stochastic processes are separated by Lemma
3.1, we only need to compute $N\times N$ matrices $M$ times using our method.
This greatly reduces the unnecessary calculation, reduces the computer storage
space, and improves the operation efficiency. See Figure 1 for the comparison
of actual operation efficiency.

Figure 1: The plot of the number of grid points versus logarithmic of CPU
times (seconds) for European (1) and barrier (1) options under rough Heston
model. Here we set $M=N$. The blue lines in 1 and 1 are obtained by Algorithm
4.1 and 4.2, respectively. The red lines are obtained by Algorithm 5.1.  Table
4: European call option under the RSLV models.

$\varepsilon$ | Relative Errors |
---|---|---
R-H |  R-$4/2$ |  R-$\alpha$-H | R-SABR | R-H-SABR | R-Q-SLV
$10^{-4}$ | 4.71e-3 | 9.98e-3 | 5.76e-3 | 8.78e-3 | 4.21e-3 | 5.86e-3
$10^{-5}$ | 2.84e-3 | 9.68e-3 | 2.56e-3 | 3.81e-3 | 2.86e-3 | 4.53e-3
$10^{-6}$ | 2.07e-3 | 7.92e-3 | 1.22e-3 | 1.75e-3 | 2.32e-3 | 3.97e-3
$10^{-7}$ | 1.75e-3 | 7.55e-3 | 6.66e-4 | 8.84e-4 | 2.09e-3 | 3.44e-3
$10^{-8}$ | 1.61e-3 | 6.57e-3 | 4.34e-4 | 5.25e-4 | 1.99e-3 | 2.98e-3
benchmark | 6.0545 | 0.0362 | 6.0001 | 4.9269 | 6.0018 | 6.0000

  * *

Here “R-H”, “R-$4/2$”, “R-$\alpha$-H”, “R-SABR”, “R-H-SABR”, “R-Q-SLV”
represent “rough Heston”, “rough $4/2$ model”, “rough $\alpha$-Hyper”, “rough
SABR”, “rough Heston-SABR”, “rough quadratic SLV” models respectively.
Benchmarks are obtained by fast simulation algorithm based on the Monte Carlo
method in Ma and Wu (2021) with $10^{5}$ simple paths, and they take about
$610.86$ seconds on average. The results in the table are calculated with
$N=M=100$ via Algorithm 5.1, and take only $0.18$ seconds on average.

Figure 2: The logarithmic plot of the number of grid points versus relative
error for European options under rough stochastic local volatility models via
Algorithm 5.1. Figures 2, 2, 2, 2, 2, 2 correspond to rough Heston, $4/2$,
$\alpha$-Hyper, SABR, Heston-SABR, quadratic SLV models, respectively. The red
lines are obtained by fixing $M=90$, and then increasing the grid points of
$N$ from $50$ to $90$. The blue lines are obtained by fixing $N=90$ and
changing $M$. Here we take the prices obtained by $M=N=100$,
$\varepsilon=10^{-10}$ as the benchmarks to calculate relative errors.  Table
5: Barrier call option under the RSLV models.

$\varepsilon$ | Relative Errors |
---|---|---
R-H |  R-$4/2$ |  R-$\alpha$-H | R-SABR | R-H-SABR | R-Q-SLV
$10^{-4}$ | 4.21e-3 | 7.52e-3 | 3.54e-3 | 1.51e-3 | 2.84e-3 | 9.34e-3
$10^{-5}$ | 2.86e-3 | 6.18e-3 | 1.68e-3 | 7.42e-4 | 2.07e-3 | 5.76e-3
$10^{-6}$ | 2.32e-3 | 5.63e-3 | 9.07e-4 | 4.21e-4 | 1.75e-3 | 2.56e-3
$10^{-7}$ | 2.08e-3 | 5.40e-3 | 5.86e-4 | 2.88e-4 | 1.61e-3 | 1.22e-3
$10^{-8}$ | 1.99e-3 | 5.30e-3 | 4.53e-4 | 2.32e-4 | 1.56e-3 | 6.66e-4
benchmark | 6.0492 | 0.0345 | 5.9753 | 4.8099 | 6.0000 | 5.9814

  * *

Here ”R-H”, ”R-$4/2$”, ”R-$\alpha$-H”, ”R-SABR”, ”R-H-SABR”, ”R-Q-SLV”
represent ”rough Heston”, ”rough $4/2$ model”, ”rough $\alpha$-Hyper”, ”rough
SABR”, ”rough Heston-SABR”, ”rough rough quadratic SLV” models respectively.
Benchmarks are obtained by fast simulation algorithm based on Monte Carlo
method in Ma and Wu (2021) with $10^{5}$ simple paths, and they take about
$800.13$ seconds on average. The results in the table are calculated with
$N=M=100$ via Algorithm 5.1, and they take only $0.20$ seconds on average.

Figure 3: The logarithmic plot of the number of grid points versus relative
error for barrier options under rough stochastic local volatility models via
Algorithm 5.1. Figures 3, 3, 3, 3, 3, 3 correspond to rough Heston, $4/2$,
$\alpha$-Hyper, SABR, Heston-SABR, quadratic SLV models, respectively. The red
lines are obtained by fixing $M=90$, and then increasing the grid points of
$N$ from $50$ to $90$. The blue ones are obtained by fixing $N=90$ and
changing $M$. Here we take the prices obtained by $M=N=100$,
$\varepsilon=10^{-10}$ as the benchmarks to calculate relative errors.

Table 4 shows that with the decrease of $\varepsilon$, the European options
prices under RSLV models will converge to the benchmarks. For barrier options,
there are similar results, which are listed in Table 5. These results show the
accuracy of the semimartingale and CTMC approximation algorithm. From Figures
2 and 3, we see that for a fixed $\varepsilon$, increasing the number of CTMC
grids will make the relative error between the calculated results and the
benchmark gradually decrease. It is worth mentioning that, for fixed
$N=M=100$, the average CPU times to calculate the European and barrier option
prices are respectively $0.18$ and $0.20$ seconds. This shows the very high
efficiency of our algorithm.

Table 6: American call option under the RSLV models.

$\varepsilon$ | Relative Errors |
---|---|---
R-H | R-$4/2$ | R-$\alpha$-H | R-SABR | R-H-SABR | R-Q-SLV
$10^{-4}$ | 9.82e-3 | 9.99e-3 | 8.53e-3 | 9.25e-3 | 7.41e-3 | 8.88e-3
$10^{-5}$ | 6.73e-3 | 9.36e-3 | 6.25e-3 | 8.66e-3 | 5.55e-3 | 7.21e-3
$10^{-6}$ | 4.21e-3 | 8.95e-3 | 4.33e-3 | 7.99e-3 | 4.21e-3 | 6.85e-3
$10^{-7}$ | 3.75e-3 | 8.50e-3 | 3.28e-3 | 7.21e-3 | 3.48e-3 | 6.01e-3
$10^{-8}$ | 2.88e-3 | 8.12e-3 | 2.79e-3 | 6.54e-3 | 3.01e-3 | 5.75e-3
benchmark | 6.0635 | 0.0418 | 6.1111 | 6.0000 | 6.4410 | 7.1658

  * *

Here “R-H”, “R-$4/2$”, “R-$\alpha$-H”, “R-SABR”, “R-H-SABR”, “R-Q-SLV”
represent “rough Heston”, “rough $4/2$ model”, “rough $\alpha$-Hyper”, “rough
SABR”, “rough Heston-SABR”, “rough rough quadratic SLV” models respectively.
Benchmarks obtained by fast simulation algorithm based on least squares Monte
Carlo method with $10^{5}$ simple paths, and they take about $875.26$ seconds
on average. The results in the table are calculated with $N=M=100$ via
Algorithm 4.3, and they take $91.24$ seconds on average.

Figure 4: The logarithmic plot of the number of grid points versus relative
error for American options under rough stochastic local volatility models via
Algorithm 4.3. Figures 4, 4, 4, 4, 4, 4 correspond to rough Heston, $4/2$,
$\alpha$-Hyper, SABR, Heston-SABR, quadratic SLV models, respectively. The red
lines are obtained by fixing $n=M=90$, and then increasing the grid points of
$N$ from $50$ to $90$. The blue ones are obtained by fixing $n=N=90$ and
changing $M$, and the green ones are obtained by fixing $N=M=90$ and changing
$n$. Here we take the prices obtained by $n=N=M=100$, $\varepsilon=10^{-10}$
as the benchmarks to calculate relative errors.

######  Example 5.2 (American option)

In this example, we use Algorithm 4.3 to calculate the price of American
options under the six models listed in the Table 2.

Table 6 shows that with the decrease of $\varepsilon$, the American options
prices under RSLV models will converge to the benchmark. This demonstrates the
accuracy of the semimartingale and CTMC approximation algorithm. In addition
to showing the convergence results of CTMC grid points $N$ and $M$, Figure 4
also shows the influence of the change in grid points on the relative error in
the time direction. These results show the convergence of Algorithm 4.3. Since
the value of the American option depends on the optimal stopping time, the
fast Algorithm 5.1 is not available. However, the average CPU time to
calculate American option prices is $91.24$ seconds, which is significantly
less than the least squares Monte Carlo method. Moreover, compared with Table
4 in Goudenège et al. (2020), which spend a significant amount of time to
price American options under the simple Bergomi model, Algorithm 4.3 is
clearly more effective and widely applicable for pricing American options
under general rough stochastic local volatility models.

##  6 Conclusions

In this paper, we first propose a perturbation approach to approximate the
rough stochastic local volatility model by a perturbed stochastic local
volatility model, which is a semimartingale. After further expressing it as a
Markovian form, we propose a continuous-time Markov chain approximation
approach to derive the semi-explicit expressions of European, barrier and
American options prices. The approximate expression obtained by this method is
proved to converge to the solution of the original option pricing problem
under the RSLV model. Numerical experiments demonstrate the accuracy and very
high efficiency of the method for several RSLV models, including rough Heston
and rough SABR models, etc.

## References

  * Abi Jaber and El Euch (2019a) Abi Jaber, E., Larsson, M. and Pulido, S. (2019). Affine volterra processes. The Annals of Applied Probability, 29, 3155–3200.
  * Abi Jaber and El Euch (2019b) Abi Jaber, E. and El Euch, O. (2019). Markovian structure of the Volterra Heston model. Statistics and Probability Letters, 149, 63–72.
  * Abi Jaber and El Euch (2019c) Abi Jaber, E. and El Euch, O. (2019). Multifactor approximation of rough volatility models. SIAM Journal on Financial Mathematics, 10, 309–349.
  * Alfonsi and Kebaier, (2021) Alfonsi, A. and Kebaier, A. (2021). Approximation of Stochastic Volterra Equations with kernels of completely monotone type. arXiv preprint arXiv:2102.13505.
  * Bayer et al. (2020) Bayer, C., Ben Hammouda, C. and Tempone, R. (2020). Hierarchical adaptive sparse grids and quasi-Monte Carlo for option pricing under the rough Bergomi model. Quantitative Finance, 20, 1457–1473.
  * Bayer and Breneis (2021) Bayer, C. and Breneis, S. (2021). Markovian approximations of stochastic Volterra equations with the fractional kernel. arXiv preprint arXiv:2108.05048.
  * Bayer et al. (2016) Bayer, C., Friz, P. and Gatheral, J. (2016). Pricing under rough volatility. Quantitative Finance, 16, 887–904.
  * Bayer et al. (2020) Bayer, C., Friz, P.K., Gassiat, P., Martin, J. and Stemper, B. (2020). A regularity structure for rough volatility. Mathematical Finance, 30, 782–832.
  * Cai et al. (2015) Cai, N., Song, Y. and Kou, S. (2015). A general framework for pricing Asian options under Markov processes. Operations Research, 63, 540–554.
  * Callegaro et al. (2021) Callegaro, G., Grasselli, M. and Pages, G. (2021). Fast hybrid schemes for fractional Riccati equations (rough is not so tough). Mathematics of Operations Research, 46, 221–254.
  * Carmona et al. (2000) Carmona, P., Coutin, L. and Montseny, G. (2000). Approximation of some Gaussian processes. Statistical Inference for Stochastic Processes, 3, 161–171.
  * Cohen and Elliott (2015) Cohen, S.N. and Elliott, R.J. (2015). Stochastic Calculus and Applications (Vol. 2). New York: Birkhäuser.
  * Cui et al. (2018) Cui, Z., Kirkby, J.L. and Nguyen, D. (2018). A general valuation framework for SABR and stochastic local volatility models. SIAM Journal on Financial Mathematics, 9, 520–563.
  * Cui et al. (2019) Cui, Z., Kirkby, J.L. and Nguyen, D. (2019). Continuous-time Markov chain and regime switching approximations with applications to options pricing. In Modeling, stochastic control, optimization, and applications (pp. 115-146). Springer, Cham.
  * Da Fonseca and Martini (2016) Da Fonseca, J. and Martini, C. (2016). The $\alpha$-hypergeometric stochastic volatility model. Stochastic Processes and their Applications, 126, 1472–1502.
  * Dandapani et al. (2021) Dandapani, A., Jusselin, P. and Rosenbaum, M. (2021). From quadratic Hawkes processes to super-Heston rough volatility models with Zumbach effect. Quantitative Finance, 1–13.
  * Dung (2011) Dung, N.T. (2011). Semimartingale approximation of fractional Brownian motion and its applications. Computers and Mathematics with Applications, 61, 1844–1854.
  * Durrett (2018) Durrett, R. (2018). Stochastic Calculus: A Practical Introduction. CRC press.
  * El Euch et al. (2018) El Euch, O., Fukasawa, M. and Rosenbaum, M. (2018). The microstructural foundations of leverage effect and rough volatility. Finance and Stochastics, 22, 241–280.
  * El Euch et al. (2019) El Euch, O., Gatheral, J. and Rosenbaum, M. (2019). Roughening Heston. Risk, 84–89.
  * El Euch and Rosenbaum (2019) El Euch, O. and Rosenbaum, M. (2019). The characteristic function of rough Heston models. Mathematical Finance, 29, 3–38.
  * Eriksson and Pistorius (2015) Eriksson, B. and Pistorius, M.R. (2015). American option valuation under continuous-time markov chains. Advances in Applied Probability, 47, 378–401.
  * Ethier (2009) Ethier, S.N. and Kurtz, T.G. (2009). Markov Processes: Characterization and Convergence. John Wiley and Sons.
  * Forde et al. (2021a) Forde, M., Smith, B. and Viitasaari, L. (2021). Rough volatility, CGMY jumps with a finite history and the Rough Heston model-small-time asymptotics in the $k\sqrt{t}$ regime. Quantitative Finance, 21, 541–563.
  * Forde et al. (2021b) Forde, M., Gerhold, S. and Smith, B. (2021). Small-time, large-time, and asymptotics for the Rough Heston model. Mathematical Finance, 31, 203–241.
  * Forde and Zhang (2017) Forde, M. and Zhang, H. (2017). Asymptotics for rough stochastic volatility models. SIAM Journal on Financial Mathematics, 8, 114–145.
  * Fukasawa et al. (2019) Fukasawa, M., Takabatake, T. and Westphal, R. (2019). Is volatility rough?. arXiv preprint arXiv:1905.04852.
  * Gatheral et al. (2018) Gatheral, J., Jaisson, T. and Rosenbaum, M. (2018). Volatility is rough. Quantitative Finance, 18, 933–949.
  * Goudenège et al. (2020) Goudenège, L., Molent, A. and Zanette, A. (2020). Machine learning for pricing American options in high-dimensional Markovian and non-Markovian models. Quantitative Finance, 20, 573–591.
  * Grasselli (2017) Grasselli, M. (2017). The 4/2 stochastic volatility model: a unified approach for the Heston and the 3/2 model. Mathematical Finance, 27, 1013–1034.
  * Guennoun et al. (2018) Guennoun, H., Jacquier, A., Roome, P. and Shi, F. (2018). Asymptotic behavior of the fractional Heston model. SIAM Journal on Financial Mathematics, 9, 1017–1045.
  * Hagan et al. (2002) Hagan, P.S., Kumar, D., Lesniewski, A.S. and Woodward, D.E. (2002). Managing smile risk. Wilmott Magazine, 1, 84–108.
  * Harms (2019) Harms, P. (2019). Strong convergence rates for Markovian representations of fractional Brownian motion. arXiv preprint arXiv:1902.01471.
  * Heston (1993) Heston, S.L. (1993). A closed-form solution for options with stochastic volatility with applications to bond and currency options. The Review of Financial Studies, 6, 327–343.
  * Horvath et al. (2020) Horvath, B., Jacquier, A. and Tankov, P. (2020). Volatility options in rough volatility models. SIAM Journal on Financial Mathematics, 11, 437–469.
  * Hull and White (1987) Hull, J. and White, A. (1987). The pricing of options on assets with stochastic volatilities. The Journal of Finance, 42, 281–300.
  * Jaisson and Rosenbaum (2016) Jaisson, T. and Rosenbaum, M. (2016). Rough fractional diffusions as scaling limits of nearly unstable heavy tailed Hawkes processes. The Annals of Applied Probability, 26, 2860–2882.
  * Kushner and Dupuis (2001) Kushner, H.J. and Dupuis, P.G. (2001). Numerical Methods for Stochastic Control Problems in Continuous Time. Springer Science and Business Media.
  * Lewis (2000) Lewis A. (2000), Option Valuation under Stochastic Volatility, Finance Press.
  * Lipton (2002) Lipton, A. (2002). The vol smile problem. Risk Magazine, 15, 61–65.
  * Livieri et al. (2018) Livieri, G., Mouti, S., Pallavicini, A. and Rosenbaum, M. (2018). Rough volatility: evidence from option prices. IISE transactions, 50, 767–776.
  * Ma and Wu (2021) Ma, J. and Wu, H. (2021). A fast algorithm for simulation of rough volatility models, Quantitative Finance, forthcoming, DOI:10.1080/14697688.2021.1970213
  * Ma et al. (2021) Ma, J., Yang W. and Cui Z. (2021). Convergence analysis for continuous-time Markov chain approximation of stochastic local volatility models: option pricing and Greeks. Journal of Computational and Applied Mathematics, forthcoming, DOI: 10.1016/j.cam.2021.113901
  * McCrickerd et al. (2018) McCrickerd, R. and Pakkanen, M.S. (2018). Turbocharging Monte Carlo pricing for the rough Bergomi model. Quantitative Finance, 18, 1877–1886.
  * Mijatović and Pistorius (2013) Mijatović, A. and Pistorius, M. (2013). Continuously monitored barrier options under Markov processes. Mathematical Finance, 23, 1–38.
  * Oksendal (2013) Oksendal, B. (2013). Stochastic differential equations: an introduction with applications. Springer Science and Business Media.
  * Prigent (2003) Prigent, J.L. (2003). Weak Convergence of Financial Markets. Springer, Berlin, Heidelberg.
  * Quecke (2007) Quecke, S. (2007). Efficient numerical methods for pricing American options under Lévy models. Diss. Verlag nicht ermittelbar.
  * Richard et al. (2021) Richard, A., Tan, X. and Yang, F. (2021). On the discrete-time simulation of the rough Heston model. arXiv preprint arXiv:2107.07835.
  * Song et al. (2016) Song, Y., Cai, N. and Kou, S. (2016). A unified framework for options pricing under regime switching models. Available at SSRN 3310365.
  * Stein and Stein (1991) Stein, E.M. and Stein, J.C. (1991). Stock price distributions with stochastic volatility: an analytic approach. The Review of Financial Studies, 4, 727–752.
  * Tavella and Randall (2000) Tavella, D. and Randall, C. (2000). Pricing Financial Instruments: The Finite Difference Method. Wiley, New York.
  * Van der Stoep et al. (2014) Van der Stoep, A.W., Grzelak, L.A. and Oosterlee, C.W. (2014). The Heston stochastic-local volatility model: efficient Monte Carlo simulation. International Journal of Theoretical and Applied Finance, 17, 1–30.
  * Veraar (2012) Veraar, M. (2012). The stochastic Fubini theorem revisited. Stochastics An International Journal of Probability and Stochastic Processes, 84, 543–551.
  * Xi et al. (2019) Xi, Y., Ding, K. and Ning, N. (2019). Simultaneous two-dimensional continuous-time Markov chain approximation of two-dimensional fully coupled Markov diffusion processes. Available at SSRN 3461115.